( PDF ) Métodos de elementos finitos e diferenãs finitas para problema de Helmholtz

Download PDF

ads:

Laborat´orio Nacional de Computa¸c˜ao Cient´ıﬁca

Programa de P´os Gradua¸c˜ao em Modelagem Computacional

M´etodos de Elementos Finitos e Diferen¸cas

Finitas para o Problema de Helmholtz

Por

Daniel Thomes Fernandes

PETR

OPOLIS, RJ - BRASIL

MARC¸ O DE 2009

ads:

Livros Grátis

http://www.livrosgratis.com.br

Milhares de livros grátis para download.

ETODOS DE ELEMENTOS FINITOS E DIFERENC¸ AS FINITAS

PARA O PROBLEMA DE HELMHOLTZ

Daniel Thomes Fernandes

TESE SUBMETIDA AO CORPO DOCENTE DO LABORAT

ORIO NACIONAL

DE COMPUTAC¸

AO CIENT

IFICA COMO PARTE DOS REQUISITOS NE-

CESS

ARIOS PARA A OBTENC¸

AO DO GRAU DE DOUTOR EM MODELA-

GEM COMPUTACIONAL

Aprovada por:

Abimael Fernando Dourado Loula

Fr´ed´eric Valentin, Ph.D

Alexandre Madureira, Ph.D

Eduardo Gomes Dutra do Carmo, D.Sc

Fernando Alves Rochinha, D.Sc

Saulo Pomponet Oliveira, Ph.D

Gustavo Benitez Alvarez, D.Sc

PETR

OPOLIS, RJ - BRASIL

MARC¸ O DE 2009

ads:

FERNANDES, DANIEL THOMES

F363m M´etodos de Elementos Finitos e Diferen¸cas Finitas para o Problema de

Helmholtz / Daniel Thomes Fernandes. Petrop´olis, RJ. : LNCC/MCT, 2009.

xi, 126 p. : il.24 cm

Orientador:Abimael Fernando Dourado Loula

Tese (D.Sc.) – Laborat´orio Nacional de Computa¸c˜ao Cient´ıﬁca –

LNCC/MCT, 2009.

1. M´etodo dos elementos ﬁnitos, 2. M´etodos de diﬁren¸cas ﬁnitas, 3.

Equa¸c˜ao de Helmholtz, 4. M´etodos estabilizados

I. LNCC/MCT II. T´ıtulo

CDD 518.15

A Seraphina Thomes Fernandes, minha m˜ae,

e ao Daniel Fernandes, meu pai,

com muito orgulho.

Agradecimentos

Ao Prof Abimael Loula por me orientar, pela paciˆencia, pela conﬁan¸ca, pelo

incentivo e pelos ensinamentos sobre meu trabalho e sobre a vida.

Aos meus pais por terem me permitido ter tantas oportunidades que eles n˜ao

tiveram.

Aos colegas do LNCC pelo ambiente agrad´aval de trabalho.

A FAPERJ e `a FACC pelo suporte ﬁnanceiro.

Resumo da Tese apresentada ao LNCC/MCT como parte dos requisitos ne-

cess´arios para a obten¸c˜ao do grau de Doutor em Ciˆencias (D.Sc.)

ETODOS DE ELEMENTOS FINITOS E DIFERENC¸ AS FINITAS

PARA O PROBLEMA DE HELMHOLTZ

Daniel Thomes Fernandes

Mar¸co, 2009

Orientador: Abimael Fernando Dourado Loula

E bem sabido que m´etodos cl´assicos de elementos ﬁnitos e diferen¸cas ﬁnitas

para o problema de Helmholtz apresentam efeito de polui¸c˜ao, que pode deterio-

rar seriamente a qualidade da solu¸c˜ao aproximada. Controlar o efeito de polui¸c˜ao

´e especialmente dif´ıcil quando s˜ao utilizadas malhas n˜ao uniformes. Para malhas

uniformes com elementos quadrados s˜ao conhecidos m´etodos (p. e. o QSFEM, pro-

posto por Babuˇska et al) que minimizam a polui¸c˜ao. Neste trabalho apresentamos

inicialmente dois m´etodos de elementos ﬁnitos de Petrov-Galerkin com formula¸c˜ao

relativamente simples, o RPPG e o QSPG, ambos com razo´avel robustez para cer-

tos tipos de distor¸c˜oes dos elementos. O QSPG apresenta ainda polui¸c˜ao m´ınima

para elementos quadrados. Em seguida ´e formulado o QOFD, um m´etodo de di-

feren¸cas ﬁnitas aplic´avel a malhas n˜ao estruturadas. O QOFD apresenta grande

robustez em rela¸c˜ao a distor¸c˜oes, mas requer trabalho extra para tratar problemas

n˜ao homogˆeneos ou condi¸c˜oes de contorno n˜ao essenciais. Finalmente ´e apresen-

tado um novo m´etodo de elementos ﬁnitos de Petrov-Galerkin, o QOPG, que ´e

formulado aplicando a mesma t´ecnica usada para obter a estabiliza¸c˜ao do QOFD,

obtendo assim a mesma robustez em rela¸c˜ao a distor¸c˜oes da malha, com a van-

tagem de ser um m´etodo variacionalmente consistente. Resultados num´ericos s˜ao

apresentados ilustrando o comportamento de todos os m´etodos desenvolvidos em

compara¸c˜ao com os m´etodos de Galerkin, GLS e QSFEM.

vii

Abstract of Thesis presented to LNCC/MCT as a partial fulﬁllment of the

requirements for the degree of Doctor of Sciences (D.Sc.)

FINITE ELEMENTS AND FINITE DIFFERENCE METHODS FOR

THE HELMHOLTZ EQUATION

Daniel Thomes Fernandes

March, 2009

Advisor: Abimael Fernando Dourado Loula

It is well known that classical ﬁnite elements or ﬁnite diﬀerence methods

for Helmholtz problem present pollution eﬀects that can severely deteriorate the

quality of the approximate solution. To control pollution eﬀects is especially dif-

ﬁcult on non uniform meshes. For uniform meshes of square elements pollution

eﬀects can be minimized with the Quasi Stabilized Finite Element Method (QS-

FEM) proposed by Babusˇska el al, for example. In the present work we initi-

ally present two relatively simple Petrov-Galerkin ﬁnite element methods, referred

here as RPPG (Reduced Pollution Petrov-Galerkin) and QSPG (Quasi Stabili-

zed Petrov-Galerkin), with reasonable robustness to some type of mesh distortion.

The QSPG also shows minimal pollution, identical to QSFEM, for uniform meshes

with square elements. Next we formulate the QOFD (Quasi Stabilized Finite Dif-

ference) method, a ﬁnite diﬀerence method for unstructured meshes. The QOFD

shows great robustness relative to element distortion, but requires extra work to

consider non-essential boundary conditions and source terms. Finally we present a

Quasi Optimal Petrov-Galerkin (QOPG) ﬁnite element method. To formulate the

QOPG we use the same approach introduced for the QOFD, leading to the same

accuracy and robustness on distorted meshes, but constructed based on consistent

variational formulation. Numerical results are presented illustrating the behavior

of all methods developed compared to Galerkin, GLS and QSFEM.

viii

Sum´ario

1 Introdu¸c˜ao 1

2 O Problema de Helmholtz 7

2.1 Ondas ac´usticas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.2 Problemas modelos e formula¸c˜oes variacionais . . . . . . . . . . . . 9

2.2.1 Condi¸c˜oes de Robin . . . . . . . . . . . . . . . . . . . . . . . 9

2.2.2 Condi¸c˜oes de Dirichlet . . . . . . . . . . . . . . . . . . . . . 10

2.2.3 Condi¸c˜oes mistas . . . . . . . . . . . . . . . . . . . . . . . . 11

2.3 Solu¸c˜oes particulares do problema homogˆeneo . . . . . . . . . . . . 12

2.3.1 Ondas planas . . . . . . . . . . . . . . . . . . . . . . . . . . 12

2.3.2 Ondas evanescentes . . . . . . . . . . . . . . . . . . . . . . . 13

2.3.3 Ondas n˜ao direcionais . . . . . . . . . . . . . . . . . . . . . 13

2.4 Observa¸c˜ao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3 Solu¸c˜ao num´erica por elementos ﬁnitos e diferen¸cas ﬁnitas 17

3.1 M´etodo de Galerkin . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3.2 Polui¸c˜ao . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.3 N´umero de onda discreto . . . . . . . . . . . . . . . . . . . . . . . . 20

3.3.1 An´alise unidimensional . . . . . . . . . . . . . . . . . . . . . 21

3.3.2 An´alise bidimensional . . . . . . . . . . . . . . . . . . . . . . 24

3.4 Ressonˆancia Num´erica . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.5 Aproxima¸c˜ao por diferen¸cas ﬁnitas . . . . . . . . . . . . . . . . . . 28

3.6 Galerkin m´ınimos quadrados (GLS) . . . . . . . . . . . . . . . . . . 32

3.7 M´etodos com polui¸c˜ao m´ınima . . . . . . . . . . . . . . . . . . . . . 33

3.7.1 Elementos ﬁnitos generalizados: QSFEM . . . . . . . . . . . 34

3.7.2 Galerkin, volumes ﬁnitos e m´ınimos quadrados . . . . . . . . 37

3.7.3 Galerkin e m´ınimos quadrados ponderados . . . . . . . . . . 40

3.8 Observa¸c˜oes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

4 M´etodos de Petrov-Galerkin 44

4.1 M´etodo com polui¸c˜ao reduzida (RPPG) . . . . . . . . . . . . . . . . 44

4.2 M´etodo quase estabilizado (QSPG) . . . . . . . . . . . . . . . . . . 49

4.2.1 Conformidade . . . . . . . . . . . . . . . . . . . . . . . . . . 51

4.2.2 Parˆametros de estabiliza¸c˜ao . . . . . . . . . . . . . . . . . . 54

4.3 Implementa¸c˜ao computacional . . . . . . . . . . . . . . . . . . . . . 58

4.4 Resultados num´ericos . . . . . . . . . . . . . . . . . . . . . . . . . . 59

4.4.1 Elementos quadrados . . . . . . . . . . . . . . . . . . . . . . 59

4.4.2 Elementos retangulares . . . . . . . . . . . . . . . . . . . . . 60

4.4.3 Elementos distorcidos aleatoriamente . . . . . . . . . . . . . 62

4.4.4 Ondas n˜ao direcionais . . . . . . . . . . . . . . . . . . . . . 64

4.4.5 Problema n˜ao homogˆeneo . . . . . . . . . . . . . . . . . . . 66

4.4.6 Ondas evanescentes . . . . . . . . . . . . . . . . . . . . . . . 67

4.4.7 Problema de Dirichlet - ressonˆancia num´erica . . . . . . . . 69

4.5 Observa¸c˜oes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

5 M´etodos de diferen¸cas ﬁnitas para malhas n˜ao estruturadas 72

5.1 Nota¸c˜oes e deﬁni¸c˜oes . . . . . . . . . . . . . . . . . . . . . . . . . . 73

5.2 Truncamento nulo em algumas dire¸c˜oes . . . . . . . . . . . . . . . . 76

5.3 Minimiza¸c˜ao do erro de truncamento . . . . . . . . . . . . . . . . . 79

5.4 Caso unidimensional . . . . . . . . . . . . . . . . . . . . . . . . . . 83

5.5 Caso bidimensional . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

5.5.1 Malhas n˜ao estruturadas . . . . . . . . . . . . . . . . . . . . 84

5.5.2 Restri¸c˜ao a malhas uniformes com elementos quadrados . . . 85

5.6 Implementa¸c˜ao computacional . . . . . . . . . . . . . . . . . . . . . 86

5.7 Resultados num´ericos . . . . . . . . . . . . . . . . . . . . . . . . . . 88

5.7.1 Elementos quadrados . . . . . . . . . . . . . . . . . . . . . . 88

5.7.2 Elementos distorcidos . . . . . . . . . . . . . . . . . . . . . . 90

5.7.3 Malhas n˜ao estruturadas . . . . . . . . . . . . . . . . . . . . 92

5.7.4 Problema de Dirichlet . . . . . . . . . . . . . . . . . . . . . 94

5.8 Observa¸c˜oes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

6 M´etodos de Petrov-Galerkin com fun¸c˜oes peso quase ´otimas 97

6.1 Macroelementos e fun¸c˜oes bolha . . . . . . . . . . . . . . . . . . . . 98

6.2 M´etodo quase estabilizado . . . . . . . . . . . . . . . . . . . . . . . 99

6.3 Fun¸c˜oes peso quase ´otimas . . . . . . . . . . . . . . . . . . . . . . . 100

6.4 Condi¸c˜oes de contorno . . . . . . . . . . . . . . . . . . . . . . . . . 104

6.5 Implementa¸c˜ao computacional . . . . . . . . . . . . . . . . . . . . . 105

6.6 Resultados num´ericos . . . . . . . . . . . . . . . . . . . . . . . . . . 105

6.6.1 Malhas uniformes . . . . . . . . . . . . . . . . . . . . . . . . 105

6.6.2 Malhas n˜ao uniformes. Perturba¸c˜oes aleat´orias . . . . . . . . 108

6.6.3 Malhas n˜ao uniformes. Transforma¸c˜ao biquadr´atica . . . . . 111

6.6.4 Problema de Helmholtz n˜ao homogˆeneo . . . . . . . . . . . . 113

6.6.5 Ondas n˜ao direcionais . . . . . . . . . . . . . . . . . . . . . 115

7 Conclus˜oes e futuros desenvolvimentos 118

Referˆencias Bibliogr´aﬁcas 123

Cap´ıtulo 1

Introdu¸c˜ao

A equa¸c˜ao de Helmholtz tem importantes aplica¸c˜oes em problemas lineares

de propaga¸c˜ao de ondas harmˆonicas, como por exemplo ondas ac´usticas, ondas

el´asticas e em eletromagnetismo (Ihlenburg (1998)). Muitas vezes h´a interesse

pr´atico em obter solu¸c˜oes aproximadas tanto para problemas diretos ou inversos

envolvendo essa equa¸c˜ao.

Uma caracter´ıstica do problema direto de particular relevˆancia para a apli-

ca¸c˜ao de m´etodos de elementos ﬁnitos ou diferen¸cas ﬁnitas ´e a natureza oscilat´oria

das solu¸c˜oes da equa¸c˜ao de Helmholtz. Na aplica¸c˜ao desses m´etodos ´e necess´ario

ajustar a distˆancia m´axima h entre dois pontos nodais ao n´umero de onda k, de

modo que exista um n´umero m´ınimo desses pontos dentro de cada comprimento

de onda, permitindo assim que a solu¸c˜ao aproximada possa capturar as oscila¸c˜oes

da solu¸c˜ao exata. Estudando as propriedades do espa¸co de elementos ﬁnitos ﬁca

expl´ıcito que de fato, para diferentes valores de k, uma regra do tipo “kh constante”

´e suﬁciente para controlar o erro de interpola¸c˜ao.

No entanto, no caso dos m´etodos de elementos ﬁnitos baseados na formula¸c˜ao

cl´assica de Galerkin, os testes computacionais e a an´alise num´erica revelam que

manter kh constante n˜ao ´e suﬁciente quando k ´e grande. Apesar de o m´etodo

de Galerkin ter convergˆencia em rela¸c˜ao a h assintoticamente ´otima na norma

do espa¸co H

, quando se trata de valores altos de k o reﬁnamento da malha ne-

cess´ario para iniciar essa convergˆencia ´otima ´e impratic´avel. Esse comportamento

n˜ao robusto em rela¸c˜ao a k ´e conhecido como efeito de polui¸c˜ao e est´a relacionado

`a diferen¸ca entre o n´umero de onda k da solu¸c˜ao exata e o n´umero de onda da

solu¸c˜ao aproximada k

. Diﬁculdade da mesma natureza ocorre com m´etodos de di-

feren¸cas ﬁnitas. Valores altos de k s˜ao comuns em aplica¸c˜oes, tornando importante

a quest˜ao de como contornar o efeito de polui¸c˜ao.

Em problemas unidimensionais, o m´etodo GLS (Galerkin-M´ınimos Quadra-

dos) (Harari e Hughes (1991)) elimina a diferen¸ca k −k

, o que torna esse m´etodo

livre de polui¸c˜ao. Em duas ou mais dimens˜oes sabemos (Babuˇska e Sauter (1997))

que ´e imposs´ıvel eliminar totalmente o efeito de polui¸c˜ao. Por exemplo, o m´etodo

GLS permite fazer com que o n´umero de onda discreto coincida com o n´umero

de onda exato para ondas planas em uma dire¸c˜ao escolhida (Thompson e Pinsky

(1995)), mas em geral o erro de polui¸c˜ao continua sendo da mesma ordem que no

m´etodo de Galerkin. O melhor que se pode ter s˜ao os chamados m´etodos com

polui¸c˜ao m´ınima, nos quais o erro de polui¸c˜ao ´e de ordem bastante reduzida. O

m´etodo quase estabilizado (QSFEM) ´e descrito por Babuˇska et al. (1995) como

um m´etodo de elementos ﬁnitos generalizado com polui¸c˜ao m´ınima. Nesse m´etodo

os coeﬁcientes que deﬁnem o operador de Helmholtz aproximado s˜ao escolhidos de

modo que a diferen¸ca k −k

seja minimizada em certo sentido. Conforme proposto

por Babuˇska et al. (1995), o QSFEM ´e restrito a malhas uniformes com elementos

quadrados e, como n˜ao ´e constru´ıdo a partir de uma formula¸c˜ao variacional, o

QSFEM pode ser mais bem descrito como um m´etodo de diferen¸cas ﬁnitas.

M´etodos variacionalmente consistentes e com propriedades de estabiliza¸c˜ao

equivalentes `as do QSFEM foram obtidos posteriormente de v´arias formas, in-

cluindo o m´etodo baseado em res´ıduos proposto por Oberai e Pinsky (2000),

m´etodos de elementos ﬁnitos descont´ınuos (Alvarez et al. (2006)) e (Loula et al.

(2007)). Dois m´etodos tamb´em baseados em formula¸c˜oes variacionais residuais de

Galerkin ser˜ao apresentados neste trabalho, um combinando res´ıduos de Galer-

kin, volumes ﬁnitos e m´ınimos quadrados (Se¸c˜ao 3.7.2) e outro usando res´ıduos

de Galerkin e de m´ınimos quadrados ponderados (Se¸c˜ao 3.7.3). Esses m´etodos

tˆem propriedades de estabiliza¸c˜ao equivalentes `as do m´etodo QSFEM para malhas

uniformes com elementos quadrados, permenecendo em aberto a estabiliza¸c˜ao com

malhas n˜ao uniformes.

Todos os m´etodos citados acima possuem parˆametros livres que s˜ao ajustados

para malhas uniformes com elementos quadrados de lado h. Em tais malhas ´e

poss´ıvel, atrav´es da an´alise de dispers˜ao, obter a rela¸c˜ao entre o n´umero de onda

discreto e o n´umero de onda exato e com isso minimizar a diferen¸ca entre os dois

escolhendo parˆametros apropriados. Para malhas mais gerais a an´alise de dispers˜ao

n˜ao se aplica diretamente, sendo comum na pr´atica calcular um comprimento m´edio

h para cada elemento e ent˜ao adotar parˆametros de estabiliza¸c˜ao calculados como

se os elementos fossem quadrados. Entretanto, os parˆametros de estabiliza¸c˜ao

assim obtidos est˜ao longe dos valores ´otimos quando se trata de malhas distorcidas.

Dada a grande sensibilidade da solu¸c˜ao aproximada a varia¸c˜oes no n´umero de onda

discreto, essas estrat´egias baseadas em parˆametros ´otimos obtidos para malhas

uniformes n˜ao funcionam muito bem.

Podemos aﬁrmar que n˜ao existem m´etodos de elementos ﬁnitos lagrangianos

lineares ou bilineares precisos e robustos para an´alise num´erica do problema de

Helmholtz em malhas n˜ao uniformes e n˜ao estruturadas, para n´umeros de onda

elevados. De acordo com Zienkiewicz esse ´e um problema em aberto (Zienkiewicz

(2000)). Alternativamente, foram desenvolvidos m´etodos h´ıbridos usando fun¸c˜oes

n˜ao polinomiais, especialmente ondas planas (Farhat et al. (2003)). Ainda visando

estabilidade, tamb´em foram propostos m´etodos de elementos ﬁnitos com bolhas n˜ao

polinomiais (Harari e Gosteev (2007)) incluindo o RBF ou Residual Free Bubbles

(Franca et al. (1997)).

Outra quest˜ao igualmente relevante para a solu¸c˜ao num´erica do problema de

Helmholtz ´e a ressonˆancia. Dependendo das condi¸c˜oes de contorno e do dom´ınio,

podem existir n´umeros de onda para os quais o problema original admite inﬁnitas

solu¸c˜oes. Como o n´umero de onda discreto ´e diferente do original, ´e poss´ıvel ocorrer

ressonˆancia na solu¸c˜ao num´erica que n˜ao existe na solu¸c˜ao exata. Essa ressonˆancia

num´erica pode deteriorar seriamente a qualidade da solu¸c˜ao aproximada e ´e mais

grave quanto maior for o n´umero de onda. Esse problema tamb´em ´e amenizado

quando o erro no n´umero de onda ´e reduzido.

Objetivo da Tese

De forma simples e direta, o objetivo desta tese ´e formular m´etodos de ele-

mentos ﬁnitos e de diferen¸cas ﬁnitas para o problema de Helmholtz que sejam

quase est´aveis, no sentido deﬁnido por Babuˇska, isto ´e, capazes de minimizar ou

pelo menos reduzir signiﬁcativamente o erro de polui¸c˜ao t´ıpico destas aproxima¸c˜oes

e razoavelmente robustos quando aplicados a malhas mais gerais do que malhas

uniformes, incluindo malhas n˜ao estruturadas.

Principais resultados

No Cap´ıtulo 4 s˜ao propostos dois m´etodos de elementos ﬁnitos, formulados

como m´etodos de Petrov-Galerkin. O espa¸co U

onde a solu¸c˜ao aproximada ´e pro-

curada ´e o espa¸co com elementos lagrangianos C

bilineares usual. No primeiro

m´etodo, denominado RPPG (Reduced Pollution Petrov-Galerkin), o espa¸co V

´e

constru´ıdo de forma bastante simples e direta, sem parˆametros dependentes da geo-

metria dos elementos ou do n´umero de onda. Apesar da simplicidade observam-se

resultados muito melhores em rela¸c˜ao ao m´etodo de Galerkin, por exemplo, bem

como em rela¸c˜ao ao GLS. No segundo m´etodo (Quasi Stabilized Petrov-Galerkin,

ou QSPG) s˜ao introduzidos parˆametros que dependem diretamente dos compri-

mentos de todos os lados dos elementos e do n´umero de onda, aumentando um

pouco a complexidade mas resultando em polui¸c˜ao m´ınima no caso de elementos

quadrados e uma maior estabilidade em elementos mais gerais. Nos dois casos as

fun¸c˜oes base nodais do espa¸co V

s˜ao constru´ıdas para ter os mesmos suportes das

fun¸c˜oes base correspondentes no espa¸co U

, garantindo assim que as matrizes glo-

bais resultantes desses m´etodos tenham sempre o mesmo perﬁl de esparsidade que

a matriz correspondente do m´etodo de Galerkin. O RPPG e o QSPG mostram-se

especialmente eﬁcientes para elementos alongados. Para distor¸c˜oes mais gerais eles

n˜ao s˜ao t˜ao rubustos, mas mantˆem um desempenho melhor do que o m´etodos de

Galerkin ou GLS.

No Cap´ıtulo 5 ´e proposto um m´etodo de diferen¸cas ﬁnitas aplic´avel a malhas

n˜ao uniformes e n˜ao estruturadas, para dom´ınios bidimensionais ou tridimensio-

nais. Este m´etodo, que denominamos QOFD (Quase Optimal Finite Diﬀerence

method), ´e obtido atrav´es de um novo crit´erio usado na formula¸c˜ao do operador

discreto: a minimiza¸c˜ao do res´ıduo do operador aplicado a ondas planas em todas

as dire¸c˜oes poss´ıveis. Nesse crit´erio toda a geometria dos elementos ´e levada em

conta naturalmente. Quando todos os elementos s˜ao quadrados iguais, o erro no

n´umero de onda ´e da mesma ordem que nos m´etodos com polui¸c˜ao m´ınima. Expe-

rimentos num´ericos indicam que para malhas distorcidas os ganhos de estabilidade

s˜ao muito signiﬁcativos.

No Cap´ıtulo 6 o crit´erio usado no m´etodo QOFD ´e aplicado na contru¸c˜ao do

espa¸co das fun¸c˜oes peso para um m´etodo de elementos ﬁnitos de Petrov-Galerkin.

O objetivo dessa formula¸c˜ao ´e aliar os bons resultados do m´etodo QOFD para

malhas n˜ao estruturadas `a facilidade dos m´etodos variacionalmente consistentes

em tratar termos de fonte e condi¸c˜oes de contorno mais gerais. Trata-se de uma

formulac˜ao variacional conforme de Petrov-Galerkin nos espa¸cos de elementos ﬁni-

tos lagragianos usuais de classe C

, lineares ou bilineares, associados a espa¸cos de

fun¸c˜oes peso tamb´em lagrangianas geradas por combina¸c˜oes lineares de bolhas po-

linomiais deﬁnidas em macroelementos. Essa nova formula¸c˜ao de Petrov-Galerkin

denominamos QOPG (Quasi Optimal Petrov-Galerkin), em analogia a QOFD.

A an´alise num´erica de m´etodos de elementos ﬁnitos para o problema de

Helmholtz ´e, em geral, uma quest˜ao em aberto. Existem an´alises rigorosas para

problemas unidimensionais (Ihlenburg e Babuˇska (1995a), Ihlenburg e Babuˇska

(1997)), e expectativas de que os resultados destas an´alises podem ser aplica-

dos a problemas bi ou tri dimensionais. Por outro lado, apoiados em in´umeros

exemplos e contra-exemplos que ilustram as diﬁculdades inerentes ao problema

(polui¸c˜ao e ressonˆancia num´erica, por exemplo), existe um grande n´umero de es-

tudos num´ericos de convergˆencia visando caracterizar propriedades de estabilidade

e precis˜ao de novas formula¸c˜oes que s˜ao frequentemente propostas para superar

estas diﬁculdades.

Dentre todas as formula¸c˜oes aqui estudadas, no contexto de m´etodos de ele-

mentos ﬁnitos lagrangianos lineares ou bilineares, a formula¸c˜ao de Petrov-Galerkin

do Cap´ıtulo 6 (o QOPG) ´e a que apresenta melhores propriedades de estabilidade,

precis˜ao e robustez a distor¸c˜oes de malhas em um amplo conjunto de testes de

convergˆencia.

Cap´ıtulo 2

O Problema de Helmholtz

Neste cap´ıtulo a equa¸c˜ao de Helmholtz com condi¸c˜oes de contorno mais co-

muns ´e apresentada brevemente, inicialmente como um problema de ac´ustica e em

seguida na forma de problemas modelos que ser˜ao tratados neste trabalho.

2.1 Ondas ac´usticas

Uma das aplica¸c˜oes da equa¸c˜ao de Helmholtz mais frequente e relevante ´e

o estudo de propaga¸c˜ao de ondas sonoras. Nesse tipo de problema as vari´aveis

envolvidas (press˜ao, velocidade, densidade) s˜ao pequenas varia¸c˜oes de um estado

est´atico, o que torna poss´ıvel linearizar e simpliﬁcar as equa¸c˜oes gerais que gover-

nam o comportamento dessas grandezas e chegar `a equa¸c˜ao da onda

∆φ −

∂

∂t

= 0, (2.1)

onde φ pode ser tanto o campo de press˜ao quanto um potencial de velocidade e c

´e a velocidade do som no meio.

Para ondas harmˆonicas no tempo da forma

φ(x, t) = Re



u(x)e

−iωt

 

i =

√

−1



(2.2)

com frequˆencia ω > 0 deduzimos que o campo complexo u, que representa a

componete espacial do campo φ, satisfaz a equa¸c˜ao de Helmholtz

− ∆u − k

u = 0, (2.3)

com k

= ω

Em problemas de espalhamento ac´ustico uma onda incidente conhecida u

´e espalhada por um obst´aculo D ⊂ R

resultando em uma onda espalhada u

Nesse caso o campo total u = u

+ u

satisfaz a equa¸c˜ao de Helmholtz em R

 D.

Quando D ´e um obst´aculo sound-hard, a velocidade normal se anula no contorno,

ou seja, a onda total u satisfaz condi¸c˜oes de contorno de Neumman

∂u

∂n

= 0 em ∂D. (2.4)

Quando o obst´aculo ´e do tipo sound-soft o excesso de press˜ao se anula em ∂D,

resultando em condi¸c˜oes de Dirichlet

u = 0 em ∂D. (2.5)

Uma situa¸c˜ao mais geral ´e quando a velocidade normal ´e proporcional ao excesso

de press˜ao no contorno, nesse caso temos condi¸c˜oes de Robin

∂u

∂n

+ iλu = 0 em ∂D, (2.6)

para uma constante λ > 0. Finalmente, para garantir a unicidade da solu¸c˜ao, ´e

introduzida a condi¸c˜ao de radia¸c˜ao de Sommerfeld

lim

r→∞



∂u

∂r

− iku



= 0 (2.7)

onde r = |x|.

No contexto de elementos ﬁnitos, problemas exteriores como o descrito acima

s˜ao tratados introduzindo condi¸c˜oes de contorno absorventes aproximadas em um

contorno artiﬁcial em torno do obst´aculo. Esse tipo de problema n˜ao ser´a tratado

neste trabalho, visto que o nosso foco ´e na discretiza¸c˜ao do operador de Helmholtz,

que, como ser´a visto ao longo deste trabalho, apresenta por si s´o diﬁculdades

consider´aveis.

Assim, nos limitaremos apenas a problemas interiores com condi¸c˜oes de Di-

richlet e/ou Robin.

2.2 Problemas modelos e formula¸c˜oes variacionais

Os m´etodos que ser˜ao apresentados neste trabalho ser˜ao testados para o

problema de Helmholtz com trˆes tipos de condi¸c˜oes de contorno. Em todos os

casos consideramos um dom´ınio aberto limitado Ω ⊂ R

com contorno ∂Ω Lipschitz

cont´ınuo.

Os m´etodos de elementos ﬁnitos s˜ao formulados a partir da forma variacional

do problema tratado. Para os casos considerados neste trabalho temos sempre a

forma variacional abstrata:

• Encontrar o campo u ∈ U tal que

a(u, v) = f(v) para todo v ∈ V, (2.8)

onde a forma sesquilinear a(u, v), o funcional antilinear f(v) e os conjuntos U e

V variam de acordo com as condi¸c˜oes de contorno. A seguir ser˜ao apresentados

resumidamente algumas situa¸c˜oes consideradas neste trabalho.

2.2.1 Condi¸c˜oes de Robin

Na forma cl´assica procuramos u satisfazendo

−∆u − k

u = f em Ω, (2.9)

∂u

∂n

+ iku = g em ∂Ω. (2.10)

A forma variacional (2.8) associada ao problema (2.9)-(2.10) ´e obtida com

a(u, v) =



Ω

∇u∇v − k

uv dx + ik



∂Ω

uv ds, (2.11)

f(v) =



Ω

fv dx +



∂Ω

gv ds, (2.12)

U = V = H

(Ω), (2.13)

onde v denota o conjugado complexo de v.

Para f ∈ H

−1

(Ω) e g ∈ H

−

(∂Ω), o problema variacional acima admite uma

´unica solu¸c˜ao u (Melenk (1995), Proposi¸c˜ao 8.1.3) que satisfaz a estimativa

∇u

(Ω)

+ ku

(Ω)

≤ C(Ω, k)



f

−1

(Ω)

+ g

−

(Ω)



. (2.14)

Uma caracter´ıstica importante desta formula¸c˜ao merece destaque. A dis-

sipa¸c˜ao gerada pelas condi¸c˜oes de contorno de Robin previne a ressonˆancia, evi-

tando assim instabilidades devidas a ressonˆancia num´erica. O mesmo n˜ao ocorre

com condi¸c˜oes de contorno de Dirichlet ou Neumann.

2.2.2 Condi¸c˜oes de Dirichlet

O problema de Helmholtz com condi¸c˜oes de contorno de Dirichlet

−∆u − k

u = f em Ω, (2.15)

u = g em ∂Ω, (2.16)

´e equivalente a forma fraca (2.8) com:

a(u, v) =



Ω

∇u∇v − k

uv dx, (2.17)

f(v) =



Ω

fv dx, (2.18)

U =



u ∈ H

(Ω) : u = g em ∂Ω



, (2.19)

V =



v ∈ H

(Ω) : v = 0 em ∂Ω



. (2.20)

Um aspecto importante do problema de Dirichlet ´e que quando k

´e um

autovalor do operador −∆ com g = 0, n˜ao h´a unicidade da solu¸c˜ao, ocorrendo o

fenˆomeno conhecido como ressonˆancia. Por exemplo, se Ω for o retˆangulo (0, a) ×

(0, b), ent˜ao os autovalores s˜ao dados por

= π













, (2.21)

com autofun¸c˜oes associadas

(x, y) = sen



mπx



sen



nπx



, (2.22)

onde m, n ∈ {1, 2, 3, . . . }. Assim, se u ´e uma solu¸c˜ao do problema (2.15)-(2.16)

com k

igual a um dos λ

, ent˜ao toda fun¸c˜ao da forma u + αu

com qualquer

constante α tamb´em ´e uma solu¸c˜ao.

No contexto dos m´etodos num´ericos pode ocorrer de o n´umero de onda da

solu¸c˜ao aproximada e os autovalores do problema discreto serem bastante dife-

rentes do n´umero de onda exato e dos autovalores originais. Assim, mesmo que

k n˜ao cause ressonˆancia, pode acontecer de observarmos ressonˆancia na solu¸c˜ao

aproximada, problema conhecido como ressonˆancia num´erica. Esse assunto ser´a

retomado e exempliﬁcado no Cap´ıtulo 3.

2.2.3 Condi¸c˜oes mistas

Uma outra situa¸c˜ao que ser´a testada numericamente ´e o problema com

condi¸c˜oes mistas:

−∆u − k

u = f em Ω, (2.23)

u = g em Γ

, (2.24)

∂u

∂n

+ iku = r em Γ

, (2.25)

onde Γ

∪ Γ

= ∂Ω e Γ

∩ Γ

= ∅. Nesse caso a formula¸c˜ao variacional (2.8) ´e

obtida com

a(u, v) =



Ω

∇u∇

v − k

uv dx + ik



uv ds, (2.26)

f(v) =



Ω

fv dx +



rv ds, (2.27)

U =



u ∈ H

(Ω) : u = g em Γ



, (2.28)

V =



v ∈ H

(Ω) : v = 0 em Γ



. (2.29)

2.3 Solu¸c˜oes particulares do problema homogˆeneo

Apresentamos a seguir algumas solu¸c˜oes particulares que ser˜ao consideradas

neste trabalho, no caso em que f = 0.

2.3.1 Ondas planas

Dado um vetor unit´ario σ ∈ R

, a onda plana na dire¸c˜ao σ ´e dada por

w(x) = e

ikσ·x

. (2.30)

Considerando que

∆w = div(∇w)

= div(ikwσ)

= ik(w div σ + ∇w · σ)

= ik(0 + (ikwσ) · σ)

= i

σ · σw

= −k

w, (2.31)

conclu´ımos que tais ondas planas satisfazem o problema de Helmholtz homogˆeneo.

No caso bidimensional, a dire¸c˜ao pode ser escrita como σ = (cos θ, sen θ).

Substituindo na express˜ao (2.30) obtemos a onda plana na dire¸c˜ao θ:

w(x, y) = e

ik(x cos θ+y sen θ)

(2.32)

= cos



k(x cos θ + y sen θ)



+ i sen



k(x cos θ + y sen θ)



. (2.33)

2.3.2 Ondas evanescentes

Dados

α > k, (2.34)

β =

√

− k

, (2.35)

θ ∈ R, (2.36)

ent˜ao a fun¸c˜ao u : R

→ C deﬁnida por

u(x, y) = e

−β(x sen θ−y cos θ)

iα(x cos θ+y sen θ)

(2.37)

satisfaz a equa¸c˜ao de Helmholtz homogˆenea, como pode ser veriﬁcado diretamente

calculando as derivadas parciais.

Tais ondas tˆem comportamento oscilat´orio na dire¸c˜ao θ e exponencial na

dire¸c˜ao perpendicular a θ.

2.3.3 Ondas n˜ao direcionais

Para veriﬁcar o desempenho dos m´etodos num´ericos quando a solu¸c˜ao n˜ao ´e

uma onda com uma dire¸c˜ao bem deﬁnida, consideraremos solu¸c˜oes constru´ıdas da

forma que segue.

Em coordenadas polares (r, θ), a equa¸c˜ao de Hemholtz homogˆenea ´e

∂

∂r



∂u

∂r



∂

∂θ

+ k

u = 0. (2.38)

Buscando uma solu¸c˜ao por separa¸c˜ao de vari´aveis da forma

u(r, θ) = R(r)T (θ), (2.39)

com T peri´odica para assegurar regularidade, obtemos

T +

dθ

R + k

RT = 0. (2.40)

Multiplicando por r/(RT ), resulta



+ k





dθ



= 0. (2.41)

Considerando que um dos termos acima depende somente de r e o outro somente de

θ, ent˜ao a igualdade s´o pode ser satisfeita se ambos forem constantes. Resolvendo

dθ

= −n

(2.42)

para uma constante n, obtemos

T (θ) = A cos(nθ) + B sen(nθ). (2.43)

Para que u seja cont´ınua, T deve ser peri´odica com per´ıodo 2π, logo n deve ser

inteiro. J´a a solu¸c˜ao geral da equa¸c˜ao



+ k



− n

= 0, (2.44)

para r = 0, ´e dada por

R(r) = CJ

(kr) + DY

(kr), (2.45)

onde J

e Y

s˜ao as fun¸c˜oes de Bessel de ordem n, de primeiro e segundo tipo

respectivamente. As solu¸c˜oes que ser˜ao testadas s˜ao obtidas com as constantes

A = 1, (2.46)

B = 0, (2.47)

C = 1, (2.48)

D = i (2.49)

e com v´arios valores de n ∈ {0, 1, 2, ···}. Ou seja

u(r, θ) = cos(nθ)H

(1)

(kr), (2.50)

onde H

(1)

= J

+ iY

´e a fun¸c˜ao de Hankel de primeiro tipo e ordem n. Y

(r) ´e

singular em r = 0, portanto consideraremos apenas dom´ınios Ω tais que 0 /∈ Ω.

Fazendo as substitui¸c˜oes

r =



+ y

, (2.51)

cos θ =



+ y

(2.52)

e lembrando que cos(nθ) = T

(cos θ), onde T

´e o polinˆomio de Chebyshev de

primeiro tipo e ordem n, ent˜ao obtemos em coordenadas cartesianas:

u(x, y) = T





+ y



(1)





+ y



. (2.53)

2.4 Observa¸c˜ao

An´alise num´erica e estudos de convergˆencia de aproxima¸c˜oes por elemen-

tos ﬁnitos ou por diferen¸cas ﬁnitas mostram que as propriedades de estabilidade,

convergˆencia e precis˜ao destas aproxima¸c˜oes est˜ao intimamente relacionadas `as cor-

respondentes discretiza¸c˜oes do operador de Helmholtz. Assim, tem sido observado

que os efeitos de polui¸c˜ao e ressonˆancia num´erica s˜ao muito sens´ıveis a ordem des-

sas aproxima¸c˜oes. Por outro lado, a estabiliza¸c˜ao das aproxima¸c˜oes de um modo

geral e, em particular, das aproxima¸c˜oes de baixa ordem requer conhecimento de

solu¸c˜oes do problema homogˆeneo na determina¸c˜ao de parˆametros de estabiliza¸c˜ao

(como no GLS, por exemplo) ou na pr´opria constru¸c˜ao do operador de Helmholtz

discreto (como no QSFEM). Neste trabalho frequentemente ser˜ao utilizadas ondas

planas, solu¸c˜oes do problema de Helmholtz, na determina¸c˜ao de parˆametros de es-

tabiliza¸c˜ao em formula¸c˜oes de elementos ﬁnitos de Petrov-Galerkin e na constru¸c˜ao

de aproxima¸c˜oes por difere¸cas ﬁnitas.

Cap´ıtulo 3

Solu¸c˜ao num´erica por elementos ﬁnitos e

diferen¸cas ﬁnitas

Neste cap´ıtulo ser´a apresentada uma breve revis˜ao da literatura sobre m´e-

todos de elementos ﬁnitos e de diferen¸cas ﬁnitas para o problema de Helmholtz.

Inicialmente ser´a apresentado o m´etodo cl´assico de elementos ﬁnitos lagrangianos

lineares ou bilineares baseado na formula¸c˜ao de Galerkin, discutidas suas propri-

edades ressaltando as suas limita¸c˜oes. Alguns m´etodos de elementos ﬁnitos esta-

bilizados baseados em formula¸c˜oes residuais, como o GLS, ser˜ao tamb´em revistos

bem como m´etodos de elementos ﬁnitos generalizados, como o QSFEM que mais

se adapta a estrutura dos m´etodos de diferen¸cas ﬁnitas.

3.1 M´etodo de Galerkin

Seja M

= {Ω

, . . . , Ω

} uma parti¸c˜ao regular de Ω em Ne elementos ﬁnitos

n˜ao degenerados tais que

Ω

∩ Ω



= ∅ para e = e



e (3.1)

Ω ∪ ∂Ω =



e=1

(Ω

∪ ∂Ω

). (3.2)

Nos m´etodos de elementos ﬁnitos tratados neste trabalho ser˜ao considerados ape-

nas elementos lagrangianos C

lineares em 1D ou bilineares em 2D (elementos

quadrilaterais). O espa¸co de elementos ﬁnitos gerado por esses elementos ser´a in-

dicado por S

(Ω). As fun¸c˜oes base nodais de S

associadas aos pontos nodais x

ser˜ao denotadas por φ

Deﬁnindo apropriadamente conjuntos U

⊂ S

e V

⊂ S

, a solu¸c˜ao aproxi-

mada pelo m´etodo de Galerkin ´e obtida resolvendo o problema variacional:

• Encontrar u

∈ U

tal que

a(u

, v

) = f(v

) ∀ v

∈ V

, (3.3)

onde a(u, v) e f(v) s˜ao as formas deﬁnidas no Cap´ıtulo 2 para cada tipo de condi¸c˜ao

de contorno associada `a equa¸c˜ao de Helmholtz.

Para o m´etodo de Galerkin, em todos os casos considerados no Cap´ıtulo 2

temos que

= S

∩ V. (3.4)

No caso do conjunto U

, para condi¸c˜oes de Robin temos

= S

. (3.5)

Para o problema de Dirichlet, em geral n˜ao ´e poss´ıvel que as fun¸c˜oes do espa¸co

satisfa¸cam as condi¸c˜oes de contorno exatamente. Assim, para o problema de

Dirichlet temos que

= {u

∈ S

: u

= g

em ∂Ω}, (3.6)

onde g

´e a interpolante de g em ∂Ω.

3.2 Polui¸c˜ao

Para problemas el´ıpticos sim´etricos a aproxima¸c˜ao de Galerkin ´e a melhor

aproxima¸c˜ao na norma da energia. Mas para o problema de Helmholtz esse n˜ao

´e o caso. Para valores altos de k o operador de Helmholtz se torna indeﬁnido e

o m´etodo de Galerkin apresenta s´erios problemas, como instabilidade e polui¸c˜ao

(Ihlenburg e Babuˇska (1995a), Babuˇska et al. (1995), Babuˇska e Sauter (1997),

Ihlenburg (1998)). Uma an´alise completa do m´etodo de Galerkin para o caso

unidimensional ´e feita em (Ihlenburg e Babuˇska (1995a)) onde a seguinte estimativa

´e provada para o erro relativo na seminorma do H

|u − u

|u|

≤ C

kh + C

, kh < 1, (3.7)

com C

e C

independentes de k e h.

O primeiro termo na estimativa acima corresponde ao erro de interpola¸c˜ao,

que ´e da mesma ordem do erro da melhor aproxima¸c˜ao no espa¸co de elementos

ﬁnitos usado. O segundo corresponde `a polui¸c˜ao num´erica. Para diferentes va-

lores de k, manter kh constante ´e suﬁciente para manter o erro de interpola¸c˜ao

constante. Mas para manter o erro de polui¸c˜ao controlado ´e preciso manter k

suﬁcientemente pequeno (Douglas et al. (1993)). Ou seja, para valores crescentes

de k o m´etodo de Galerkin ´e cada vez menos capaz de bem aproveitar a capacidade

de aproxima¸c˜ao do espa¸co de elementos ﬁnitos, o que faz com que a aproxima¸c˜ao

de Galerkin se torne cada vez mais distante da melhor aproxima¸c˜ao.

A Figura 3.1 compara a convergˆencia da interpolante e da aproxima¸c˜ao de

Galerkin para um problema unidimensional com k = 100, com condi¸c˜oes de Di-

richlet em x = 0 e Robin em x = 1. Podemos notar que para h suﬁcientemente

pequeno o erro de u

se aproxima do erro da interpolante na norma do H

, mas

isso ocorre apenas para h < 1/1000, o que teria um custo exageradamente alto

para problemas bidimensionais, ou ainda maior para problemas tridimensionais.

Em (Ihlenburg e Babuˇska (1995b)) ´e provada a seguinte estimativa na norma

(0, 1):

u − u



u

≤ (C

+ C

k)k

. (3.8)

Como se pode observar por essa estimativa, na norma L

o efeito da polui¸c˜ao no

m´etodo de Galerkin n˜ao ´e eliminado com o reﬁnamento, como ocorre na norma H

(equa¸c˜ao (3.7)). Podemos observar isso na Figura 3.1, o gr´aﬁco da convergˆencia

em L

da aproxima¸c˜ao de Galerkin nunca se aproxima do gr´aﬁco da interpolante,

mantendo sempre um gap, que depende do n´umero de onda k.

Essa situa¸c˜ao torna-se ainda mais cr´ıtica com a poss´ıvel ocorrˆencia de res-

sonˆancias num´ericas, dependendo do tipo de condi¸c˜ao de contorno.

2 2.5 3 3.5 4

−log(h)

-2.5

-2

-1.5

-1

-0.5

log



|u−u

|u|



Interpolante

Galerkin

2 2.5 3 3.5 4

−log(h)

-5

-4

-3

-2

-1

log



u−u



u



Interpolante

Galerkin

Figura 3.1: Convergˆencia para o problema unidimensional, k = 100, condi¸c˜oes de

contorno de Dirichlet em x = 0 Robin em x = 1.

Para o caso bidimensional n˜ao existe uma an´alise num´erica completa para o

m´etodo de Galerkin. Mas ´e de se esperar que o efeito de polui¸c˜ao seja pelo menos

t˜ao cr´ıtico quando no caso unidimensional, o que de fato tem sido veriﬁcado nos

estudos num´ericos.

3.3 N´umero de onda discreto

Uma caracter´ıstica estreitamente relacionada ao efeito de polui¸c˜ao ´e a dife-

ren¸ca entre o n´umero de onda k da solu¸c˜ao exata e o n´umero de onda k

da solu¸c˜ao

aproximada. Para melhor ilustrar essa quest˜ao, consideramos inicialmente o caso

unidimensional. Em seguida o caso bidimensional ser´a tamb´em analisado.

3.3.1 An´alise unidimensional

Consideremos o problema unidimensional



+ k

u = 0 em (0, 1). (3.9)

Suponhamos que as condi¸c˜oes de contorno e o n´umero de onda k estejam devida-

mente deﬁnidos de modo que a solu¸c˜ao seja ´unica. A solu¸c˜ao geral desse problema

´e dada pela express˜ao

u(x) = c

ikx

+ c

−ikx

, (3.10)

onde c

e c

s˜ao determinados de acordo com as condi¸c˜oes de contorno.

Se considerarmos uma discretiza¸c˜ao uniforme do dom´ınio em elementos de

comprimento h, com pontos nodais x

= ih, ent˜ao as equa¸c˜oes associadas aos n´os

interiores tˆem todas a forma

i−1

+ 2Su

+ Ru

i+1

= 0, (3.11)

onde u

s˜ao os valores nodais da solu¸c˜ao aproximada. Procurando por uma solu¸c˜ao

da forma

= e

, (3.12)

ent˜ao obtemos, para os n´os interiores, as equa¸c˜oes

k(x

−h)

+ 2Se

+ Re

k(x

+h)

= 0. (3.13)

Dividindo a equa¸c˜ao anterior por e

obtemos

−i

+ 2S + Re

= 0, (3.14)

ou ainda,

2R cos(

kh) + 2S = 0, (3.15)

de onde obtemos duas solu¸c˜oes

k = ±k

(3.16)

onde k

´e o chamado n´umero de onda discreto:

arccos



−



. (3.17)

Assim, temos a forma geral da solu¸c˜ao aproximada

(x) = ˜c

+ ˜c

−ik

, (3.18)

com um n´umero de onda k

possivelmente diferente do n´umero de onda exato k.

No m´etodo de Galerkin temos

R = −1 −

(kh)

, (3.19)

S = 1 −

(kh)

. (3.20)

Nesse caso a equa¸c˜ao (3.17) tem solu¸c˜ao real desde que kh <

√

12. Sendo assim,

expandindo (3.17) em s´erie de Taylor em torno de kh = 0 obtemos

k − k

(kh)

+ O



(kh)



. (3.21)

A Figura 3.2 ilustra o efeito dessa diferen¸ca para crescentes valores do n´umero de

onda, com o erro de interpola¸c˜ao constante (kh = 1). Foram usadas condi¸c˜oes de

Dirichlet em x = 0 e Robin em x = 1. A diferen¸ca no n´umero de onda n˜ao causa

tanto impacto quando k = 10, mas quando k = 80, a aproxima¸c˜ao de Galerkin

chega a ﬁcar com a fase invertida para x pr´oximo de 1. Assim, nesse estudo ﬁca

claro que o problema causado pela diferen¸ca k

−k ´e mais grave para n´umeros de

ondas mais altos.

Parece razo´avel tomar como crit´erio na formula¸c˜ao de m´etodos estabilizados

a redu¸c˜ao ou minimiza¸c˜ao da diferen¸ca k − k

. De fato, para o problema unidi-

-1.5

-1

-0.5

0.5

1.5

0 0.2 0.4 0.6 0.8 1

Solucao u para k = 10.0

Re(u)

Galerkin, h = 1/10

-1.5

-1

-0.5

0.5

1.5

0 0.2 0.4 0.6 0.8 1

Solucao u para k = 40.0

Re(u)

Galerkin, h = 1/40

-1.5

-1

-0.5

0.5

1.5

0 0.2 0.4 0.6 0.8 1

Solucao u para k = 80.0

Re(u)

Galerkin, h = 1/80

Figura 3.2: Aproxima¸c˜ao de Galerkin comparada com a solu¸c˜ao exata para dife-

rentes valores de k, mantendo a resolu¸c˜ao da malha constante e portanto o erro de

interpola¸c˜ao controlado. Condi¸c˜oes de Dirichlet em x = 0 e Robin em x = 1 para

evitar ressonˆancia num´erica.

mensional, se um m´etodo resulta em um n´umero de onda discreto igual ao original,

ent˜ao ´e poss´ıvel discretizar as condi¸c˜oes de contorno e termo n˜ao homogˆeneo f de

modo que a polui¸c˜ao seja eliminada (Babuˇska e Sauter (1997)). Por exemplo, o

m´etodo GLS para problemas homogˆeneos tem essa propriedade no caso unidimen-

sional.

3.3.2 An´alise bidimensional

Considere uma discretiza¸c˜ao do dom´ınio consistindo em uma malha uniforme

com elementos quadrados. Por quest˜ao de simetria e invariˆancia em rela¸c˜ao a

transla¸c˜ao, as equa¸c˜oes associadas a um n´o interior qualquer (x

, y

) tem a forma

geral

− h, y

+ h) + A

, y

+ h) + A

+ h, y

+ h)

+ A

− h, y

) + A

, y

) + A

+ h, y

)

+ A

− h, y

− h) + A

, y

− h) + A

+ h, y

− h) = 0.

(3.22)

Se procuramos por solu¸c˜oes que sejam ondas planas discretas da forma

, y

) = e

cos θ+y

sen θ)

(3.23)

com um n´umero de onda discreto k

, ent˜ao, substituindo (3.23) em (3.22), encon-

tramos a seguinte equa¸c˜ao n˜ao linear que determina k

em fun¸c˜ao de k para cada

dire¸c˜ao θ:

+ 2A



cos(k

h cos θ) + cos(k

h sen θ)



+ 4A

cos(k

h cos θ) cos(k

h sen θ) = 0.

(3.24)

No m´etodo de Galerkin os coeﬁcientes A

, A

e A

s˜ao dados por

−

(kh)

, (3.25)

= −

−

(kh)

, (3.26)

= −

−

(kh)

. (3.27)

Substituindo as express˜oes acima em (3.24) e fazendo a expans˜ao

h =

∞



n=1

(kh)

, (3.28)

ent˜ao transformamos a equa¸c˜ao de dispers˜ao (3.24) em uma equa¸c˜ao da forma

F (kh) = 0. (3.29)

Expandindo F (kh) em em torno de kh = 0, ent˜ao a s´erie resultante tem que

ser identicamente nula em cada ordem de kh. Com isso obtemos equa¸c˜oes en-

volvendo os coeﬁcientes r

da expans˜ao (3.28) e determinamos tais r

resolvendo

essas equa¸c˜oes. Usamos o software Matematica (Wolfram Research (2005)) para

executar esses c´alculos e obter a rela¸c˜ao entre k e k

k − k

3 + cos 4θ

(kh)

+ O



(kh)



. (3.30)

A metodologia descrita acima ser´a usada adiante para obter essa rela¸c˜ao entre k e

para outras aproxima¸c˜oes.

Conforme Babuˇska e Sauter (1997), o valor m´aximo do m´odulo da diferen¸ca

acima para θ ∈ [0, 2π] pode ser visto como um indicador da qualidade de apro-

xima¸c˜ao de um m´etodo que gera um estˆencil de 9 pontos como o da equa¸c˜ao (3.22).

Isso ´e usado para formular o QSFEM, que ser´a visto na Se¸c˜ao 3.7.1.

3.4 Ressonˆancia Num´erica

Consideremos novamente o problema unidimensional



+ k

u = 0 em (0, 1), (3.31)

com condi¸c˜oes de contorno de Dirichlet

u(0) = a e u(1) = b. (3.32)

A solu¸c˜ao desse problema ´e dada por

u(x) =

a sen(k − kx) + b sen(kx)

sen k

. (3.33)

Quando sen(k) se aproxima de zero, a amplitude da solu¸c˜ao (3.33) cresce indeﬁ-

nidamente. Enquanto que para os valores de k tais que sen(k) = 0, o problema

(3.31)-(3.32) admite inﬁnitas solu¸c˜oes. Isso acontece pois para tais valores de k,

torna-se um autovalor do operador −∆ em 1D com condi¸c˜oes de Dirichlet ho-

mogˆeneas. Esses autovalores e autofun¸c˜oes s˜ao dados por

= n

, (3.34)

(x) = sen



x, (3.35)

para n ∈ {1, 2, . . . }. Essa coincidˆencia entre k

e algum λ

tem sido referida como

ressonˆancia.

O m´etodo de Galerkin com elementos lineares de comprimento h, para o

problema considerado, resulta em equa¸c˜oes associadas aos n´os interiores da forma

i−1

) + A

i+1

i−1

) + B

i+1

)) = 0, (3.36)

com os valores prescritos

(0) = a e u

(1) = b. (3.37)

A solu¸c˜ao desse problema discreto ´e dada explicitamente por

) =

a sen(k

− k

) + b sen(k

)

sen k

, (3.38)

onde k

´e o n´umero de onda discreto:

arccos



−

+ k

+ 2k



. (3.39)

Podemos notar nas equa¸c˜oes acima que, para o problema discreto, ocorre

um fenˆomeno similar `a ressonˆancia. Mas neste caso a amplitude da solu¸c˜ao tende

para inﬁnito quando sen





→ 0, o que pode ocorrer mesmo que a solu¸c˜ao exata

esteja longe de uma ressonˆancia. Os n´umeros de onda discretos que fazem com

que ocorra essa ressonˆancia na solu¸c˜ao aproximada s˜ao dados por

= nπ, para n ∈ {1, 2, . . . }. (3.40)

Substituindo em (3.39), encontramos que os n´umeros de onda k que fazem com

que o problema discreto tenha ressonˆancia tˆem a forma

−A

− 2A

cos(hnπ)

+ 2B

cos(hnπ)

, para n ∈ {1, 2, . . . }. (3.41)

N˜ao por acaso os n´umeros acima s˜ao os autovalores do problema de autovalor

generalizado



, v



) = −λ(u

, v

) ∀v

∈ V

, (3.42)

associado a formula¸c˜ao variacional do problema de Dirichlet (3.31)-(3.32).

Assim, os n´umeros de onda que causam ressonˆancia no problema discreto

podem ser diferentes dos que causam ressonˆancia no problema original e isso se

deve ao fato de esses dois problemas terem autovalores diferentes. Se tiv´essemos

= k ent˜ao esse tipo de problema n˜ao ocorreria.

Na pr´atica, as amplitudes da solu¸c˜ao exata (3.33) e aproximada (3.38) podem

ter magnitudes bem diferentes quando o n´umero de onda k se aproxima de uma

ressonˆancia do problema original ou de uma ressonˆancia do problema aproximado,

o que causa uma s´eria deteriora¸c˜ao da qualidade da solu¸c˜ao aproximada.

Em problemas bidimensionais temos o mesmo tipo de diﬁculdade, com o

agravante de que, como podemos ver na equa¸c˜ao (2.21), a distˆancia entre dois

autovalores consecutivos ´e menor do que no problema unidimensional, aumentando

as chances de ocorrer ressonˆancia num´erica.

No problema de Helmholtz com condi¸c˜oes de contorno de Robin n˜ao ocorre

ressonˆancia, j´a que o problema de autovalor associado n˜ao tem autovalores re-

ais. Podemos fazer a mesma an´alise que ﬁzemos acima para veriﬁcar que no pro-

blema discreto com condi¸c˜oes de contorno de Robin tamb´em n˜ao ocorre ressonˆancia

num´erica. Para ilustrar isso, na Figura 3.3 apresentamos dois estudos de con-

vergˆencia para o problema bidimensional homogˆeneo. A solu¸c˜ao ´e uma onda plana

(equa¸c˜ao 2.32) na dire¸c˜ao θ = 0, o dom´ınio ´e o quadrado unit´ario (0, 1) × (0, 1).

Podemos notar que para condi¸c˜oes de Dirichlet a solu¸c˜ao aproximada, al´em de ser

afetada pela polui¸c˜ao, ´e tamb´em severamente afetada por ressonˆancia num´erica,

que afeta a convergˆencia monotˆonica da aproxima¸c˜ao.

3.5 Aproxima¸c˜ao por diferen¸cas ﬁnitas

Apresentaremos resumidamente m´etodos de diferen¸cas ﬁnitas para o pro-

blema de Helmholtz homogˆeneo em malhas uniformes, com o objetivo de mostrar

que, assim como as aproxima¸c˜oes por m´etodos elementos ﬁnitos cl´assicos deriva-

dos a partir da formula¸c˜ao de Galerkin, estas aproxima¸c˜oes por diferen¸cas ﬁnitas

sofrem igualmente os efeitos de polui¸c˜ao num´erica.

Consideramos inicialmente os m´etodos de segunda ordem obtidos classica-

mente usando s´erie de Taylor. Para introduzir uma aproxima¸c˜ao por diferen¸cas

ﬁnitas para o problema unidimensional, deﬁnimos as aproxima¸c˜oes para a derivada

primeira

j+1

− u

(3.43)

2 2.2 2.4 2.6 2.8 3

−log(h)

-2

-1.5

-1

-0.5

0.5

1.5

log (|u − u

|/|u|)

Interpolante

Galerkin

Condi¸c˜oes de Robin

2 2.2 2.4 2.6 2.8 3

−log(h)

-2

-1.5

-1

-0.5

0.5

1.5

log (|u − u

|/|u|)

Interpolante

Galerkin

Condi¸c˜oes de Dirichlet

Figura 3.3: Convergˆencia em H

, onda plana na dire¸c˜ao θ = 0, k = 100. Malhas

estruturadas de 100 × 100 at´e 1000 × 1000 elementos quadrados.

− u

j−1

(3.44)

e para a derivada segunda

j+1

− 2u

+ u

j−1

(3.45)

para as quais, usando a s´erie de Taylor obtemos

u(x

) =

u(x

+ h) − u(x

)

∂u(x

)

∂x

+ O(h) (3.46)

u(x

) =

u(x

) − u(x

− h)

∂u(x

)

∂x

+ O(h) (3.47)

u(x

) =

u(x

+ h) − 2u(x

) + u(x

− h)

∂

u(x

)

∂x

+ O(h

) (3.48)

Uma aproxima¸c˜ao por diferen¸cas ﬁnitas de segunda ordem, em um ponto

inerior j, para o problema de Helmholtz homogˆeneo unidimensional ´e ent˜ao deﬁnida

como

− D

− k

= 0, (3.49)

que pode ser posta na forma expl´ıcita

j−1

+ 2Su

+ Ru

j+1

= 0, (3.50)

com

R = −1 (3.51)

S = 1 −

(kh)

. (3.52)

Este ´e um esquema consistente com erro de truncamento da ordem de h

, isto ´e, a

solu¸c˜ao exata u(x) satisfaz a equa¸c˜ao

− D

− k

= τ

= O(h

). (3.53)

A an´alise de dispers˜ao desta aproxima¸c˜ao conduz `a seguinte estimativa para

o n´umero de ondas discreto

= k +

+ O(k

). (3.54)

Similarmente, para o problema bidimensional temos o esquema

− D

j,l

− D

j,l

− k

j,l

= 0, (3.55)

que pode ser apresentado na forma exp´ıcita

j−1,l+1

+ A

j,l+1

+ A

j+1,l+1

+ A

j−1,l

+ A

uj, l + A

j+1,l

+ A

j−1,l−1

+ A

j,l−1

+ A

j+1,l−1

= 0.

(3.56)

com

= 4 − h

(3.57)

= −1 (3.58)

= 0. (3.59)

Este esquema usual de segunda ordem envolvendo cinco pontos apresenta um erro

de truncamento de segunda ordem:

− D

u(x

, y

) − D

u(x

, y

) − k

u(x

, y

) = τ

j,l

= O(h

), (3.60)

e erro de fase relativo:

k − k

= −

3 + cos 4θ

(kh)

+ O



(kh)



, (3.61)

obtido como descrito na Se¸c˜ao 3.3.2.

Como podemos observar pelas estimativas obtidas para os n´umeros de onda

discretos, estas aproxima¸c˜oes por diferen¸cas ﬁnitas apresentam erros de fase de

mesma ordem que os correspondente aos m´etodos de elementos ﬁnitos de Galerkin

anteriormente apresentados.

Portanto, para estas aproxima¸c˜oes devemos esperar os mesmos efeitos de

polui¸c˜ao observados nas aproxima¸c˜oes de Galerkin. Aproxima¸c˜oes de diferen¸cas

ﬁnitas de ordens mais elevadas tˆem sido propostas para reduzir a polui¸c˜ao. Singer e

Turkel (1998) apresentam duas aproxima¸c˜oes de quarta ordem para o problema de

Helmholtz bidimensional sobre um estˆencil de nove pontos em malhas uniformes,

uma baseada na aproxima¸c˜ao de Pad´e e a outra constu´ıda com a utiliza¸c˜ao da

pr´opria equa¸c˜ao de Helmhotz. A t´ıtulo de compara¸c˜ao, para este segundo esquema

de quarta ordem proposto por Singer e Turkel, obtivemos a seguite estimativa para

o erro de fase relativo

k − k

5 − cos 4θ

2880

(kh)

+ O



(kh)



. (3.62)

Esquemas de sexta ordem tamb´em tˆem sido propostos (Sutmann (2007)).

Para estˆenceis de nove pontos, o QSFEM proposto por Babuˇska et al. (1995), e

que tamb´em ser´a apresentado neste cap´ıtulo, ´e o que apresenta polui¸c˜ao m´ınima.

3.6 Galerkin m´ınimos quadrados (GLS)

No m´etodo GLS (Harari e Hughes (1991), Thompson e Pinsky (1995)) s˜ao

acrescentados `a formula¸c˜ao de Galerkin termos constru´ıdos a partir do res´ıduo

da equa¸c˜ao (2.3) no interior dos elementos. Formalmente o m´etodo consiste em

encontrar u

∈ U

tal que

a(u

, v

) + τ a

, v

) = f(v

) + τ f

), ∀v

∈ V

, (3.63)

onde

, v

) =



Ω

(−∆u

− k

)(−∆v

− k

)dΩ (3.64)

) =



Ω

(−∆v

− k

)f dΩ (3.65)

Ω = Ω 



∂Ω

. (3.66)

O parˆametro τ pode ser ajustado para que o erro no n´umero de onda discreto

seja eliminado em malhas uniformes de elementos quadrados quando a solu¸c˜ao ´e

uma onda plana em uma dire¸c˜ao θ

escolhida. Normalmente a dire¸c˜ao escolhida ´e

, o que resulta no parˆametro de estabiliza¸c˜ao

τ =



1 − 6

4 − cos s

− cos s

− 2 cos s

cos s

(2 + cos s

)(2 + cos s



, (3.67)

com

= kh cos

e s

= kh sen

. (3.68)

Mas, para outras dire¸c˜oes θ = θ

, o erro no n´umero de onda discreto ´e da

mesma ordem que no m´etodo de Galerkin:

k − k

cos(4θ)

(kh)

+ O



(kh)



, (3.69)

o que pode afetar severamente a qualidade da aproxima¸c˜ao.

Outra diﬁculdade ´e que τ depende do produto kh, mas em elementos n˜ao

quadrados n˜ao existe uma escolha ´unica para o parˆametro h. A Figura 3.4 ilustra

quatro escolhas poss´ıveis e seus efeitos em uma malha com elementos retangulares.

3π

-1.5

-1

-0.5

log(u − u

/ u)

Interpolante

h = min{h

, h

}

h = max{h

, h

}

h =

√

Area

Figura 3.4: Dependˆencia do erro do m´etodo GLS em rela¸c˜ao a dire¸c˜ao θ da onda

plana, para diferentes deﬁni¸c˜oes do parˆametro m´edio h. Malha uniforme com

50 × 100 elementos retangulares, k = 50, condi¸c˜oes de Robin.

3.7 M´etodos com polui¸c˜ao m´ınima

Inspirados na an´alise de dispers˜ao apresentada anteriormente, Babuˇska et al.

(1995) desenvolveram uma aproxima¸c˜ao para o problema de Helmholtz homogˆeneo,

por eles denominado QSFEM (Quasi Stabilized Finite Element Method), que visa

minimizar a diferen¸ca k −k

. A seguir apresentamos resumidamente esse m´etodo.

3.7.1 Elementos ﬁnitos generalizados: QSFEM

O QSFEM ´e descrito por Babuˇska et al. (1995) como um m´etodo de ele-

mentos ﬁnitos generalizado para o problema de Helmholtz com polui¸c˜ao m´ınima.

No entanto, diferentemente de como s˜ao formulados normalmente os m´etodos de

elementos ﬁnitos, o QSFEM ´e deﬁnido algebricamente e n˜ao constru´ıdo a partir

de uma formula¸c˜ao variacional. Para deﬁnir o QSFEM, sup˜oe-se que o dom´ınio

seja discretizado por uma malha uniforme formada por elementos quadrados e que

toda equa¸c˜ao associada a um ponto nodal no interior do dom´ınio tenha a forma

− h, y

+ h) + G

, y

+ h) + G

+ h, y

+ h)

+ G

− h, y

) + G

, y

) + G

+ h, y

)

+ G

− h, y

− h) + G

, y

− h) + G

+ h, y

− h) = 0,

(3.70)

onde G

, G

e G

dependem apenas do produto kh e u

∈ S

(Ω) ´e a solu¸c˜ao

aproximada, dada como uma combina¸c˜ao linear das fun¸c˜oes base nodais φ

(x) =



)φ

(x). (3.71)

Outra forma de interpretar a equa¸c˜ao de dispers˜ao (3.24) ´e encontrada em

(Babuˇska e Sauter (1997)). Aplicando a transformada de Fourier na equa¸c˜ao de

Helmholtz homogˆenea ∆u + k

u = 0 em R

encontramos que a transformada

ˆu(s

, s

) satisfaz

+ s

− k

)ˆu(s

, s

) = 0, (3.72)

de onde conclu´ımos que o suporte de ˆu est´a contido no c´ırculo de raio k

{(s

, s

) ∈ R

; s

+ s

= k

}. (3.73)

Mas se considerarmos uma malha uniforme em todo o R

e uma fun¸c˜ao u

da forma

(3.71) satisfazendo a equa¸c˜ao (3.70) em todos os pontos nodais, ent˜ao o suporte

da transformada de Fourier de u

n˜ao est´a contido num c´ırculo e sim no conjunto

dos zeros da equa¸c˜ao de dispers˜ao, ou seja, na curva

{(s

, s

) ∈ R

; G

+2G



cos(hs

)+cos(hs

)



+4G

cos(hs

) cos(hs

) = 0}. (3.74)

O QSFEM ´e constru´ıdo (Babuˇska et al. (1995)) escolhendo os coeﬁcientes

, G

e G

que minimizam a distˆancia m´axima entre os c´ırculo (3.73) e a curva

(3.74). O estˆencil assim obtido faz com que essas duas curvas se interceptem em

16 pontos da forma (s

, s

) = (cos θ, sen θ) para

θ =

+ n

; n = 0, 1, . . . , 15. (3.75)

Outra forma de obter exatamente os mesmo coeﬁcientes ´e predeﬁnir G

= 4

e resolver o sistema de equa¸c˜oes em G

e G

+ 2G



cos(hR

) + cos(hR

)



+ 4G

cos(hR

) cos(hR

) = 0, (3.76)

+ 2G



cos(hS

) + cos(hS

)



+ 4G

cos(hS

) cos(hS

) = 0, (3.77)

onde

= k cos

, (3.78)

= k sen

, (3.79)

= k cos

3π

, (3.80)

= k sen

3π

. (3.81)

O resultado ´e

= 2

− c

+ s

) − c

+ s

)

, (3.82)

+ s

− c

− s

+ s

) − c

+ s

)

, (3.83)

com

= cos



kh cos



, (3.84)

= cos



kh sen



, (3.85)

= cos



kh cos

3π



, (3.86)

= cos



kh sen

3π



. (3.87)

Com esses coeﬁcientes, procedendo como descrito na Se¸c˜ao 3.3.2 obtemos:

k − k

= −

cos 8θ

774144

(kh)

+ O



(kh)



, (3.88)

ou seja, o erro relativo no n´umero de onda para o QSFEM ´e de sexta ordem em

kh, isso ´e o melhor que se pode obter para um estˆencil de 9 pontos da forma (3.70).

Essa ´ultima forma de determinar o estˆencil ´otimo foi posteriormente usada

para obter parˆametros de estabiliza¸c˜ao para outros m´etodos, por exemplo, os pro-

postos por Oberai e Pinsky (2000), Loula et al. (2007), Harari e Gosteev (2007),

do Carmo et al. (2008).

Apesar de n˜ao se observar visualmente polui¸c˜ao nos experimentos feitos com

o QSFEM, sabe-se (Babuˇska e Sauter (1997), Teorema 5.6) que ´e imposs´ıvel elimi-

nar totalmente o efeito de polui¸c˜ao em problemas bidimensionais. N˜ao se conhece

uma estimativa da ordem do erro em L

ou H

devido `a polui¸c˜ao para o QSFEM

ou para outros m´etodos com estˆenceis equivalentes em dom´ınio limitados. Na ver-

dade uma an´alise num´erica completa para tais m´etodos n˜ao foi apresentada, isso

inclui os m´etodos que apresentaremos adiante.

Como o QSFEM n˜ao ´e constru´ıdo a partir de uma formula¸c˜ao variacional,

n˜ao ´e poss´ıvel aplic´a-lo diretamente a malhas mais gerais ou a problemas n˜ao

homogˆeneos ou com condi¸c˜oes de contorno mais gerais. M´etodos variacionalmente

consistentes que geram estˆenceis de 9 pontos equivalentes ao do QSFEM quando

restritos a malhas uniformes de elementos quadrados foram propostos de v´arias

formas. Como exemplo apresentamos a seguir dois desses m´etodos (Fernandes e

Loula (2006)), desenvolvidos durante os primeiros est´agios deste trabalho.

3.7.2 Galerkin, volumes ﬁnitos e m´ınimos quadrados

No m´etodo GLSFV (Galerkin + Least Squares + Finite Volume), um se-

gundo parˆametro livre ´e introduzido adicionando-se outro termo de res´ıduo ao

m´etodo GLS.

Considere uma malha sobreposta `a malha de elementos ﬁnitos, conforme a

Figura 3.5. Essa malha sobreposta ´e constru´ıda de modo que cada i-´esimo n´o da

malha esteja em um ´unico subconjunto Q

⊂ Ω. Por exemplo, na Figura 3.5, ao

n´o central est´a associado o quadrado cinza. Estando deﬁnidos esses conjuntos Q

podemos deﬁnir um operador linear P : S

→ L

(Ω) do seguinte modo. Para as

fun¸c˜oes bases nodais φ

, P ´e dado por

P φ

(x) =











1, se x ∈ Q

0, se x /∈ Q

(3.89)

Figura 3.5: Malha de volume ﬁnitos (linha tracejada) sobreposta `a malha de ele-

mentos ﬁnitos (linha cont´ınua).

Ent˜ao, para qualquer fun¸c˜ao v



i=1

de S

, temos que

P v



i=1

P φ

. (3.90)

O m´etodo GLSFV consiste em encontrar u

∈ U

tal que

, v

) + αB

, v

) + βB

F V

, v

)

= F

) + αF

) + βF

F V

) (3.91)

para todo v

∈ V

, onde

F V

, v

) = −



Ω

P v

dΩ −



∂Q

P v

∇u

ds, (3.92)

F V

) =



Ω

fP v

dΩ, (3.93)

∂Q =



∂Q

, (3.94)

P v

∇u

 = (P v

∇u

)

· n

+ (P v

∇u

)

−

· n

−

(3.95)

3.7.2.1 Consistˆencia

O m´etodo GLSFV ´e consistente no sentido de que a solu¸c˜ao exata u do

problema variacional (2.8) satisfaz a equa¸c˜ao (3.91). Como o m´etodo GLS ´e con-

sistente, e considerando que todo elemento v

∈ V

´e combina¸c˜ao linear de fun¸c˜oes

base nodais φ

, para veriﬁcar a consistˆencia do m´etodo GLSFV basta mostrar que

F V

(u, φ

) − F

F V

(φ

) = 0 para toda fun¸c˜ao φ

. De fato,

F V

(u, φ

) − F

F V

(φ

) = −



Ω

uP φ

dΩ −



∂Q

P φ

∇u · nds −



Ω

fP φ

dΩ

= −



u dΩ −



∂Q

∇u · nds −



f dΩ (3.96)

= −



u dΩ −



∆u dΩ −



f dΩ (3.97)



−k

u − ∆u − f dΩ = 0 (3.98)

o que prova a consistˆencia da formula¸c˜ao GLSFV.

3.7.2.2 Determina¸c˜ao dos parˆametros de estabiliza¸c˜ao

Os parˆametros de estabiliza¸c˜ao s˜ao determinados a partir da an´alise usada

para construir o QSFEM. No problema homogˆeneo, com uma malha uniforme de

elementos bilineares de lados h, as equa¸c˜oes associadas aos pontos nodais interiores

tˆem a mesma foram geral (3.70) do QSFEM, com

−

+ α

+ β



−



, (3.99)

= −

−

+ α

+ β



−



, (3.100)

= −

−

+ α

+ β



−



. (3.101)

Neste caso os coeﬁcientes G

, G

e G

dependem linearmente dos parˆametros

livres α e β.

Como no QSFEM, determinamos os parˆametros ´otimos for¸cando a equa¸c˜ao

de dispers˜ao a ser satisfeita com k

= k para dire¸c˜oes θ

e θ

devidamente escolhi-

das. Isso nos leva a um sistema de duas equa¸c˜oes lineares em duas inc´ognitas α e

β. Resolvendo esse sistema para θ

e θ

3π

encontramos:

α =

− b

), (3.102)

β =

− b

(−a

+ a

). (3.103)

Onde:

= A

+ A

+ S

) + A

, (3.104)

= A

F V

+ A

F V

+ S

) + A

F V

, (3.105)

= A

+ A

+ S

) + A

, (3.106)

(3.107)

= cos(kh cos θ

), (3.108)

= cos(kh sen θ

), (3.109)

−

, A

= −

−

, A

= −

−

, (3.110)

, A

, (3.111)

F V

= −

−

, A

F V

−

, A

F V

−

. (3.112)

(3.113)

Expandindo α(kh) e β(kh) em s´erie de Taylor em torno de kh = 0 obtemos

α =



+ α

(kh)

+ α

(kh)

+ α

(kh)

+ . . .



(3.114)

β = β

+ β

(kh)

+ β

(kh)

+ β

(kh)

+ . . . (3.115)

(3.116)

com

= 4, α

, α

137

3600

, α

10657

2304000

, (3.117)

= 2, β

= −

120

, β

3200

, β

9713

6912000

. (3.118)

Para o problema homogˆeneo de Dirichlet, discretizando o dom´ınio com uma

malha uniforme de elementos quadrados, o GLSFV produz os mesmos resultados

que o QSFEM.

3.7.3 Galerkin e m´ınimos quadrados ponderados

Outra formula¸c˜ao consistente e com as mesmas propriedades de estabiliza¸c˜ao

do m´etodo GLSFV ´e o m´etodo GLS com um segundo termo de m´ınimos quadrados

ponderado por uma fun¸c˜ao peso que varia no interior dos elementos. Esse m´etodo

consiste em encontrar u

∈ U

tal que

, v

) + αB

, v

) + βB

, v

) (3.119)

= F

) + αF

) + βF

), ∀v

∈ V

, (3.120)

onde

, v

) =



Ω

P (x, y)(−∆u

− k

)(−∆v

− k

)dΩ, (3.121)

) =



Ω

P (x, y)(−∆v

− k

)f dΩ, (3.122)

Ω = Ω −



∂Ω

. (3.123)

Em cada elemento a fun¸c˜ao de peso P ´e dada em coordenadas locais:

P (ξ, η) =

(1 − η

) (1 − ξ

)

, (3.124)

ou seja, P (ξ, η) ´e uma bolha biquadr´atica deﬁnida em cada elemento.

Como no m´etodo anterior, utilizamos a equa¸c˜ao de dispers˜ao para encontrar

os parˆametros ´otimos:

α =

(kh)



+ α

(kh)

+ α

(kh)

+ α

(kh)

+ . . .



(3.125)

β =

(kh)



+ β

(kh)

+ β

(kh)

+ β

(kh)

+ . . .



(3.126)

com

= −25, α

, α

5184

, α

2023

933120

, (3.127)

225

, β

= −

, β

= −

155

2304

, β

= −

749

82944

. (3.128)

3.8 Observa¸c˜oes

Conforme j´a mencionado, o objetivo desta tese ´e formular m´etodos de ele-

mentos ﬁnitos e de diferen¸cas ﬁnitas para o problema de Helmholtz que sejam

quase est´aveis, no sentido deﬁnido por Babuska, isto ´e, capazes de minimizar ou

pelo menos reduzir signiﬁcativamente o erro de polui¸c˜ao t´ıpico destas aproxima¸c˜oes

e razoavelmente robustos quando aplicados a malhas mais gerais do que ma-

lhas uniformes, incluindo malhas n˜ao estruturadas. A parte destacada ´e o

maior desaﬁo, pois, como vimos, para malhas uniformes de elementos quadrados,

a quest˜ao da polui¸c˜ao est´a razoavelmente bem resolvida.

Os m´etodos de elementos ﬁnitos que apresentamos neste cap´ıtulo possuem

parˆametros livres que s˜ao ajustados para malhas uniformes com elementos qua-

drados de lado h. Essa t´ecnica ´e bastante comum nos m´etodos estabilizados para

o problema de Helmholtz. No entanto, para malhas mais gerais a an´alise de dis-

pers˜ao que ´e usada para determinar os parˆametros de estabiliza¸c˜ao n˜ao se aplica

diretamente. Nesses casos pode ser calculado um comprimento m´edio h para cada

elemento e ent˜ao os parˆametros de estabiliza¸c˜ao s˜ao calculados como se os ele-

mentos fossem quadrados. Entretanto, experimentos num´ericos indicam que os

parˆametros assim obtidos est˜ao longe dos valores ´otimos quando se trata de ma-

lhas distorcidas. Dada a grande sensibilidade da solu¸c˜ao aproximada a varia¸c˜oes no

n´umero de onda discreto, estas estrat´egias baseadas em parˆametros ´otimos obtidos

para malhas uniformes n˜ao funcionam muito bem. Assim, tendo em vista nosso ob-

jetivo, temos que contornar essa dependˆencia em rela¸c˜ao a um ´unico comprimento

No pr´oximo cap´ıtulo apresentamos duas formula¸c˜oes de elementos ﬁnitos de

Petrov-Galerkin que s˜ao, pelo menos em parte, livres dessa limita¸c˜ao. Na primeira

n˜ao h´a parˆametros de estabiliza¸c˜ao dependentes de h. Na segunda, os parˆametros

dependem diretamente de todos os comprimentos dos lados dos elementos.

Para os m´etodos de diferen¸cas ﬁnitas que descrevemos brevemente neste

cap´ıtulo a situa¸c˜ao ´e ainda mais complicada, visto que esses m´etodos s˜ao formu-

lados em termos de estˆenceis que s˜ao pouco ﬂex´ıveis em rela¸c˜ao a distribui¸c˜ao dos

pontos onde a solu¸c˜ao ´e aproximada, como ´e normalmente o caso de m´etodos de

diferen¸cas ﬁnitas usuais. No Cap´ıtulo 5 apresentamos uma formula¸c˜ao diferente

nesse sentido, onde n˜ao ´e necess´aria uma distribui¸c˜ao uniforme dos pontos onde a

solu¸c˜ao ´e aproximada. Um crit´erio diferente do usado no QSFEM ser´a usado para

deﬁnir os estˆenceis ´otimos. Este novo crit´erio ser´a usado tamb´em no Cap´ıtulo 6

para formular um m´etodo de elementos ﬁnitos de Petrov-Galerkin que apresenta

grande robustez em rela¸c˜ao a distor¸c˜oes da malha.

Cap´ıtulo 4

M´etodos de Petrov-Galerkin

Neste cap´ıtulo s˜ao apresentados dois m´etodos de Petrov-Galerkin para a

equa¸c˜ao de Helmholtz. Apenas dom´ınios bidimensionais s˜ao considerados. Nos

dois m´etodos o espa¸co U

⊂ S

= span{φ

, i = 1, . . . , N} ´e o espa¸co de elementos

ﬁnitos com elementos bilineares do m´etodo de Galerkin da Se¸c˜ao 3.1, enquanto

que o espa¸co V

´e gerado por fun¸c˜oes base nodais diferentes das do espa¸co U

= span{ψ

, i = 1, . . . , N

}, (4.1)

com ψ

= φ

, como nos m´etodos de Petrov-Galerkin em geral. Cada fun¸c˜ao base

nodal ψ

ser´a constru´ıda para ter o mesmo suporte que a fun¸c˜ao φ

correspondente,

supp(ψ

) = supp(φ

), i = 1 . . . N

, (4.2)

resultando em matrizes globais esparsas com o mesmo perﬁl (os mesmos termos

nulos) da matriz do m´etodo de Galerkin.

4.1 M´etodo com polui¸c˜ao reduzida (RPPG)

O RPPG (Reduced Pollution Petrov-Galerkin (FEM)) ´e o mais simples dos

m´etodos propostos neste trabalho. Trata-se de uma vers˜ao truncada do m´etodo

QSPG da pr´oxima se¸c˜ao, mas como foi obtido primeiro, ser´a apresentado separa-

damente.

Nos m´etodos de elementos ﬁnitos, as fun¸c˜oes base dos espa¸cos U

e V

s˜ao

constru´ıdas por partes elemento por elemento. Em coordenadas locais ξ do ele-

mento de referˆencia, denotamos por φ

(ξ) a fun¸c˜ao base associada ao n´o local i .

Por exemplo, em 1D, temos duas fun¸c˜oes locais para elementos lineares:

(ξ) =

1 − ξ

, (4.3)

(ξ) =

1 + ξ

, (4.4)

onde −1 ≤ ξ ≤ 1 ´e a coordenada local do dom´ınio de referˆencia.

Em 2D, para elementos bilineares, temos

(ξ, η) =



1 − ξ



1 − η



, (4.5)

(ξ, η) =



1 + ξ



1 − η



, (4.6)

(ξ, η) =



1 + ξ



1 + η



, (4.7)

(ξ, η) =



1 − ξ



1 + η



, (4.8)

com −1 ≤ ξ ≤ 1 e −1 ≤ η ≤ 1.

O RPPG ´e deﬁnido apenas para elementos quadrilaterais. O espa¸co V

´e

constru´ıdo de forma an´aloga ao U

. As fun¸c˜oes base nodais ψ

do espa¸co V

tamb´em s˜ao deﬁnidas por partes em coordenadas locais dos elementos, de forma

semelhante `as fun¸c˜oes φ

. Deﬁnimos ψ

como o produto de fun¸c˜oes de forma

unidimensionais p

(t) e p

(t):

(ξ, η) = p

(ξ) p

(η), (4.9)

(ξ, η) = p

(ξ) p

(η), (4.10)

(ξ, η) = p

(ξ) p

(η), (4.11)

(ξ, η) = p

(ξ) p

(η). (4.12)

Onde p

e p

s˜ao deﬁnidas como

(t) =

2 − 7t + 5t

, (4.13)

(t) =

2 + 7t − 5t

. (4.14)

Os polinˆomios p

e p

, cujos gr´aﬁcos est˜ao apresentados na Figura 4.1, foram

obtidos inicialmente de forma emp´ırica. A id´eia de procur´a-los surgiu quando se

tentava modiﬁcar o esquema de integra¸c˜ao num´erica com o intuito de obter al-

guma estabiliza¸c˜ao para a formula¸c˜ao de Galerkin levando em conta as distor¸c˜oes

dos elementos naturalmente. Essa possibilidade foi testada empiricamente, modiﬁ-

cando os pesos dos pontos de integra¸c˜ao e observando os resultados, que indicaram

que certas escolhas dos pesos levavam a ganhos signiﬁcativos de estabilidade. Mas

ent˜ao notou-se que o efeito causado por uma modiﬁca¸c˜ao nos pesos dos pontos de

integra¸c˜ao poderia ser obtido tamb´em modiﬁcando o integrando e fazendo a inte-

gral exata, o que ´e perfeitamente poss´ıvel, no contexto dos m´etodos dos elementos

ﬁnitos, modiﬁcando os espa¸cos U

e/ou V

A id´eia inicial de modiﬁcar o esquema de integra¸c˜ao n˜ao foi levada adiante,

mas ´e interessante citar aqui o trabalho de Guddati e Yue (2004). Nesse trabalho

as coordenadas dos pontos de integra¸c˜ao (e n˜ao os pesos) s˜ao modiﬁcadas com o

intuito de reduzir o erro no n´umero de onda discreto. Uma conex˜ao interessante

com esse trabalho ´e que, para elementos retangulares, a matriz local resultante do

m´etodo proposto por Guddati e Yue (2004) ´e igual `a do m´etodo RPPG que aqui

propomos.

No Cap´ıtulo 3 foram apresentadas as seguintes estimativas para o erro no

n´umero de onda discreto em malhas uniformes de elementos de lado h para o

m´etodo de Galerkin

k − k

3 + cos 4θ

(kh)

+ O



(kh)



, (4.15)

−1 0 1

(t)

1−t

1+t

Figura 4.1: Gr´aﬁcos dos polinˆomios p

e p

do m´etodo RPPG e os gr´aﬁcos das

fun¸c˜oes correspondentes no m´etodo de Galerkin.

para o GLS

k − k

cos 4θ

(kh)

+ O



(kh)



(4.16)

e para o QSFEM

k − k

= −

cos 8θ

774144

(kh)

+ O



(kh)



. (4.17)

O mesmo c´alculo usado para obter as estimativas acima foi usado para obter o

seguinte resultado para o RPPG:

k − k

= −

5 + 3 cos 4θ

3840

(kh)

+ O



(kh)



. (4.18)

Comparando os quatro resultados vemos que o erro relativo no n´umero de onda do

RPPG ´e de quarta ordem em rela¸c˜ao a kh, o que ´e bem melhor comparado com

o m´etodo de Galerkin e o GLS, que s˜ao de segunda ordem. Mas esse erro ainda ´e

bem maior do que o m´ınimo obtido com o QSFEM, que ´e de sexta ordem.

Nas Figuras 4.2 e 4.3 o RPPG ´e comparado com o m´etodo de Galerkin em um

quadrado unit´ario com condi¸c˜oes de contorno de Robin e malhas uniformes com

elementos quadrados. Nota-se uma signiﬁcativa redu¸c˜ao no efeito de polui¸c˜ao. Na

norma do L

a aproxima¸c˜ao do RPPG chega a ser melhor do que a interpolante.

Outros resultados num´ericos ser˜ao apresentados adiante.

-1

-0.5

log(u − u

/ u)

Interpolante

Galerkin

RPPG

Figura 4.2: Dependˆencia do erro na norma L

em rela¸c˜ao a dire¸c˜ao θ da onda

plana, k = 100, condi¸c˜oes contorno de Robin, malha uniforme com 100 × 100

elementos quadrados.

2 2.2 2.4 2.6

−log(h)

-2.5

-2

-1.5

-1

-0.5

0.5

log (u − u

/u)

Interpolante

Galerkin

RPPG

2 2.2 2.4 2.6

−log(h)

-1.5

-1

-0.5

0.5

log (|u − u

|/|u|)

Interpolante

Galerkin

RPPG

Figura 4.3: Convergˆencia para uma onda plana com θ =

e k = 100, condi¸c˜oes

de contorno de Robin, malhas uniformes de elementos quadrados.

Uma caracter´ıstica interessante do RPPG ´e que n˜ao existem parˆametros

dependentes do n´umero de onda ou dos tamanhos dos elementos. Isso ´e interessante

do ponto de vista pr´atico pois a implementa¸c˜ao do RPPG pode ser feita com poucas

modiﬁca¸c˜oes a partir de uma implementa¸c˜ao do m´etodo de Galerkin. O m´etodo

QSPG apresentado na pr´oxima se¸c˜ao inclui parˆametros dependendes de k e dos

comprimentos dos lados dos elementos. Veremos que com isso ´e poss´ıvel obter

polui¸c˜ao m´ınima para malhas uniformes com elementos quadrados.

4.2 M´etodo quase estabilizado (QSPG)

Considere o elemento de referˆencia ilustrado na Figura 4.4. Os comprimentos

, h

e h

indicados s˜ao os obtidos a partir das coordenadas globais dos n´os.

= (−1, −1) ξ

= (1, −1)

= (1, 1)ξ

= (−1, 1)

Figura 4.4: Elemento de referˆencia.

Novamente constru´ımos o espa¸co V

usando fun¸c˜oes base nodais ψ

diferentes

das φ

do espa¸co U

. Deﬁnimos fun¸c˜oes nodais em coordenadas locais do elemento

de referˆencia da seguinte forma:

(ξ, η) = P



ξ, α(kh

), β(kh

)





η, α(kh

), β(kh

)



, (4.19)

(ξ, η) = P



ξ, α(kh

), β(kh

)





η, α(kh

), β(kh

)



, (4.20)

(ξ, η) = P



ξ, α(kh

), β(kh

)





η, α(kh

), β(kh

)



, (4.21)

(ξ, η) = P



ξ, α(kh

), β(kh

)





η, α(kh

), β(kh

)



, (4.22)

onde P

(t, α, β) e P

(t, α, β) s˜ao os polinˆomios de grau 3 em t

(t, α, β) = α + βt +



− α



−



+ β



, (4.23)

(t, α, β) = α − βt +



− α





+ β



, (4.24)

com parˆametros α e β a ser deﬁnidos.

Os polinˆomios P

(t, α, β) e P

(t, α, β) satisfazem, independentemente dos

parˆametros α e β, as seguintes condi¸c˜oes:

(−1, α, β) = 1, (4.25)

(1, α, β) = 0, (4.26)

(−1, α, β) = 0, (4.27)

(1, α, β) = 1, (4.28)

(t, α, β) = P

(−t, α, β). (4.29)

A dependˆencia em rela¸c˜ao a α e β poderia ser de outras formas, a escolha

particular das equa¸c˜oes 4.23 e 4.24 implica que

(0, α, β) = P

(0, α, β) = α e (4.30)



(0, α, β) = −P



(0, α, β) = β. (4.31)

Se escolhermos parˆametros constantes α =

e β = −

, recuperamos o

m´etodo de Galerkin. O RPPG ´e obtido escolhendo α =

e β = −

. No caso

do QSPG, os parˆametros α e β dependem de kh, a deﬁni¸c˜ao precisa ser´a dada na

Se¸c˜ao 4.2.2. Independentemente de como eles s˜ao escolhidos, as fun¸c˜oes do espa¸co

resultante s˜ao cont´ınuas. Isso ´e veriﬁcado com mais cuidado a seguir.

4.2.1 Conformidade

Por constru¸c˜ao, as fun¸c˜oes globais ψ

s˜ao cont´ınuas no interior dos elementos.

A continuidade nas fronteiras dos elementos ´e consequˆencia das condi¸c˜oes impostas

nas equa¸c˜oes (4.25)-(4.29).

As fun¸c˜oes locais ψ

s˜ao constru´ıdas de modo que as formas que elas assumem

em algum lado do elemento dependem apenas do produto do comprimento desse

lado pelo n´umero de onda, conforme est´a expl´ıcito nas equa¸c˜oes seguintes, que s˜ao

veriﬁcadas diretamente a partir das deﬁni¸c˜oes:

Lado 1:











(ξ, η) = P

(ξ, α(kh

), β(kh

))

(ξ, η) = P

(−ξ, α(kh

), β(kh

))

(ξ, η) = 0

(4.32)

Lado 2:











(ξ, η) = 0

(ξ, η) = P

(η, α(kh

), β(kh

))

(ξ, η) = P

(−η, α(kh

), β(kh

)),

(ξ, η) = 0

(4.33)

Lado 3:











(ξ, η) = 0

(ξ, η) = P

(−ξ, α(kh

), β(kh

))

(ξ, η) = P

(ξ, α(kh

), β(kh

))

(4.34)

Lado 4:











(ξ, η) = P

(η, α(kh

), β(kh

))

(ξ, η) = 0

(ξ, η) = P

(−η, α(kh

), β(kh

))

(4.35)

Al´em disso as fun¸c˜oes globais ψ

(x, y) n˜ao dependem da numera¸c˜ao local dos

n´os. De fato, partindo de um elemento com coordenadas locais (ξ, η), ent˜ao uma

renumera¸c˜ao dos n´os R = (1, 2, 3, 4) → (4, 1, 2, 3) resulta em um novo sistema de

coordenadas locais (ξ



, η



) = (η, −ξ); se um ponto global (x, y) tem coordenadas

(ξ, η) no elemento original, ent˜ao esse mesmo ponto tem coordenadas (ξ



, η



) no

elemento com a nova numera¸c˜ao. A mesma renumera¸c˜ao aplica-se para os compri-

mentos dos lados, ou seja, h

= h



, h

= h



, h

= h



e h

= h



. A fun¸c˜ao global ψ

associada ao n´o 1 na numera¸c˜ao local original ´e dada por

ψ(x, y) = ψ(x(ξ, η), y(ξ, η)) = ψ

(ξ, η) (4.36)

= P

(ξ, α(kh

), β(kh

)) · P

(η, α(kh

), β(kh

)). (4.37)

Com os n´os renumerados, temos que

ψ(x, y) = ψ(x(ξ



, η



), y(ξ



, η



)) = ψ

(ξ



, η



) (4.38)

= P

(ξ



, α(kh



), β(kh



)) · P

(η



, α(kh



), β(kh



)) (4.39)

= P

(η, α(kh

), β(kh

)) · P

(−ξ, α(kh

), β(kh

)) (4.40)

= P

(η, α(kh

), β(kh

)) · P

(ξ, α(kh

), β(kh

)) (4.41)

= P

(ξ, α(kh

), β(kh

)) · P

(η, α(kh

), β(kh

)), (4.42)

que ´e o mesmo resultado encontrado em (4.37).

O mesmo pode ser veriﬁcado para as fun¸c˜oes globais associadas aos outros

n´os. Como qualquer renumera¸c˜ao v´alida dos n´os pode ser obtida por aplica¸c˜oes

iterativas de R, ent˜ao conclu´ımos que as fun¸c˜oes globais n˜ao mudam quando uma

renumera¸c˜ao local dos n´os ´e feita.

Essa independˆencia em rela¸c˜ao a numera¸c˜ao dos n´os nos permite veriﬁcar a

continuidade analisando, sem perda de generalidade, apenas a situa¸c˜ao ilustrada

na Figura 4.5. Nessa ﬁgura x

, . . . , x

s˜ao n´os globais aos quais est˜ao associadas

as fun¸c˜oes base globais ψ

, . . . , ψ

respectivamente. O segmento x

pode ser

parametrizado com t ∈ [−1, 1] por

x(t) =

1 − t

1 + t

. (4.43)

Com a transforma¸c˜ao usual entre coordenadas locais e globais, no elemento e um

ponto x(t) tem coordenadas locais

ξ(t) =

1 − t

1 + t

= (1, t), (4.44)

enquanto que no elemento e



o mesmo ponto tem coordenadas locais



(t) =

1 − t



1 + t



= (−1, t). (4.45)

1 2

′

1 2

Figura 4.5: Elementos vizinhos.

A continuidade ao longo do segmento x

segue das equa¸c˜oes (4.33) e (4.35).

De fato, temos que as fun¸c˜oes globais ψ

, ψ

e ψ

se anulam ao longo do

segmento considerado. Isso garante a continuidade dessas fun¸c˜oes, j´a que, por

constru¸c˜ao, esse segmento ´e parte da fronteira de seus suportes. Quanto `as fun¸c˜oes

e ψ

, no lado do elemento e temos que

(x(t)) = ψ

(ξ(t)) = ψ

(1, t) = P

(t, α(x

− x

), β(x

− x

)), (4.46)

(x(t)) = ψ

(ξ(t)) = ψ

(1, t) = P

(−t, α(x

− x

), β(x

− x

)). (4.47)

o mesmo resultado ocorre no lado do elemento e



(x(t)) = ψ

(ξ



(t)) = ψ

(−1, t) = P

(t, α(x

− x

), β(x

− x

)), (4.48)

(x(t)) = ψ

(ξ



(t)) = ψ

(−1, t) = P

(−t, α(x

− x

), β(x

− x

)), (4.49)

o que nos permite concluir que ψ

e ψ

s˜ao cont´ınuas ao longo do segmento consi-

derado.

4.2.2 Parˆametros de estabiliza¸c˜ao

Para completar a formula¸c˜ao do QSPG falta apenas deﬁnir os parˆametros de

estabiliza¸c˜ao α(kh) e β(kh).

Em um elemento geral, temos na verdade oito parˆametros independentes,

α(kh

), α(kh

), β(kh

) e β(kh

), sendo dois as-

sociados a cada lado do elemento. Conforme a se¸c˜ao anterior, para assegurar a

continuidade atrav´es da fronteira dos elementos, ´e preciso calcular esses parˆemtros

de forma consistente para todos os elementos, se usarmos uma deﬁni¸c˜ao de α(kh) e

β(kh) para um elemento, ent˜ao ´e preciso usar a mesma deﬁni¸c˜ao em todos outros.

Em particular, a deﬁni¸c˜ao usada para elementos quadrados tem que ser usada

tamb´em nos demais elementos.

O crit´erio adotado para deﬁnir α(kh) e β(kh) ´e que, quando restrito a malhas

uniformes com elementos quadrados, o QSPG tenha polui¸c˜ao m´ınima no mesmo

sentido que o QSFEM.

E dessa forma de determinar esses parˆametos que vem a

denomina¸c˜ao “QSPG”: Quasi Stabilized Petrov-Galerkin (FEM).

No caso de um elemento quadrado temos que

= h

= h, (4.50)

logo os oito parˆametros acabam reduzindo-se a apenas dois:

α(kh

) = α(kh

) = α(kh) e (4.51)

β(kh

) = β(kh

) = β(kh). (4.52)

Assim como no QSFEM, bem como em todos os m´etodos de elementos ﬁnitos

bilineares descritos anteriormente, o estˆencil do QSPG com malhas de elementos

quadrados pode ser representado pela matriz













. (4.53)

Os coeﬁcientes A

, A

e A

associados ao QSPG s˜ao

225



30 − ah



, (4.54)

= −

225



30a − 15b + abh



, (4.55)

= −

900



60 + bh



(4.56)

com

a = 2 + 5α(kh) − β(kh), (4.57)

b = 1 + 10α(kh) + 2β(kh). (4.58)

Na Se¸c˜ao 3.7.1 os coeﬁcientes do estˆencil do QSFEM foram obtidos resol-

vendo a equa¸c˜ao de dispers˜ao (3.24) em rela¸c˜ao a A

e A

com k

= k para

as duas dire¸c˜oes θ =

e θ =

3π

. No QSPG temos que A

, A

e A

dependem dos

dois parˆametros livres α(kh) e β(kh), ent˜ao podemos proceder como no caso do

QSFEM mas resolvendo em rela¸c˜ao a α(kh) e β(kh). Assim, α(kh) e β(kh) s˜ao

determinados resolvendo o sistema de equa¸c˜oes n˜ao lineares

(α, β) + 2R

(α, β) + 4S

(α, β) = 0, (4.59)

(α, β) + 2R

(α, β) + 4S

(α, β) = 0, (4.60)

onde

= cos



kh cos



+ cos



kh sen



, (4.61)

= cos



kh cos



cos



kh sen



, (4.62)

= cos



kh cos

3π



+ cos



kh sen

3π



, (4.63)

= cos



kh cos

3π



cos



kh sen

3π



. (4.64)

As equa¸c˜oes (4.59)-(4.60) podem ser resolvidas substituindo (4.54)-(4.56) em

(4.59)-(4.60), resolvendo em rela¸c˜ao a a e b e ﬁnalmente usando (4.57)-(4.58) para

encontrar

α(kh) =

−5 + a + 2b

, (4.65)

β(kh) =

3 + a − 2b

. (4.66)

Encontramos assim 4 solu¸c˜oes distintas.

As duas primeiras solu¸c˜oes s˜ao

a = 0 e b = 0, (4.67)

a =

e b = −

. (4.68)

Estas s˜ao solu¸c˜oes esp´urias, pois implicam em A

= A

= 0, sendo ent˜ao

descartadas.

Para escrever as outras duas solu¸c˜oes, usamos as seguintes express˜oes auxi-

liares:

δ = R

+ (S

− S

)

+ R

− R

+ S

), (4.69)

 = S

(1 − R

) + S

− 1), (4.70)

γ = R

− R

− S

+ S

. (4.71)

Assumimos que δ > 0,  > 0 e γ > 0, que ´e verdadeiro para 0 < kh ≤ 1 (mas n˜ao

s´o nesse intervalo), o que nos permite fazer algumas simpliﬁca¸c˜oes.

A terceira solu¸c˜ao pode ent˜ao ser escrita como

a = 15

1 +



√

, (4.72)

b = 30

1 +

√

. (4.73)

Com ajuda do software Mathematica (Wolfram Research (2005)) foi poss´ıvel calcu-

lar as seguintes expans˜oes em s´eries para os parˆametros resultantes dessa solu¸c˜ao:

α(kh) = −1 −

640

(kh)

−

73728

(kh)

−

409600

(kh)

+ O



(kh)



, (4.74)

β(kh) = −

(kh)

−

(kh)

1152

(kh)

+ O



(kh)



. (4.75)

Podemos notar nessas express˜oes que P

(0) = α(kh) → −1 quando h → 0. Mas

o termo

(kh)

na expans˜ao de β(kh) implica que |P



(0)| cresce indeﬁnidamente

quando h → 0.

A quarta solu¸c˜ao ´e dada por

a = 15

1 −



√

, (4.76)

b = −30

1 −

√

, (4.77)

de onde obtemos

α(kh) =

640

(kh)

73728

(kh)

409600

(kh)

+ O



(kh)



, (4.78)

β(kh) = −

(kh)

−

1152

(kh)

552960

(kh)

+ O



(kh)



. (4.79)

Nesta ´ultima solu¸c˜ao tanto α(kh) quanto β(kh) convergem para valores ﬁnitos

quando h → 0, por isso os parˆametros escolhidos para o QSPG s˜ao esses.

Para o QSPG, o erro no n´umero de onda para malhas uniformes com ele-

mentos quadrados ´e idˆentico ao do QSFEM.

Se usarmos apenas o termo independente de kh em cada uma das s´eries

(4.78) e (4.79), ou seja, α =

e β = −

, ent˜ao reca´ımos no RPPG. Assim, quando

h → 0, o RPPG se aproxima do QSPG, isso pode ser observado nos experimentos

num´ericos que apresentamos adiante.

4.3 Implementa¸c˜ao computacional

O RPPG ´e bastante simples de ser implementado a partir de uma imple-

menta¸c˜ao do m´etodo de Galerkin. Mas dois cuidados devem ser tomados. O

primeiro ´e aumentar o n´umero de pontos de integra¸c˜ao num´erica, visto que s˜ao

usados polinˆomios de grau mais alto. Nos experimentos que apresentamos a seguir,

foram usados 3 × 3 pontos de integra¸c˜ao de Gauss-Legendre. Outro ponto a ser

observado ´e que a matriz global do RPPG, diferentemente do m´etodo de Galer-

kin, n˜ao ´e sim´etrica em geral, logo o m´etodo usado para resolver o sistema linear

associado deve ser apropriado para esse tipo de matriz mais geral.

As observa¸c˜oes acima tamb´em valem para o QSPG, que ainda requer um

trabalho extra. A forma exata das fun¸c˜oes base em coordenadas locais varia de

elemento para elemento. Assim, ´e necess´ario recalcular essas fun¸c˜oes para cada

elemento. Isso gera um custo computacional maior, mas n˜ao afeta seriamente o

custo total do m´etodo que, para malhas com muitos graus de liberdade, ´e dominado

pelo custo de resolver o sistema global ﬁnal.

4.4 Resultados num´ericos

4.4.1 Elementos quadrados

Nestes testes usamos malhas uniformes com elementos quadrados. O dom´ınio

´e o quadrado unit´ario Ω = (0, 1) × (0, 1). O n´umero de onda ´e k = 100. Usamos

condi¸c˜oes de contorno de Robin

∂u

∂n

+ iku = g em ∂Ω, (4.80)

onde g ´e deﬁnida apropriadamente em cada caso para que a solu¸c˜ao exata seja a

pretendida.

Na Figura 4.6 comparamos o desempenho dos m´etodos de Galerkin, GLS,

RPPG e QSPG para solu¸c˜oes exatas que s˜ao ondas planas em v´arias dire¸c˜oes.

usada uma malha com 100×100 elementos. Vemos que, na norma do L

, a solu¸c˜ao

pelo QSPG ´e melhor que a interpolante. O RPPG tem tamb´em bom desempenho

comparado com a interpolante. O GLS exibe muita varia¸c˜ao, funcionando muito

bem para a dire¸c˜ao

mas n˜ao muito para

. Os resultados produzidos pelo m´etodo

de Galerkin s˜ao muito afetados pelo efeito de polui¸c˜ao em todas as dire¸c˜oes.

-1

-0.5

log(u − u

/ u)

Interpolante

Galerkin

GLS

RPPG

QSPG

Figura 4.6: Dependˆencia do erro em rela¸c˜ao a dire¸c˜ao θ da onda plana. Malha

com 100 × 100 elementos quadrados, k = 100, condi¸c˜oes de Robin.

Podemos observar na Figura 4.6 que a dire¸c˜ao

´e a que apresenta maiores

erros na norma L

. Na Figura 4.7 apresentamos resultados de estudos de con-

vergˆencia para uma onda plana nessa dire¸c˜ao. Para o RPPG o efeito de polui¸c˜ao

´e rapidamente reduzido pelo reﬁnamento da malha. Para o GLS isso ocorre mais

lentamente e apenas na norma do H

. Em L

o gr´aﬁco do GLS ﬁca sempre bem

distante do gr´aﬁco da interpolante, o mesmo ocorre, de forma mais acentuada, com

o m´etodo de Galerkin. Como referido anteriormente, o m´etodo de Galerkin n˜ao ´e

capaz de eliminar completamente a polui¸c˜ao na norma L

, mesmo quando h → 0.

2 2.2 2.4 2.6

−log(h)

-1.5

-1

-0.5

0.5

log (|u − u

|/|u|)

Interpolante

Galerkin

GLS

RPPG

QSPG

2 2.2 2.4 2.6

−log(h)

-2.5

-2

-1.5

-1

-0.5

0.5

log (u − u

/u)

Interpolante

Galerkin

GLS

RPPG

QSPG

Figura 4.7: Convergˆencia em H

e L

, k = 100, condi¸c˜oes de contorno de Robin,

a solu¸c˜ao exata ´e uma onda plana na dire¸c˜ao

4.4.2 Elementos retangulares

Os estudos feitos na se¸c˜ao anterior s˜ao repetidos aqui para elementos re-

tangulares. O dom´ınio ´e um quadrado unit´ario. S˜ao usadas malhas com n × 2n

elementos retangulares, ou seja, em todo elemento o comprimento mede o dobro

da altura.

Este tipo de elemento causa problemas em m´etodos que tˆem parˆametros de-

pendentes de um comprimento m´edio h em cada elemento. Quando se trata de

um elemento muito alongado, n˜ao existe uma escolha razo´avel para tal compri-

mento m´edio. No caso do RPPG, n˜ao h´a esse problema, pois n˜ao h´a parˆametros

dependentes da geometria dos elementos. J´a o QSPG lida naturalmente com essa

situa¸c˜ao, visto que os parˆametros de estabiliza¸c˜ao dependem explicitamente dos

comprimentos de todos os lados do elemento e n˜ao de um ´unico comprimento

m´edio.

-1.5

-1

-0.5

log(u − u

/ u)

Interpolante

Galerkin

GLS

RPPG

QSPG

Figura 4.8: Dependˆencia do erro em rela¸c˜ao a dire¸c˜ao θ da onda plana. Malha

com 100 × 200 elementos retangulares, k = 100, condi¸c˜oes de Robin.

Para a Figura 4.8 foi usada uma malha com 100 × 200 elementos. Observa-

se agora que a dire¸c˜ao

n˜ao ´e a de maior erro para todos os m´etodos, uma

consequˆencia do fato de que nessa dire¸c˜ao a malha ´e mais reﬁnada. Novamente

o RPPG e o QSPG exibem um desempenho muito bom, veriﬁcado novamente no

estudo de convergˆencia, nas normas L

e H

, apresentado na Figura 4.9.

2.2 2.4 2.6

−log(h)

-1.4

-1.2

-1

-0.8

-0.6

-0.4

log (|u − u

|/|u|)

Interpolante

Galerkin

GLS

RPPG

QSPG

2.2 2.4 2.6

−log(h)

-3

-2.5

-2

-1.5

-1

-0.5

log (u − u

/u)

Interpolante

Galerkin

GLS

RPPG

QSPG

Figura 4.9: Convergˆencia em H

e L

, k = 100, malhas com n × 2n elementos

retangulares, condi¸c˜oes de contorno de Robin, a solu¸c˜ao exata ´e uma onda plana

na dire¸c˜ao

4.4.3 Elementos distorcidos aleatoriamente

Nesta se¸c˜ao s˜ao realizados testes com elementos distorcidos. A distor¸c˜ao ´e ob-

tida substituindo as coordenadas (x

, y

) dos n´os interiores de uma malha formada

originalmente por elementos quadrados por coordenadas (x



, y



) aleatoriamente

perturbadas da seguinte forma:



= x

+ r

h, (4.81)



= y

+ s

h, (4.82)

onde r

e s

s˜ao sequˆencias de n´umeros aleat´orios distribu´ıdos uniformemente no

intervalo [−0.24, 0.24]. A Figura 4.10 ilustra uma malha com 20 × 20 elementos

gerada dessa forma.

Os resultados deste estudo est˜ao apresentados nas Figuras 4.11 e 4.12. Po-

demos notar que o RPPG e o QSPG produzem resultados um pouco piores do que

a interpolante, mas ainda muito melhores do que o GLS e o m´etodo de Galerkin.

Figura 4.10: Malha com 20 × 20 elementos distorcidos.

E interessante notar que neste estudo o RPPG funcionou melhor do que o QSPG.

-1

-0.8

-0.6

-0.4

-0.2

log(u − u

/ u)

Interpolante

Galerkin

GLS

RPPG

QSPG

Figura 4.11: Dependˆencia do erro em rela¸c˜ao a dire¸c˜ao θ da onda plana. Malha

com 100 × 100 elementos distorcidos, k = 100, condi¸c˜oes de Robin.

2 2.2 2.4 2.6

−log(h)

-1.5

-1

-0.5

0.5

log (|u − u

|/|u|)

Interpolante

Galerkin

GLS

RPPG

QSPG

2 2.2 2.4 2.6

−log(h)

-2.5

-2

-1.5

-1

-0.5

0.5

log (u − u

/u)

Interpolante

Galerkin

GLS

RPPG

QSPG

Figura 4.12: Convergˆencia em H

e L

, k = 100, elementos distorcidos, condi¸c˜oes

de contorno de Robin, a solu¸c˜ao exata ´e uma onda plana na dire¸c˜ao

4.4.4 Ondas n˜ao direcionais

Nesta se¸c˜ao apresentamos estudos onde as solu¸c˜oes n˜ao s˜ao ondas com di-

re¸c˜oes bem deﬁnidas. Para isso usamos as solu¸c˜oes descritas na Se¸c˜ao 2.3.3, ou

seja,

(x, y) = T





+ y



(1)





+ y



, (4.83)

onde H

(1)

´e a fun¸c˜ao de Hankel de primeiro tipo e ordem n e T

´e o polinˆomio de

Chebyshev de primeiro tipo e ordem n.

O dom´ınio Ω ´e o anel com raios interno r

e externo r

= 1:

Ω =



(x, y) ∈ R

;



+ y

< 1



. (4.84)

S˜ao usadas condi¸c˜oes de contorno de Dirichlet no contorno interno e de Robin

no contorno externo. O n´umero de onda ´e k = 100.

As malhas s˜ao geradas com n×10n elementos, onde n ´e o n´umero de divis˜oes

no sentido radial e 10n ´e o n´umero de divis˜oes no sentido angular. Uma malha

t´ıpica gerada dessa forma, com 12 × 120 elementos ´e mostrada na Figura 4.13.

Os resultados s˜ao apresentados na Figura 4.14, que corresponde a solu¸c˜ao

(4.83) com n = 0 e na Figura 4.15, onde n = 8. O m´etodo de Galerkin ´e bastante

afetado pela polui¸c˜ao, o GLS tamb´em mas em menor grau. As aproxima¸c˜oes do

QSPG s˜ao melhores do que a interpolante na norma do L

, o RPPG ´e um pouco

afetado pela polui¸c˜ao, mas para malhas mais reﬁnadas converge t˜ao bem quanto

o QSPG.

Figura 4.13: Exemplo de malha para um dom´ınio anelar.

2.2 2.4 2.6 2.8

−log(h)

-1.5

-1

-0.5

0.5

log (|u − u

|/|u|)

Interpolante

Galerkin

GLS

RPPG

QSPG

2.2 2.4 2.6 2.8

−log(h)

-2.5

-2

-1.5

-1

-0.5

0.5

log (u − u

/u)

Interpolante

Galerkin

GLS

RPPG

QSPG

Figura 4.14: Convergˆencia em H

e L

, solu¸c˜ao sem uma dire¸c˜ao deﬁnida (equa¸c˜ao

(4.83) com n = 0), k = 100, condi¸c˜oes de contorno de Dirichlet e Robin.

2.2 2.4 2.6 2.8

−log(h)

-1.5

-1

-0.5

0.5

log (|u − u

|/|u|)

Interpolante

Galerkin

GLS

RPPG

QSPG

2.2 2.4 2.6 2.8

−log(h)

-2.5

-2

-1.5

-1

-0.5

0.5

log (u − u

/u)

Interpolante

Galerkin

GLS

RPPG

QSPG

Figura 4.15: Convergˆencia em H

e L

, solu¸c˜ao sem uma dire¸c˜ao deﬁnida (equa¸c˜ao

(4.83) com n = 8), k = 100, condi¸c˜oes de contorno de Dirichlet e Robin.

4.4.5 Problema n˜ao homogˆeneo

Os m´etodos propostos neste cap´ıtulo foram constru´ıdo tendo como crit´erio a

estabiliza¸c˜ao para o problema homogˆeneo. Nesta se¸c˜ao fazemos um teste simples

para o caso nao homogˆeneo. Deﬁnimos o termo f na equa¸c˜ao

− ∆u − k

u = f (4.85)

de forma que a solu¸c˜ao exata seja dada por

u(x, y) = sen(πx) sen(πy). (4.86)

As condi¸c˜oes de contorno s˜ao de Robin. O dom´ınio ´e o quadrado unit´ario

e as malhas s˜ao as malhas distorcidas geradas exatamente como nos estudos da

Se¸c˜ao 4.4.3.

Na Figura 4.16 apresentamos um estudo de convergˆencia com k = 10. Na se-

minorma do H

todos os m´etodos produzem resultados equivalentes `a interpolante.

Na norma L

os resultados de todos os m´etodos s˜ao melhores que a interpolante,

O RPPG e O QSPG s˜ao os m´etodos que produzem melhores resultados.

A Figura 4.17 corresponde ao mesmo estudo com k = 100. Na seminorma do

novamente todos os m´etodos produzem resultados equivalentes `a interpolante.

Na norma L

notamos que as aproxima¸c˜oes do RPPG, do QSPG e de Galerkin

tˆem os erros inicialemnte pr´oximos ao erro da aproxima¸c˜ao do GLS. Mas com o

reﬁnamento da malhas o RPPG e o QSPG se afastam, produzindo aproxima¸c˜oes

melhores. O m´etodo de Galerkin tem o comportamento semelhante e ´e o que tem

os melhores resultados.

E importante observar que, como a solu¸c˜ao deste problema

depende fundamentalmente do termo de fonte, os efeitos de polui¸c˜ao n˜ao ocorrem

em nenhum desses m´etodos, inclusive no m´etodo de Galerkin.

4.4.6 Ondas evanescentes

Neste estudo consideramos a seguinte solu¸c˜ao exata do problema homogˆeneo

u(x, y) = e

−β(x sen θ−y cos θ)

iα(x cos θ+y sen θ)

, (4.87)

1 1.5 2

−log(h)

-2.2

-2

-1.8

-1.6

-1.4

-1.2

-1

log (|u − u

|/|u|)

Interpolante

Galerkin

GLS

RPPG

QSPG

1 1.5 2

−log(h)

-4.5

-4

-3.5

-3

-2.5

-2

log (u − u

/u)

Interpolante

Galerkin

GLS

RPPG

QSPG

Figura 4.16: Convergˆencia em H

e L

, problema n˜ao homogˆeneo, k = 10

2 2.2 2.4 2.6

−log(h)

-2.5

-2.4

-2.3

-2.2

-2.1

-2

log (|u − u

|/|u|)

Interpolante

Galerkin

GLS

RPPG

QSPG

2 2.2 2.4 2.6

−log(h)

-5.2

-5

-4.8

-4.6

-4.4

-4.2

-4

-3.8

log (u − u

/u)

Interpolante

Galerkin

GLS

RPPG

QSPG

Figura 4.17: Convergˆencia em H

e L

, problema n˜ao homogˆeneo, k = 100

com

k = 60, (4.88)

α = 61, (4.89)

β =

√

− k

= 11, (4.90)

para v´arias dire¸c˜oes θ entre 0 e 2π.

O dom´ınio ´e o quadrado unit´ario, discretizado por uma malha uniforme de

elementos quadrados.

Na Figura 4.18 vemos que, na norma do L

, as aproxima¸c˜oes do RPPG e

do QSPG s˜ao quase sempre melhores do que a aproxima¸c˜ao do GLS e em muitas

dire¸c˜oes melhor do que a interpolante.

3π

π 5π

3π

7π

2π

-2

-1.5

-1

-0.5

log(u − u

/ u)

Interpolante

Galerkin

GLS

RPPG

QSPG

Figura 4.18: Dependˆencia do erro na norma L

em rela¸c˜ao a dire¸c˜ao θ da onda

evanescente. Malha com 100 × 100 elementos quadrados, k = 60, condi¸c˜oes de

Robin.

4.4.7 Problema de Dirichlet - ressonˆancia num´erica

Por ﬁm, apresentamos um exemplo com condi¸c˜oes de Dirichlet onde a res-

sonˆancia num´erica deteriora seriamente as solu¸c˜oes aproximadas. O estudo ´e o

mesmo feito na Se¸c˜ao 4.4.3, mas em vez de condi¸c˜oes de Robin, consideramos

agora condi¸c˜oes de Dirichlet.

As Figuras 4.19 e 4.20 mostram os resultados. Como podemos notar, o

RPPG e o QSPG tˆem alguma vantagem sobre os outros m´etodos. No gr´aﬁco

da convergˆencia em L

e H

, podemos notar que o erro do RPPG e do QSPG

come¸cam a aumentar depois de certo reﬁnamento. Isso ocorre pois os autovalores

do problema discreto dependem de h e com o reﬁnamento da malha, um desses

autovalores se aproxima do n´umero de onda k. Se a malha continuasse sendo

reﬁnada, o erro continuaria a aumentar at´e passar essa ressonˆancia, e s´o ent˜ao

esses m´etodos voltariam a convergir.

E importante ressaltar que a ressonˆancia num´erica ´e inevit´avel para este

problema, visto que nenhum m´etodo pode eliminar totalmente a diferen¸ca entre o

n´umero de onda discreto e o exato.

-1

log(u − u

/ u)

Interpolante

Galerkin

GLS

RPPG

QSPG

Figura 4.19: Dependˆencia do erro em rela¸c˜ao a dire¸c˜ao θ da onda plana. Malha

com 100 × 100 elementos distorcidos, k = 100, condi¸c˜oes de Dirichlet.

4.5 Observa¸c˜oes

Uma limita¸c˜ao do QSPG ´e que ele ´e aplic´avel somente a elementos quadri-

laterais, que nem sempre s˜ao os mais convenientes para discretizer dom´ınios mais

complexos. Uma estabiliza¸c˜ao semelhante foi tentada para elementos triangulares,

mas os ganhos n˜ao foram t˜ao expressivos como os obtidos com o QSPG.

Apesar da relativa simplicidade das formula¸c˜oes do RPPG e do QSPG, esses

m´etodos s˜ao razoavelmente robustos em rela¸c˜ao a distor¸c˜oes dos elementos, princi-

palmente para elementos alongados, como veriﬁcamos na Se¸c˜ao 4.4.2. No entanto,

para certos tipos de elementos, como os mostrados na Se¸c˜ao 4.4.3, persiste alguma

polui¸c˜ao, que n˜ao causa tanta instabilidade para o problema de Robin, mas que

2 2.2 2.4 2.6

−log(h)

-1.5

-1

-0.5

0.5

1.5

log (|u − u

|/|u|)

Interpolante

Galerkin

GLS

RPPG

QSPG

2 2.2 2.4 2.6

−log(h)

-2

-1

log (u − u

/u)

Interpolante

Galerkin

GLS

RPPG

QSPG

Figura 4.20: Convergˆencia em H

e L

, k = 100, elementos distorcidos, condi¸c˜oes

de contorno de Dirichlet, a solu¸c˜ao exata ´e uma onda plana na dire¸c˜ao

causa s´erios efeitos causados por ressonˆancia num´ericas para o problema de Diri-

chlet. Assim, buscar m´etodos mais robustos a distor¸c˜oes continua sendo o objetivo

principal dos pr´oximos cap´ıtulos.

Cap´ıtulo 5

M´etodos de diferen¸cas ﬁnitas para

malhas n˜ao estruturadas

No QSFEM o crit´erio usado para deﬁnir os coeﬁcientes de um estˆencil de

diferen¸cas ﬁnitas para o operador de Helmholtz foi a minimiza¸c˜ao da diferen¸ca entre

o n´umero de onda exato e o discreto. A an´alise usada para chegar a esse crit´erio

requer uma malha peri´odica, o que diﬁculta a extens˜ao dos resultados obtidos para

malhas mais gerais. Al´em disso, na pr´atica, o estˆencil ´e determinado for¸cando o

n´umero de onda discreto a ser igual ao exato para ondas planas em certas dire¸c˜oes

´otimas. Isso tamb´em traz diﬁculdades quando se trata de elementos distorcidos

pois nesses casos n˜ao h´a uma escolha evidente para tais dire¸c˜oes ´otimas. Essas

diﬁculdades na escolha das dire¸c˜oes ´otimas tornam-se ainda maiores com malhas

n˜ao estruturadas onde o n´umero de coeﬁcientes no estˆencil varia de um n´o para

outro da malha.

Neste cap´ıtulo ser´a formulado um novo crit´erio para a constru¸c˜ao de um

esquema de diferen¸cas ﬁnitas. Mais precisamente, esse crit´erio consiste em mini-

mizar o erro de truncamento para ondas planas levando em conta igualmente todas

as dire¸c˜oes poss´ıveis. A formula¸c˜ao obtida ´e geral o suﬁciente para ser aplicada

a problemas em qualquer n´umero de dimens˜oes e para qualquer distribui¸c˜ao de

pontos onde a solu¸c˜ao ´e aproximada.

5.1 Nota¸c˜oes e deﬁni¸c˜oes

Apenas o problema homogˆeneo com condi¸c˜oes de contorno de Dirichlet ser´a

considerado inicialmente:

−∆u − k

u = 0 em Ω, (5.1)

u = g em Γ. (5.2)

Dado um conjunto ﬁnito X, ser´a usada a seguinte nota¸c˜ao para o n´umero de

elementos de X:

|X| = “N´umero de elementos de X”. (5.3)

Seja N = {x

, . . . , x

|N|−1

} ⊂ Ω ∪ Γ um conjunto ﬁnito de pontos onde a

solu¸c˜ao do problema (5.1)-(5.2) ser´a aproximada. O conjunto dos pontos de N

que est˜ao no interior do dom´ınio ser´a denotado por I :

I := N ∩ Ω, (5.4)

e o conjunto dos pontos de N que est˜ao na fronteira por B :

B := N ∩Γ. (5.5)

Para todo i ∈ {0, . . . , |N| − 1}, ser´a considerado como bem deﬁnido um

conjunto A

⊂ {0, . . . , |N|−1}. Tais conjuntos determinam os chamados “pontos

adjacentes ao ponto x

”. Se j ∈ A

, ent˜ao dizemos que x

´e um ponto adjacente

a x

. Os conjuntos A

podem ser deﬁnidos de v´arias formas diferentes, a ´unica

condi¸c˜ao assumida por enquanto ´e que i ∈ A

para todo i ∈ {0, . . . , |N| − 1}.

Dado um ponto x ∈ N, seu ´ındice ´e denotado por ind (x), ou seja, para todo

i ∈ {0, . . . , |N| − 1},

ind (x

) := i. (5.6)

Por extens˜ao, para um conjunto de pontos X ⊂ N, ind (X) ´e o conjunto dos

´ındices dos pontos de X. Por exemplo, ind (N) = {0, . . . , |N| − 1}.

Uma aproxima¸c˜ao por diferen¸cas ﬁnitas para o operador de Helmholtz apli-

cado a um campo u(x) ´e deﬁnida em cada ponto x

em fun¸c˜ao dos valores que u

assume nos pontos adjacentes a x

u(x

) :=



j∈A

u(x

). (5.7)

Os coeﬁcientes S

podem ser deﬁnidos de v´arias formas, resultando em diferentes

m´etodos. Tradicionalmente s˜ao usados desenvolvimentos em s´erie de Taylor para

deﬁni-los. Em geral eles dependem, para cada x

∈ N, das distˆancias entre todos

os pares de pontos adjacentes a x

A partir do operador acima deﬁnimos a aproxima¸c˜ao por diferen¸cas ﬁnitas

para a solu¸c˜ao do problema (5.1)-(5.2) como sendo a solu¸c˜ao do sistema de equa¸c˜oes

lineares



j∈A

= 0 para todo i ∈ ind (I), (5.8)

onde U

´e uma poss´ıvel aproxima¸c˜ao para a solu¸c˜ao u(x

), com U

prescritos para

∈ B e inc´ognitos para x

∈ I. Claro que esta solu¸c˜ao aproximada s´o est´a de

fato bem deﬁnida se o sistema de equa¸c˜oes acima admitir uma ´unica solu¸c˜ao, o

que depende dos dados do problema original, dos coeﬁcientes S

e da escolha dos

pontos adjacentes A

Os S

formam na verdade a matriz global associada ao sistema de equa¸c˜oes li-

neares que ´e resolvido para obter a solu¸c˜ao aproximada. Para um dado i, deﬁnimos

como sendo a lista ordenada dos coeﬁcientes S

para j ∈ A

. Em analogia aos

m´etodos de diferen¸cas ﬁnitas usuais, chamamos cada S

de um “estˆencil”, mesmo

sabendo que o sentido original da palavra se perde, aﬁnal aqui cada estˆencil pode

ser diferente dos demais.

Ao aplicarmos o operador discreto L

`a solu¸c˜ao exata u do problema obtemos

o erro de truncamento local τ

em cada ponto x

:= L

u(x

) =



j∈N

u(x

). (5.9)

Em particular, o erro de truncamento local para uma onda plana na dire¸c˜ao σ ser´a

denotado por:

(σ) := L

w(σ, x

) =



j∈N

w(σ, x

), (5.10)

onde w(σ, x) ´e a onda plana na dire¸c˜ao σ deﬁnida em (2.30):

w(σ, x) = e

ikσ·x

. (5.11)

Exemplo 5.1 O m´etodo de diferen¸cas ﬁnitas cl´assico, baseado em aproxima¸c˜oes

por s´eries de Taylor, para o problema unidimensional em Ω = (0, 1) pode ser

formulado usando as deﬁni¸c˜oes acima com

N = {0, h, . . . , 1 − h, 1}, (5.12)

= ih, (5.13)

= {i − 1, i, i + 1} ∩ ind (N) , (5.14)

B = {0, 1}, (5.15)

I = N  B (5.16)

e, para todo i ∈ ind (I),

i−1

, (5.17)

= k

−

, (5.18)

i+1

. (5.19)

Exemplo 5.2 No caso bidimensional, com Ω = (0, 1) × (0, 1), um conjunto de

(m + 1) × (n + 1) pontos uniformemente distribu´ıdos pode ser constru´ıdo como

N = {(j∆x, l∆y); j ∈ {0, . . . , m}, l ∈ {0, . . . , n}}, (5.20)

onde ∆x = 1/m e ∆y = 1/n. Para um m´etodo de diferen¸cas ﬁnitas com estˆenceis

de 9 pontos, os pontos adjacentes seriam

= ind ({x

+ (j∆x, l∆y); j, l ∈ {−1, 0, 1}}) ∩ ind (N) . (5.21)

Adiante, o m´etodo QSFEM ser´a formulado usando os conjuntos A

acima.

O principal resultado deste cap´ıtulo ´e uma forma sistem´atica para deﬁnir os

coeﬁcientes S

correspondentes a um conjunto arbitr´ario de pontos N e conjun-

tos de pontos adjacentes A

tamb´em arbitr´arios. Isso ser´a descrito na Se¸c˜ao 5.3.

Na pr´oxima se¸c˜ao ser´a descrita uma primeira tentativa de atingir esse objetivo,

inspirada nos resultados do QSFEM.

5.2 Truncamento nulo em algumas dire¸c˜oes

Dado um ponto interior x

, ´e poss´ıvel obter a lista de coeﬁcientes S

, tal que

(σ

) = 0 para |A

| − 1 diferentes dire¸c˜oes σ

. Para isso basta resolver o sistema

de equa¸c˜oes



j∈A

w(σ

, x

) = 0, s = 1, 2, . . . , n

− 1, (5.22)

com o coeﬁciente predeﬁnido

= 1. (5.23)

Assim, ﬁxando dire¸c˜oes σ

, σ

, . . ., σ

, com M = min

i∈ind(I)

(|A

|) − 1, ´e poss´ıvel

obter um esquema de diferen¸cas tal que τ

(σ

) = 0 para todo x

∈ I e para

todas as dire¸c˜oes σ

, s = 1, . . . , M, o que, consequentemente, produz uma solu¸c˜ao

aproximada nodalmente exata se a solu¸c˜ao exata for uma onda plana em uma das

dire¸c˜oes escolhidas σ

Para o problema unidimensional, qualquer solu¸c˜ao ´unica do problema ho-

mogˆeneo pode ser escrita como uma combina¸c˜ao linear das ondas planas w

(x) =

w(1, x) = exp(ikx) e w

(x) = w(−1, x) = exp(−ikx). Assim, se escolhermos S

que resultem em aproxima¸c˜oes nodalmente exatas para w

and w

, ent˜ao obtemos

solu¸c˜oes nodalmente exatas para qualquer solu¸c˜ao do problema homogˆeneo. Por

exemplo, para Ω = (0, 1), N = {x

, . . . , x

}, com x

< x

i+1

e A

= {i − 1, i, i +

1} ∩ {0, . . . , n}, obtemos tais coeﬁcientes ´otimos resolvendo as equa¸c˜oes (5.22) e

(5.23) com σ

= 1 e σ

= −1:

i−1

) + S

i+1

i−1

) = 0 (5.24)

i−1

) + S

i+1

i−1

) = 0 (5.25)

= 1 (5.26)

cuja solu¸c˜ao ´e

i−1

= −

sen k(x

i+1

− x

)

sen k(x

i+1

− x

i−1

)

, (5.27)

= 1, (5.28)

i+i

= −

sen k(x

− x

i−1

)

sen k(x

i+1

− x

i−1

)

. (5.29)

No caso bidimensional, para Ω = (0, 1) × (0, 1), usando os conjuntos do

Exemplo 5.2 com m = n e ∆x = ∆y = h, o m´etodo QSFEM, discutido no Cap´ıtulo

3, pode ser obtido usando o esquema descrito acima com as dire¸c˜oes ﬁxadas

σ = (cos θ, sen θ), (5.30)

para os oito valores de θ no conjunto



+ n

3π

+ n

; n = 0, 1, 2, 3



. (5.31)

Resolvendo as equa¸c˜oes (5.22) e (5.23) para essas dire¸c˜oes encontramos

= 1, (5.32)

ind(x

+(h,0))

= S

ind(x

+(−h,0))

= S

ind(x

+(0,−h))

= S

ind(x

+(0,h))

= A, (5.33)

ind(x

+(h,h))

= S

ind(x

+(−h,h))

= S

ind(x

+(−h,−h))

= S

ind(x

+(h,−h))

= B, (5.34)

onde

A =

− c

+ s

) − c

+ s

)

, (5.35)

B =

+ s

− c

− s

+ s

) − c

+ s

)

, (5.36)

= cos



kh cos



, (5.37)

= cos



kh sen



, (5.38)

= cos



kh cos

3π



, (5.39)

= cos



kh sen

3π



. (5.40)

O estˆencil assim obtido ´e idˆentico ao do QSFEM, diferindo apenas pelo fator

uma vez que no QSFEM foi ﬁxado S

= 4.

Observamos que o QSFEM ´e restrito a estˆenceis de 9 pontos e malhas unifor-

mes com elementos quadrados. A alternativa apresentada acima ´e, em princ´ıpio,

aplic´avel a conﬁgura¸c˜oes mais gerais, inclu´ındo malhas n˜ao uniformes e n˜ao es-

truturadas. Entretanto, para uma discretiza¸c˜ao mais geral, uma escolha para as

dire¸c˜oes ´otimas ou quase ´otimas n˜ao ´e evidente, sendo uma quest˜ao deixada em

aberto.

Uma poss´ıvel alternativa ao esquema descrito acima seria considerar mais

dire¸c˜oes em cada ponto do que as |A

| − 1 propostas, permitindo assim uma dis-

tribui¸c˜ao mais densa e uniforme das dire¸c˜oes privilegiadas. Mas isso tornaria o

sistema (5.22) mal posto, visto que ter´ıamos mais equa¸c˜oes do que inc´ognitas.

Insistindo nesta id´eia, poder´ıamos substituir as equa¸c˜oes exatas (5.22) por um

problema de m´ınimos quadrados, ou seja, encontrar, para cada i ∈ ind (I), os

coeﬁcientes S

que minimizam

J(S

) :=

2π



s=1

|τ

(σ

(5.41)

2π



s=1









j∈A

w(σ

, x

)





k∈A

w(σ

, x

)







(5.42)



j∈A



k∈A





s=1

2π

w(σ

, x

− x

)



, (5.43)

com S

prescrito. O fator constante

2π

foi inclu´ıdo por conveniˆencia, ele n˜ao

interfere na minimiza¸c˜ao. Deﬁnindo os coeﬁcientes dessa forma, o n´umero n de

dire¸c˜oes se torna arbitr´ario, restando na pr´atica a quest˜ao de como distribuir essas

dire¸c˜oes. Uma possibilidade seria distribu´ı-las uniformemente, ou seja, fazer σ



cos

2π(s−1)

n−1

, sen

2π(s−1)

n−1



. Neste caso a express˜ao (5.43) converge para um valor que

depende apenas de x

− x

quando n aumenta, j´a que

lim

n→∞

2π



s=1

w(σ

, x

− x

) =



2π



(cos θ, sen θ), x

− x



dθ. (5.44)

Devido ao grande esfor¸co computacional de se incluir muitas dire¸c˜oes no pro-

blema de minimiza¸c˜ao proposto acima, essa id´eia n˜ao foi levada adiante. Mas a

situa¸c˜ao limite, com inﬁnitas dire¸c˜oes distribu´ıdas uniformemente, ser´a desenvol-

vida e estendida para d dimens˜oes na pr´oxima se¸c˜ao.

5.3 Minimiza¸c˜ao do erro de truncamento

Nesta se¸c˜ao propomos um outro crit´erio para determinar S

: a minimiza¸c˜ao,

com respeito a S

, do erro de truncamento local τ

(σ) na norma L

(Σ). O conjunto

Σ ´e dado por

Σ :=



σ ∈ R

: |σ| = 1



, (5.45)

ou seja, ´e o conjunto de todas as dire¸c˜oes σ poss´ıves em d dimens˜oes (o c´ırculo

unit´ario em duas dimens˜oes ou a esfera unit´aria em trˆes dimens˜oes).

Mais explicitamente, para cada n´o i deﬁnimos o funcional de m´ınimos qua-

drados

J(S

) = τ

(σ)

(Σ)



|τ

(σ)|

dσ, (5.46)

ent˜ao predeﬁnimos S

= 1 e escolhemos S

para j ∈ A

 {i} que minimizam

J(S

). O diferencial dσ est´a apenas indicando que a vari´avel do integrando ´e σ.

A integral ´e calculada usando a medida de arco no c´ırculo unit´ario para d = 2, a

medida de ´area na esfera unit´aria para d = 3 ou seus an´alogos em mais dimens˜oes,

de modo que cada dire¸c˜ao tenha o mesmo peso na integra¸c˜ao.

Substituindo (5.10) em (5.46),

J(S

) =





m∈A



n∈A

w(σ, x

)w(σ, x

)dσ (5.47)



m∈A



n∈A



w(σ, x

)w(σ, x

)dσ (5.48)



m∈A



n∈A



ikσ·(x

−x

)

dσ (5.49)



m∈A



n∈A

, (5.50)

onde



ikσ·(x

−x

)

dσ (5.51)

= (2π)

−1

(kh

)

(kh

)

−1

. (5.52)

Minimizando (5.50) com a restri¸c˜ao S

= 1, concluimos que, para cada x

∈ I, S

s˜ao determinados resolvendo o sistema de |A

|−1 equa¸c˜oes com |A

|−1 inc´ognitas



n∈A

= 0 ∀m ∈ A

 {i}. (5.53)

Adiante vamos nos referir ao sistema de equa¸c˜oes acima pela express˜ao compacta

W S

= 0.

Chamamos o m´etodo assim obtido de “Quasi Optimal Finite Diﬀerence

method”, ou pela sigla QOFD.

Antes de resolver o sistema de equa¸c˜oes (5.53) ´e necess´ario calcular a integral

(5.51). Isso ´e feito atrav´es da f´ormula expl´ıcita (5.52), deduzida como segue.

Escrevendo x

− x

= h

com |r

| = 1 e h

= |x

− x

|, ent˜ao



ikh

·σ

dσ



cos(kh

· σ)dσ + i



sen(kh

· σ)dσ.

(5.54)

Uma integral da forma



f(r · σ)dσ (5.55)

pode ser simpliﬁcada usando a mudan¸ca de vari´aveis σ = Qu, onde Q ´e uma

rota¸c˜ao tal que r = Q(1, 0, . . . , 0). Assim,

r · σ = r · Qu = Q

r · u = Q

−1

r · u = (1, 0, . . . , 0) · u, (5.56)

Q(Σ) = Σ, (5.57)

dσ = du, (5.58)

e (5.55) se torna



f((1, 0, . . . , 0) · u)du, (5.59)

que n˜ao depende de r.

Para d = 2 podemos usar coordenadas polares para obter:



f((1, 0, . . . , 0) · u)du =



2π

f(cos(θ))dθ, (5.60)

Para d ≥ 3, aplicando a transforma¸c˜ao de coordenadas esf´ericas para d di-

mens˜oes

= cos(θ

)

= sen(θ

) cos(θ

)

= sen(θ

) sen(θ

) cos(θ

)

···

n−1

= sen(θ

) ···sen(θ

n−2

) cos(θ

n−1

)

= sen(θ

) ···sen(θ

n−2

) sen(θ

n−1

)

du = sen

n−2

(θ

) sen

n−3

(θ

) ···sen(θ

n−1

)dθ

···dθ

n−1

obtemos



f((1, 0, . . . , 0) · u)du



2π

···



f(cos θ

) sen

n−2

(θ

) sen

n−3

(θ

) ···sen(θ

n−1

)dθ

dθ

···dθ

n−1





f(cos θ

) sen

n−2

(θ

)dθ

 



sen

n−3

(θ

)dθ



···





2π

sen(θ

n−1

)dθ

n−1



Para f(s) = cos(kh

s) ou f(s) = sen(kh

s) as integrais acima podem ser

calculadas usando o software Mathematica. Os resultado s˜ao:



cos(kh

· u)du = (2π)

−1

(kh

)

(kh

)

−1

e (5.61)



sen(kh

· u)du = 0. (5.62)

´e a fun¸c˜ao de Bessel de primeiro tipo e ordem s. Com esse resultado conclu´ımos

que W ´e uma matriz real, logo os coeﬁcientes S

tamb´em o s˜ao.

Para d = 2,

= 2πJ

(kh

) (5.63)

= 2π



1 −

(kh

)

(kh

)

−

(kh

)

2304

+ . . .



, (5.64)

e para d = 3,

= (2π)

(kh

)

√

(5.65)

= 4π

sen kh

(5.66)

= 4π



1 −

(kh

)

(kh

)

120

−

(kh

)

5040

+ . . .



. (5.67)

5.4 Caso unidimensional

A integral d-dimensional (5.46) n˜ao faz sentido para o problema unidimensi-

onal. Mas ´e interessante notar que o resultado geral ﬁnal (5.61) pode ser aplicado

para d = 1, resultando em

= (2π)

−

(kh

)

(kh

)

−

(5.68)

= 2 cos(k(x

− x

)). (5.69)

Usando as equa¸c˜oes gerais (5.53) e considerando novamente uma discretiza¸c˜ao

unidimensional usual, obtemos os seguintes coeﬁcientes:

i−1

= −

sen k(x

i+1

− x

)

sen k(x

i+1

− x

i−1

)

, (5.70)

= 1, (5.71)

i+1

= −

sen k(x

− x

i−1

)

sen k(x

i+1

− x

i−1

)

. (5.72)

Este ´e exatamente o resultado em (5.27)-(5.29), ou seja, os coeﬁcientes acima

produzem solu¸c˜oes nodalmente exatas para o problema homogˆeneo, mesmo para

malhas n˜ao uniformes.

A matriz global resultante n˜ao ´e sim´etrica em geral. A simetria ´e recuperada

para malhas uniformes, com x

i+1

− x

= h. Neste caso os coeﬁcientes (5.70) e

(5.72) simpliﬁcam para

i−1

= −

2 cos(kh)

, (5.73)

i+1

= −

2 cos(kh)

. (5.74)

5.5 Caso bidimensional

5.5.1 Malhas n˜ao estruturadas

Dada uma malha de elementos ﬁnitos, obtemos facilmente os conjuntos N,

B e I deﬁnidos no in´ıcio deste cap´ıtulo. Outra informa¸c˜ao de interesse s˜ao os con-

juntos E

dos elementos a que cada ponto nodal x

pertence. Para as necessidades

deste cap´ıtulo, podemos considerar que um elemento e ´e o conjunto do pontos

nodais de e. Sendo assim, podemos construir os conjuntos A

da seguinte forma:

= ind





e∈E



, (5.75)

ou seja, A

´e o conjunto dos ´ındices dos n´os pertencentes a pelo menos um elemento

que cont´em x

Uma outra possibilidade que ser´a testada adiante ´e aumentar os conjuntos

. Por exemplo, se chamarmos os pontos com ´ındice em A

de vizinhos de x

ent˜ao podemos deﬁnir A



como o conjunto de ´ındices dos vizinhos de x

mais os

vizinhos dos vizinhos de x

. J´a A



poderia incluir tamb´em os vizinhos dos vizinhos

dos vizinhos. Seguindo essa id´eia, deﬁnimos recursivamente os seguintes conjuntos:

= ind





e∈E



, (5.76)

n+1



j∈A

. (5.77)

A Figura 5.1 ilustra esse processo de adicionar recursivamente camadas de pontos

nodais at´e chegar ao conjunto A

Figura 5.1: 4 camadas de pontos adjacentes a um dado n´o.

Observa¸c˜ao 5.1 O termo “malhas n˜ao estruturadas” que aparece no t´ıtulo deste

cap´ıtulo ´e devido somente a esta forma particular de construir os conjuntos A

que ´e v´alida para qualquer tipo de malha. Poder´ıamos tamb´em deﬁnir os pontos

adjacentes sem o uso de malhas, mas essa possibilidade n˜ao ser´a explorada aqui.

Observa¸c˜ao 5.2 No contexto dos m´etodos de elementos ﬁnitos, temos que

a(φ

, φ

) = 0 ∀ j /∈ A

, (5.78)

onde φ

s˜ao as fun¸c˜oes bases nodais. Ou seja, se adotarmos A

= A

, ent˜ao os

coeﬁcientes S

para o m´etodo QOFD deﬁnem uma matriz esparsa com o mesmo

perﬁl da matriz global do m´etodo de Galerkin com elementos lineares ou bilineares.

5.5.2 Restri¸c˜ao a malhas uniformes com elementos quadrados

Quando todos os elementos s˜ao quadrados iguais com lados medindo h, con-

siderando

= A

, (5.79)

ent˜ao obtemos um m´etodo de diferen¸cas ﬁnitas com um estˆencil de 9 pontos que

se repete em cada um dos pontos x

, ou seja,

ind(x

+(mh,nh))

= S

ind(x

+(mh,nh))

(5.80)

para todos i, j ∈ I e m, n ∈ {−1, 0, 1}. Com ajuda do Mathematica ´e poss´ıvel

calcular esses coeﬁcientes explicitamente, mas o resultado ´e longo demais para ser

escrito ou usado em computa¸c˜oes. Mais interessante ´e usar esses coeﬁcientes para

calcular a rela¸c˜ao entre o n´umero de onda discreto e o exato, o que permite uma

compara¸c˜ao entre o QOFD e o QSFEM. Procedendo como na Se¸c˜ao 3.3.2, obtemos

para o QOFD:

− k



5 + cos 2θ



cos 8θ

774144

(kh)

+ O((kh)

). (5.81)

Relembrando que no QSFEM

− k

cos 8θ

774144

(kh)

+ O((kh)

), (5.82)

conclu´ımos que em malhas uniformes o m´etodo QOFD ´e diferente do QSFEM, mas

o erro de fase ´e da mesma ordem, sendo multiplicado por um fator



5+cos 2θ



que

varia entre 1 e

. Na Figura 5.2 o m´etodo QOFD ´e comparado com o QSFEM,

n˜ao ´e poss´ıvel distinguir os resutados da interpolande, do QSFEM e do QOFD.

5.6 Implementa¸c˜ao computacional

A implementa¸c˜ao do QOFD em malhas de elementos ﬁnitos requer um pre-

processamento da malha para identiﬁcar os conjuntos dos pontos adjacentes A

de modo que essa informa¸c˜ao possa ser recuperada de forma eﬁciente quando ne-

cess´ario.

Para isso ´e constru´ıda, para cada ponto nodal N, uma lista indicando os

elementos a que N pertence. Esse processo ´e bastante simples, basta percorrer

-1.07

-1.065

-1.06

-1.055

-1.05

log(u − u

/ u)

Interpolante

QOFD

QSFEM

Figura 5.2: Dependˆencia do erro em rela¸c˜ao a θ, k = 150, 150 × 150 elementos

quadrados, malha uniforme, elementos quadrados, condi¸c˜oes de Dirichlet.

elemento por elemento, incluindo cada elemento E na lista de elementos que cont´em

N para todos os n´os N de E.

Tendo assim constru´ıda essa lista de elementos a que cada n´o pertence, torna-

se f´acil identiﬁcar os pontos adjacentes A

de um dado ponto x

sempre que ne-

cess´ario. Basta visitar todos os elementos E a que x

pertence e incluir os n´os de

E em A

Sabendo construir A

para todo i, ent˜ao A

pode ser constru´ıdo a partir dos

conjuntos A

, para j ∈ A

e assim recursivamente para A

, etc.

Na pr´atica a mem´oria e o tempo de processamento extras necess´arios para o

processo descrito acima n˜ao causa um impacto signiﬁcativo na aplica¸c˜ao do QOFD.

Muito mais cr´ıtico ´e o c´alculo das matrizes locais W e a solu¸c˜ao dos sistemas

lineares locais W S

= 0 para a determina¸c˜ao dos estˆenceis S

para todos os n´os x

As matrizes locais W s˜ao em geral mal condicionadas e a situa¸c˜ao piora

quando kh → 0 ou quando o n´umero de pontos adjacentes aumenta. Por isso,

para resolver os sistemas locais W S

= 0 quando o n´umero de pontos adjacentes a

´e grande foi necess´ario usar aritm´etica de ponto ﬂutuante com precis˜ao maior

do que a normal do hardware utilizado, o que aumenta signiﬁcativamente o custo

computacional do QOFD. Devemos notar no entanto que esses problemas locais

s˜ao independentes uns dos outros, podendo ser resolvidos paralelamente de forma

eﬁcientemente em sistemas com mais de um processador. Al´em disso, o custo

computacional da solu¸c˜ao do sistema global ﬁnal ´e de ordem mais alta em rela¸c˜ao

ao n´umero de pontos nodais do que o custo total de resolver os sistemas locais.

Assim, pelo menos em teoria, para problemas muito grandes, o custo ﬁnal acabaria

sendo dominado pela solu¸c˜ao do sistema de equa¸c˜oes global.

5.7 Resultados num´ericos

Apresentamos nesta se¸c˜ao experimentos num´ericos com o m´etodo QOFD. O

erro da solu¸c˜ao aproximada ser´a medido de quatro formas. Na norma L

e na

seminorma H

, o erro medido corresponde ao erro de uma solu¸c˜ao aproximada

dada por

(x) =



i ∈ ind(N)

(x), (5.83)

onde U

´e a solu¸c˜ao aproximada pontual obtida com o QOFD e φ

s˜ao as fun¸c˜oes

bases nodais do espa¸co de elementos ﬁnitos lagrangianos lineares ou bilineares S

As outras duas normas usadas s˜ao as normas discretas l

e l

∞

, os erros nessas

normas s˜ao deﬁnidos como

u − u





|N|



i ∈ ind(N)

− u(x

, (5.84)

u − u



∞

= max

i ∈ ind(N)

− u(x

)|, (5.85)

onde u ´e a solu¸c˜ao exata.

5.7.1 Elementos quadrados

Neste estudo veriﬁcamos a convergˆencia do m´etodo QOFD com elementos

quadrados.

Os seguintes dados foram usados:

• Dom´ınio: quadrado unit´ario (0, 1) × (0, 1).

• Solu¸c˜ao exata: onda plana na dire¸c˜ao θ = 20

◦

, com k = 50.

• Malhas estruturadas com n × n elementos quadrados, com n variando de

50 at´e 800.

Na Figura 5.3 n˜ao se nota diferen¸cas entre o m´etodo QOFD e a interpolante,

medindo o erro nas normas L

e H

2 2.5

−log(h)

-3.5

-3

-2.5

-2

-1.5

log (u − u

/u)

Interpolante

QOFD

2 2.5

−log(h)

-1.8

-1.6

-1.4

-1.2

-1

-0.8

-0.6

log (|u − u

|/|u|)

Interpolante

QOFD

Figura 5.3: Convergˆencia em L

e H

, onda plana na dire¸c˜ao θ = 20

◦

, com k = 50,

elementos quadrados.

Na Figura 5.4 a convergˆencia ´e veriﬁcada nas normas discretas l

e l

∞

. Nessas

normas ´e poss´ıvel observar o efeito do mal condicionamento das matrizes W sobre

o comportamento assint´otico das aproxima¸c˜oes. O m´etodo QOFD 16 ´e implemen-

tado usando precis˜ao de ponto ﬂutuante de aproximadamente 16 casas decimais na

solu¸c˜ao dos sistemas de equa¸c˜oes locais W S = 0. Podemos observar instabilidade

a partir de kh ≈ 0.18. Aumentando a precis˜ao para 19 casas decimais, a instabi-

lidade ocorre para malhas mais reﬁnadas. Com 32 casas decimais n˜ao observamos

instabilidade. A convergˆencia observada ´e de sexta ordem.

Aumentando o n´umero de onda para k = 200, convergˆencia em L

e H

ainda ocorre normalmente conforme vemos na Figura 5.5.

2 2.5

−log(h)

-11

-10

-9

-8

-7

-6

-5

-4

log (|u − u

)

QOFD 16

QOFD 19

QOFD 32

2 2.5

−log(h)

-11

-10

-9

-8

-7

-6

-5

-4

log (|u − u

∞

)

QOFD 16

QOFD 19

QOFD 32

Figura 5.4: Convergˆencia em l

e l

∞

, onda plana na dire¸c˜ao θ = 20

◦

, k = 50,

elementos quadrados.

2.4 2.6 2.8

−log(h)

-2.2

-2

-1.8

-1.6

-1.4

-1.2

log (u − u

/u)

Interpolante

QOFD

2.4 2.6 2.8

−log(h)

-1.2

-1.1

-1

-0.9

-0.8

-0.7

-0.6

log (|u − u

|/|u|)

Interpolante

QOFD

Figura 5.5: Convergˆencia em L

e H

, onda plana na dire¸c˜ao θ = 20

◦

, k = 200,

elementos quadrados.

5.7.2 Elementos distorcidos

Os experimentos da se¸c˜ao anterior s˜ao repetidos aqui com elementos distor-

cidos. A distor¸c˜ao ´e obtida substituindo as coordenadas (x

, y

) dos n´os interiores

por coordenadas (x



, y



) aleatoriamente perturbadas da seguinte forma:



= x

+ r

h, (5.86)



= y

+ s

h, (5.87)

onde r

e s

s˜ao sequˆencias de n´umeros aleat´orios distribu´ıdos uniformemente no

intervalo [−0.2, 0.2].

Na Figura 5.6 n˜ao se nota diferen¸cas entre o m´etodo QOFD e a interpolante.

2 2.5

−log(h)

-3

-2.5

-2

-1.5

-1

log (u − u

/u)

Interpolante

QOFD

2 2.5

−log(h)

-1.6

-1.4

-1.2

-1

-0.8

-0.6

log (|u − u

|/|u|)

Interpolante

QOFD

Figura 5.6: Convergˆencia em L

e H

, onda plana na dire¸c˜ao θ = 20

◦

, k = 50,

elementos distorcidos.

Nas normas discretas a ordem de convergˆencia diminui de 6 no caso de ele-

mentos quadrados para 3 no caso dos elementos distorcidos testados nesta se¸c˜ao,

conforme observamos na Figura 5.7. Como o erro n˜ao atinge valores t˜ao pequenos

como no estudo anterior, n˜ao ´e poss´ıvel observar o efeito do mal condicionamento

da matriz W .

Aumentando o n´umero de onda para k = 200, a partir de kh ≈ 0.6, n˜ao h´a

diferen¸cas signﬁcativas entre a convergˆencia em L

da interpolante e a convergˆencia

da solu¸c˜ao pelo m´etodo QOFD (Fig. 5.8).

2 2.5

−log(h)

-5

-4.5

-4

-3.5

-3

-2.5

-2

log (|u − u

)

QOFD 16

QOFD 19

QOFD 32

2 2.5

−log(h)

-4.5

-4

-3.5

-3

-2.5

-2

-1.5

log (|u − u

∞

)

QOFD 16

QOFD 19

QOFD 32

Figura 5.7: Convergˆencia em l

e l

∞

, onda plana na dire¸c˜ao θ = 20

◦

, k = 50,

elementos distorcidos.

2.4 2.6 2.8

−log(h)

-2

-1.5

-1

-0.5

log (u − u

/u)

Interpolante

QOFD

2.4 2.6 2.8

−log(h)

-1

-0.8

-0.6

-0.4

-0.2

log (|u − u

|/|u|)

Interpolante

QOFD

Figura 5.8: Convergˆencia em L

e H

, onda plana na dire¸c˜ao θ = 20

◦

, k = 50,

elementos distorcidos.

5.7.3 Malhas n˜ao estruturadas

Neste estudo o m´etodo QOFD ´e testado em malhas n˜ao estruturadas. O

dom´ınio considerado ´e um c´ırculo centrado na origem com diˆametro igual a 1.

Nos testes de convergˆencia v´arias malhas diferentes s˜ao geradas usando o software

triangle (Shewchuk (1996)), com a ´area m´axima dos elementos cada vez menor. O

parˆametro h que aparece nos gr´aﬁcos de convergˆencia ´e deﬁnido como

h =

√

N´umero de elementos

. (5.88)

A Figura 5.9 ilustra uma das malhas gerada pelo algoritmo usado.

Figura 5.9: Malha n˜ao estruturada.

O m´etodo QOFD ´e testado usando os conjuntos A

, A

e A

deﬁnidos na

Se¸c˜ao 5.5.1, nos gr´aﬁcos as trˆes possibilidades correspondem respectivamente `as

etiquetas QOFD 1L, QOFD 2L e QOFD 3L.

Os gr´aﬁcos na Figura 5.10 correspondem a um teste de convergˆencia com

k = 50 com a solu¸c˜ao sendo uma onda plana na dire¸c˜ao θ = 20

◦

. Em H

convergˆencia ´e igual `a interpolante nos trˆes casos. Em L

, usando apenas uma

camada de elementos o m´etodo QOFD n˜ao converge t˜ao bem quanto a interpolante,

mas com 2 ou 3 camadas a diferen¸ca desaparece.

Nas normas discretas l

e l

∞

(Fig. 5.11) notamos que o m´etodo QOFD tem

um desempenho muito melhor do que o m´etodo de Galerkin. Mas pelos gr´aﬁcos

2 2.2 2.4 2.6

−log(h)

-2.5

-2

-1.5

-1

-0.5

log (u − u

/u)

Interpolante

Galerkin

QOFD 1L

QOFD 2L

QOFD 3L

2 2.2 2.4 2.6

−log(h)

-1

-0.5

log (|u − u

|/|u|)

Interpolante

Galerkin

QOFD 1L

QOFD 2L

QOFD 3L

Figura 5.10: Convergˆencia em L

e H

, k = 50, malha n˜ao estruturada.

n˜ao ´e poss´ıvel estimar uma ordem de convergˆencia, visto que o erro n˜ao decai

monotonicamente com o reﬁnamento da malha. Acreditamos que isso se deve

principalmente ao mal condicionamento das matrizes W , mas ´e poss´ıvel tamb´em

que haja ressonˆancia num´erica. Esses gr´aﬁcos servem tamb´em para real¸car o efeito

do aumento do n´umero de camadas de elementos. Com apenas uma camada, o erro

do m´etodo QOFD ´e pelo menos uma ordem de grandeza menor do que o do m´etodo

de Galerkin. Aumentando para duas camadas, obtemos erros da ordem de 10

−6

para as malhas mais reﬁnadas. Para trˆes camadas o erro chega a ser t˜ao pequeno

quanto 10

−11

, mas n˜ao diminui de forma est´avel quando a malha ´e reﬁnada.

5.7.4 Problema de Dirichlet

Por ﬁm, repetimos aqui um exemplo com condi¸c˜oes de Dirichlet feito para

os m´etodos RPPG e QSPG na Se¸c˜ao 4.4.7.

As Figuras 5.12 e 5.13 mostram os resultados. O m´etodo QOFD sofre bem

menos com ressonˆancia num´erica do que o QSPG. Essa maior robustez em rela¸c˜ao

a ressonˆancia num´erica mesmo com o tipo de distor¸c˜ao da malha considerado neste

essemplo foi observada tamb´em em outros experimentos. Isso indica que o n´umero

2 2.2 2.4 2.6

−log(h)

-12

-10

-8

-6

-4

-2

log (|u − u

)

Galerkin

QOFD 1L

QOFD 2L

QOFD 3L

2 2.2 2.4 2.6

−log(h)

-12

-10

-8

-6

-4

-2

log (|u − u

∞

)

Galerkin

QOFD 1L

QOFD 2L

QOFD 3L

Figura 5.11: Convergˆencia em l

e l

∞

, k = 50, malha n˜ao estruturada.

de onda discreto do m´etodo QOFD ´e mais preciso do que o do QSPG para esses

tipos de malha.

-1

-0.5

0.5

1.5

log(u − u

/ u)

Interpolante

QOFD

QSPG

Figura 5.12: Dependˆencia do erro em rela¸c˜ao a dire¸c˜ao θ da onda plana. Malha

com 100 × 100 elementos distorcidos, k = 100, condi¸c˜oes de Dirichlet.

2 2.2 2.4 2.6

−log(h)

-1.5

-1

-0.5

0.5

1.5

log (|u − u

|/|u|)

Interpolante

QOFD

QSPG

2 2.2 2.4 2.6

−log(h)

-2

-1

log (u − u

/u)

Interpolante

QOFD

QSPG

Figura 5.13: Convergˆencia em H

e L

, k = 100, elementos distorcidos, condi¸c˜oes

de contorno de Dirichlet, a solu¸c˜ao exata ´e uma onda plana na dire¸c˜ao

5.8 Observa¸c˜oes

A an´alise num´erica da formula¸c˜ao de diferen¸cas ﬁnitas proposta (QOFD) em

malhas n˜ao estruturadas, em princ´ıpio, poderia ser desenvolvida usando o teorema

de equivalˆencia de Lax. Nesse sentido, pelo menos para malhas uniformes, pode-

mos estimar o erro de truncamento. Resta ainda provar a estabilidade, que normal-

mente ´e a quest˜ao mais dif´ıcil. Independentemente da an´alise num´erica, os estudos

de convergˆencia apresentados atestam um comportamento bastante robusto desta

formula¸c˜ao a distor¸c˜oes da malha, mesmo com valores bastante elevados do n´umero

de onda k.

O tratamento do problema n˜ao homogˆeneo e de condi¸c˜oes de contorno n˜ao

essenciais com as formula¸c˜oes deste cap´ıtulo n˜ao ´e uma extens˜ao imediata do que

est´a feito. Essas quest˜oes seriam em princ´ıpio muito mais simples se os m´etodos

apresentados fossem baseados em uma formula¸c˜ao variacional. O pr´oximo cap´ıtulo

trata exatamente disso. Ser´a apresentado um m´etodo de elementos ﬁnitos de

Petrov-Galerkin que ´e contru´ıdo satisfazendo o mesmo crit´erio usado para de-

terminar os esquemas de diferen¸cas ﬁnitas deste cap´ıtulo.

Cap´ıtulo 6

M´etodos de Petrov-Galerkin com

fun¸c˜oes peso quase ´otimas

Conforme j´a adiantado no ﬁnal do cap´ıtulo anterior, o objetivo deste cap´ıtulo

´e aplicar o crit´erio apresentado na Se¸c˜ao 5.3 na formula¸c˜ao de um m´etodo de ele-

mentos ﬁnitos de Petrov-Galerkin variacionalmente consistente e com propriedades

de estabilidade e precis˜ao equivalentes ao m´etodo de diferen¸cas ﬁnitas QOFD.

Uma primeira id´eia seria adaptar o m´etodo QSPG do Cap´ıtulo 4. Mas,

conforme discutido na Se¸c˜ao 4.2.1, para assegurar a conformidade do QSPG ´e

necess´ario ajustar seus parˆametros para elementos quadrados, o que ´e justamente

a restri¸c˜ao que estamos tentando evitar. Por esse motivo os m´etodos apresentados

neste cap´ıtulo s˜ao constru´ıdos de forma diferente dos m´etodos do Cap´ıtulo 4.

Para facilitar as referˆencias ao cap´ıtulo anterior, ser´a considerado inicial-

mente o problema de Helmholtz sem termo de fonte e com condi¸c˜oes de contorno

de Dirichlet. Al´em disso ser˜ao usadas as mesmas nota¸c˜oes para os conjuntos dos

pontos nodais, pontos interiores e de contorno: respectivamente N, I e B.

O espa¸co U

, onde a solu¸c˜ao aproximada ´e buscada, ´e espa¸co de elementos

ﬁnitos lagrangianos C

bilineares em 2D. Novamente denotamos por φ

as fun¸c˜oes

bases nodais do espa¸co U

O espa¸co V

´e gerado a partir de fun¸c˜oes base nodais modiﬁcadas ψ

= span{ψ

; i ∈ I}. (6.1)

Cada ψ

ser´a constru´ıda para ter o mesmo suporte da fun¸c˜ao φ

correspon-

dente, o que garante que a matriz global resultante tenha o mesmo padr˜ao de

esparsidade da matriz do m´etodo de Galerkin.

6.1 Macroelementos e fun¸c˜oes bolha

Chamamos de macroelemento M

o suporte de uma fun¸c˜ao base nodal global

. Cada M

cobre os pontos nodais x

para todo j ∈ A

, onde A

= A

, conforme

deﬁnido na Se¸c˜ao 5.5.1.

Deﬁnimos o que chamamos de “bolhas” no macroelemento M

como sendo

as seguintes |A

| fun¸c˜oes:











se i = a

se i = a

(6.2)

para todo a ∈ |A

|. Por constru¸c˜ao, b

(x) = 0 se x ∈ ∂M

, uma vez que φ

(x) = 0

para x ∈ ∂M

A bolha central b

, que na verdade ´e a pr´opria φ

, tem como suporte todo

o macroelemento M

. Dentre as demais bolhas, algumas s˜ao restritas a um ´unico

elemento e as restantes s˜ao restritas a dois elementos, conforme ilustrado na Figura

6.1.

Cada fun¸c˜ao base nodal ψ

do espa¸co V

´e constru´ıda como uma combina¸c˜ao

linear das bolhas associadas ao macroelemento correspondente, ou seja, deﬁnimos,

para todo i ∈ N,



a∈A

. (6.3)

Assim, para que o espa¸co V

ﬁque completamente deﬁnido, ´e preciso deter-

minar para cada ponto nodal x

os |A

| parˆametros λ

. Nas Se¸c˜oes 6.2 e 6.3 ser˜ao

consideradas duas formas para determinar esses parˆametros.

supp(b

) cobrindo um ´unico elemento.

supp(b

) cobrindo dois elementos.

Figura 6.1: Em azul, supp(φ

). Em amarelo, supp(φ

). E em azul/amarelo

supp(φ

) = supp(b

6.2 M´etodo quase estabilizado

Depois de montado o sistema global de equa¸c˜oes do m´etodo de elementos

ﬁnitos formulado com o espa¸co V

constru´ıdo na se¸c˜ao anterior, as equa¸c˜oes asso-

ciadas aos pontos interiores x

∈ N podem ser expressas da mesma forma que na

equa¸c˜ao (5.8):



j∈A

= 0 para todo i ∈ ind (I), (6.4)

onde U

´e o grau de liberdade associado ao ponto x

= a(φ

, ψ

) =



a∈A

a(φ

, b

). (6.5)

Vale notar que S

nada mais ´e do que a matriz global associada ao m´etodo que

estamos propondo. Deﬁnindo, para cada ponto x

, a matriz

= a(φ

, b

), (6.6)

ent˜ao podemos expressar os coeﬁcientes S

como



a∈A

. (6.7)

Considerando ainda o vetor Λ

formado pelos parˆametros λ

e o vetor S

dos

coeﬁcientes S

, ent˜ao obtemos a representa¸c˜ao mais compacta

= D

. (6.8)

Essa dependˆencia linear entre os coeﬁcientes S

e os parˆametros λ

´e impor-

tante por raz˜oes pr´aticas. Se modiﬁc´assemos tamb´em o espa¸co U

(por exemplo

para for¸car a matriz global ser sim´etrica) essa dependˆencia linear se perderia.

Para obter um m´etodo com polui¸c˜ao m´ınima no mesmo sentido que o QS-

FEM, basta considerar uma malha uniforme de elementos quadrados, identiﬁcar

os λ

que devem ser iguais por quest˜ao de simetria e, como foi feito no QSFEM,

resolver a equa¸c˜ao de dispers˜ao com k

= k para as dire¸c˜oes ´otimas

3π

rela¸c˜ao aos parˆametros λ

restantes.

Mas, como j´a dito antes, para malhas mais gerais n˜ao h´a uma forma evidente

de escolher dire¸c˜oes ´otimas, sendo interessante ent˜ao aplicar o crit´erio usado no

m´etodo QOFD, que ´e ﬂex´ıvel o suﬁciente para determinar os λ

em situa¸c˜oes mais

gerais.

E o que fazemos na pr´oxima se¸c˜ao.

6.3 Fun¸c˜oes peso quase ´otimas

No m´etodo QOFD, o processo de minimiza¸c˜ao usado para determinar os

estˆenceis S

pˆode ser interpretado como a minimiza¸c˜ao do erro de truncamento do

operador discreto aplicado a ondas planas. Nos m´etodos de elementos ﬁnitos n˜ao

existe o conceito de erro de truncamento local, assim, n˜ao faz tanto sentido aplicar

aqui diretamente a estrat´egia usada no QOFD. Mas podemos ver aquele processo

de minimiza¸c˜ao de forma ligeiramente diferente.

100

Dada uma onda plana w(σ, x), levando em conta que o m´etodo de Petrov-

Galerkin ´e variacionalmente consistente, que qualquer onda plana ´e uma solu¸c˜ao

do problema homogˆeneo e que as φ

consideradas acima tˆem o suporte no interior

do dom´ınio, ent˜ao para toda fun¸c˜ao peso ψ

associada as n´os interiores, w satisfaz

a(w, ψ

) = 0, (6.9)

Mas se deﬁninirmos w

(σ, x) como a interpolante de w(σ, x), ent˜ao w

n˜ao satisfaz

o problema variacional discreto. Em vez disso temos um res´ıduo

R(σ, Λ

) = a(w

, ψ

) (6.10)

para cada n´o interior i. O que estamos propondo ´e minimizar o res´ıduo R(σ, Λ

) em

rela¸c˜ao a Λ

levando em conta todas as dire¸c˜oes poss´ıveis σ. Podemos interpretar

que escolhemos as interpolantes w

(σ, x) como alvo do m´etodo. Minimizando o

res´ıduo acima esperamos que as solu¸c˜oes aproximadas para essas ondas planas se

aproximem de suas interpolantes.

Assim, de forma semelhante ao que foi feito no Cap´ıtulo 5, deﬁnimos o fun-

cional de m´ınimos quadrados

J(Λ

) =



|R(σ, Λ

σ. (6.11)

Notando que

R(σ, Λ

) = a(w

, ψ

) (6.12)

= a







j∈ ind(N)

w(σ, x

)φ

, ψ





(6.13)



j∈A

w(σ, x

)a(φ

, ψ

) (6.14)



j∈A

w(σ, x

), (6.15)

101

ent˜ao

J(Λ

) =





m∈A



n∈A

w(σ, x

)w(σ, x

)dσ (6.16)



m∈A



n∈A



w(σ, x

)w(σ, x

)dσ (6.17)



m∈A



n∈A



ikσ·(x

−x

)

dσ (6.18)



m∈A



n∈A

(6.19)

= S

· W S

(6.20)

= (D

) · W D

(6.21)

= Λ



)

W D



, (6.22)

A matriz W ´e a obtida para o m´etodo QOFD na equa¸c˜ao (5.51), no caso bidimen-

sional ela ´e dada por

= 2πJ

(kh

) (6.23)

= 2π



1 −

(kh

)

(kh

)

−

(kh

)

2304

+ . . .



, (6.24)

com

= |x

− x

|. (6.25)

Para que a minimiza¸c˜ao de J(Λ

) n˜ao tenha como solu¸c˜ao Λ

= 0, devemos

impor alguma restri¸c˜ao em Λ

. Observando a deﬁni¸c˜ao de ψ

na equa¸c˜ao (6.3), ´e

razo´avel impor

= 1, (6.26)

para que tenhamos

) = 1. (6.27)

Sendo assim, os parˆametros Λ

s˜ao ent˜ao determinados minimizando J(Λ

)

102

em rela¸c˜ao aos |A

| − 1 parˆametros livres λ

, com os parˆametros centrais prede-

ﬁnidos λ

= 1. Chamamos a formula¸c˜ao obtida com esses parˆametros de QOPG

(Quasi Optimal Petrov-Galerkin) (Loula e Fernandes (submetido)).

Observa¸c˜ao 6.1 Um ponto importante ´e o fato de J(Λ

) ser local, no sentido de

que, para um n´o ﬁxado x

, J(Λ

) depende dos |A

| − 1 parˆametros λ

associados

apenas ao n´o x

. Com isso ´e poss´ıvel minimizar J(Λ

) em cada n´o separadamente.

Se em vez de V

modiﬁc´assemos U

(ou se modiﬁc´assemos os dois), essa condi¸c˜ao

de localidade n˜ao valeria, J(Λ

) dependeria tamb´em de alguns parˆametros associ-

ados aos pontos vizinhos, o que acoplaria todas as equa¸c˜oes que determinam os

parˆametros Λ

Observa¸c˜ao 6.2 Para melhor compreender esta estrat´egia de estabiliza¸c˜ao, vamos

considerar que para qualquer onda plana w(x, y) temos a seguinte estimativa para

a sua aproxima¸c˜ao w

, em rela¸c˜ao a w

(a interpolande de w):

w

− w

 ≤ C, (6.28)

como resultado de consistˆencia (no sentido de diferen¸cas ﬁnitas) e estabilidade no

sentido variacional (de elementos ﬁnitos). Pela desigualdade do triangulo temos

w − w

 ≤ C(w − w

 + w

− w

) ≤ C(w − w

 + ) (6.29)

Logo, reduzir a diferen¸ca entre w

− w

 corresponde exatamente a reduzir o

erro de polui¸c˜ao. Essa estrat´egia ´e comum a quase todos os m´etodos estabili-

zados. De modo que ela pode ser aplicada a muitos outros problemas que de-

mandam estabiliza¸c˜oes, tais como problemas de convec¸c˜ao-difus˜ao ou convec¸c˜ao

rea¸c˜ao predominantemente convecctivos ou reativos. Dado que as solu¸c˜oes da

equa¸c˜ao de Helmholtz homogˆenea podem ser aproximadas com muita precis˜ao por

combina¸c˜oes lineares de ondas planas (Monk e Wang (1999)), esta estrat´egia n˜ao

est´a restrita a ondas planas. Resultados num´ericos que ser˜ao apresentados neste

cap´ıtulo conﬁrmam esta observa¸c˜ao.

103

6.4 Condi¸c˜oes de contorno

Com condi¸c˜oes de contorno de Dirichlet, o QOPG gera equa¸c˜oes associadas

apenas a pontos interiores. Equa¸c˜oes associadas a n´os de fronteira ocorrem quando

tratamos de condi¸c˜oes de Neumann ou Robin por exemplo. Isso deve ser levado em

conta no m´etodo em quest˜ao pois o n´umero de n´os adjacentes e portanto o n´umero

de parˆametros livres associados aos n´os de fronteira ´e menor do que os associados

aos n´os interiores; o que pode afetar a efetividade do processo de minimiza¸c˜ao

descrito acima. Na pr´atica pode-se veriﬁcar com experimentos num´ericos que o

m´etodo proposto ´e seriamente afetado quando n˜ao se d´a algum tratamento especial

aos n´os de fronteira sujeitos a condi¸c˜oes de contorno de Neumann ou Robin.

Para tratar adequadamente estas condi¸c˜oes de contorno n˜ao essenciais, duas

formas distintas foram consideradas:

(1) Uma camada de elementos ´e adicionada no exterior do dom´ınio com o

intuito de completar os estˆenceis de fronteira. Em seguida os parˆametros

s˜ao determinados levando em contas esses elementos extras. Finalmente

o m´etodo ´e completado apenas com os elementos originais.

(2) Nos n´os de fronteira os parˆametros Λ

s˜ao determinados usando o processo

de minimiza¸c˜ao semelhante ao normal mas o res´ıduo (6.10) ´e substitu´ıdo

por outro calculado pela forma sesquilinear modiﬁcada

(u, v) =



Ω

∇u∇v − k

uv dx −



∂Ω

∇u · nv ds. (6.30)

Ap´os a determina¸c˜ao dos parˆametros a forma a normal ´e usada para mon-

tar o sistema de equa¸c˜oes global.

No caso de uma malha uniforme os m´etodos (1) e (2) produzem resultados

com a mesma precis˜ao. Nos experimentos num´ericos deste cap´ıtulo foi usado o

esquema (1).

104

6.5 Implementa¸c˜ao computacional

Diferentemente do que ´e feito nos m´etodos de elementos ﬁnitos usuais, no

QOPG a matriz global deve ser montada equa¸c˜ao por equa¸c˜ao. Isso se deve ao

fato de que os parˆametros de estabiliza¸c˜ao λ

s˜ao calculados separadamente em

cada n´o, resolvendo, para cada x

os sistema local



)

W D



= 0. (6.31)

Sendo mais preciso, n˜ao ´e exatamente o sistema acima que ´e resolvido. O parˆametro

´e prescrito e portanto uma das equa¸c˜oes do sistema acima n˜ao precisa ser satis-

feita no processo de minimiza¸c˜ao descrito na Se¸c˜ao 6.3.

Uma vez calculados esses parˆametros para um n´o x

, os coeﬁcientes da

equa¸c˜ao global associada a x

s˜ao obtidos calculando

= D

. (6.32)

Do mesmo modo que no m´etodo QOFD, deve-se ter cuidado ao resolver os

sistemas locais



)

W D



= 0, dado o mal condicionamento das matrizes W .

6.6 Resultados num´ericos

Nesta se¸c˜ao apresentamos resultados de experimentos num´ericos ilustrando

o desempenho da formula¸c˜ao proposta considerando malhas uniformes e n˜ao uni-

formes.

6.6.1 Malhas uniformes

Consideramos primeiramente o m´etodo QOPG em malhas uniformes, com

condi¸c˜oes de contorno de Dirichlet e em seguida com condi¸c˜oes de Robin.

Condi¸c˜oes de contorno de Dirichlet.

Na Figura 6.2 apresentamos o erro relativo na norma do L

(Ω) para os

105

m´etodos de Galerkin, GLS e QOPG, comparados com a interpolante de uma onda

plana em diferentes dire¸c˜oes θ ∈ [0, π/2], em uma malha com 50 × 50 elementos,

com k = 50 (kh = 1).

Efeitos de polui¸c˜ao deterioram completamente a aproxima¸c˜ao de Galerkin: o

erro relativo em L

(Ω) chega a ser duas ordens de magnitude maior do que o erro

da interpolante. O erro relativo da aproxima¸c˜ao pelo GLS ´e uma ordem menor

comparado `a aproxima¸c˜ao de Galerkin, mas maior do que o erro da aproxima¸c˜ao

do QOPG..

3π

-1

-0.5

0.5

log (u − u

/u)

Interpolant

Galerkin

GLS

QOPG

Figura 6.2: Erros relativos na norma L

para ondas planas em diferentes dire¸c˜oes,

com k = 50, condi¸c˜oes de contorno de Dirichlet e uma malha com 50×50 elementos

quadrados.

Na Figura 6.3 vemos resultados de convergˆencia para uma onda plana na

dire¸c˜ao θ = π/4 com k = 50, para o erro relativo na norma do L

(Ω) e na semi-

norma do H

(Ω). Novamente observamos claramente polui¸c˜ao na convergˆencia dos

m´etodos de Galerkin e GLS. A convergˆencia para o QOPG ´e praticamente igual `a

da interpolante.

Resultados similares s˜ao apresentados para uma onda plana com k = 100 na

Figura 6.4, onde os erros relativos em L

(Ω) s˜ao mostrados para ondas planas em

diferentes dire¸c˜oes θ ∈ [0, π/2], e na Figura 6.5, onde s˜ao mostrados estudos de

convergˆencia para θ = π/4. Como esperado, para valores maiores de k, os efeitos

de polui¸c˜ao s˜ao agravados. Nesse caso, as aproxima¸c˜oes de Galerkin e do GLS n˜ao

exibem convergˆencia. Para evitar esse tipo de problema consideraremos condi¸c˜oes

106

1.6 1.8 2 2.2 2.4

−log(h)

-2.5

-2

-1.5

-1

-0.5

0.5

log (u − u

/u)

Interpolant

Galerkin

GLS

QOPG

1.6 1.8 2 2.2 2.4

−log(h)

-1.5

-1

-0.5

0.5

log (|u − u

|/|u|)

Interpolant

Galerkin

GLS

QOPG

Figura 6.3: Estudo de convergˆencia para condi¸c˜oes de contorno e Dirichlet para

uma onda plana com k = 50 and θ = π/4.

de Robin daqui por diante.

3π

-1

-0.5

0.5

log (u − u

/u)

Interpolant

Galerkin

GLS

QOPG

Figura 6.4: Erros relativos na norma L

para ondas planas em diferentes dire¸c˜oes,

com k = 100, condi¸c˜oes de contorno de Dirichlet e uma malha com 100 × 100

elementos quadrados.

Condi¸c˜oes de contorno de Robin.

Apresentamos agora alguns resultados com k = 100 para condi¸c˜oes de con-

torno de Robin tais que a solu¸c˜ao exata seja uma onda plana na dire¸c˜ao θ.

Na Figura 6.6 apresentamos a dependˆencia do erro em L

(Ω) para θ ∈

[0, π/2]. Como n˜ao h´a problemas de ressonˆancia num´erica, um comportamento

107

2 2.2 2.4 2.6

−log(h)

-2

-1

log (u − u

/u)

Interpolant

Galerkin

GLS

QOPG

2 2.2 2.4 2.6

−log(h)

-1.5

-1

-0.5

0.5

1.5

log (|u − u

|/|u|)

Interpolant

Galerkin

GLS

QOPG

Figura 6.5: Estudo de convergˆencia para condi¸c˜oes de contorno e Dirichlet para

uma onda plana com k = 100 and θ = π/4.

regular pode ser observado para as aproxima¸c˜oes de Galerkin e do GLS, mas o

efeito de polui¸c˜ao ainda est´a presente.

A Figura 6.7 apresenta resultados de convergˆencia para uma onda plana com

θ = π/4. Ainda podemos observar bastante polui¸c˜ao nas curvas de convergˆencia

dos m´etodos de Galerkin e GLS, mas menos severa do que com condi¸c˜oes de con-

torno de Dirichlet. No entanto, o erro relativo da aproxima¸c˜ao de Galerkin na

malha menos reﬁnada ´e uma ordem de magnitude maior do que o erro da inter-

polante, tanto na norma L

(Ω) quanto na seminorma do H

(Ω). A aproxima¸c˜ao

do GLS ´e mais precisa, em particular na seminorma do H

(Ω). A aproxima¸c˜ao do

QOPG exibe a mesma precis˜ao da interpolante em H

(Ω) e at´e mais precis˜ao que

a interpolante na norma do L

(Ω).

6.6.2 Malhas n˜ao uniformes. Perturba¸c˜oes aleat´orias

Nesta se¸c˜ao s˜ao realizados testes com elementos distorcidos. A distor¸c˜ao ´e ob-

tida substituindo as coordenadas (x

, y

) dos n´os interiores de uma malha formada

originalmente por elementos quadrados por coordenadas (x



, y



) aleatoriamente

108

3π

-1

-0.5

log (u − u

/u)

Interpolant

Galerkin

GLS

QOPG

Figura 6.6: Erros relativos na norma L

para ondas planas em diferentes dire¸c˜oes,

com k = 100, condi¸c˜oes de contorno de Robin e uma malha com 100×100 elementos

quadrados.

2 2.2 2.4 2.6

−log(h)

-2.5

-2

-1.5

-1

-0.5

0.5

log (u − u

/u)

Interpolant

Galerkin

GLS

QOPG

2 2.2 2.4 2.6

−log(h)

-1.5

-1

-0.5

0.5

log (|u − u

|/|u|)

Interpolant

Galerkin

GLS

QOPG

Figura 6.7: Estudo de convergˆencia para condi¸c˜oes de contorno de Robin para uma

onda plana com k = 100 and θ = π/4.

perturbadas da seguinte forma:



= x

+ r

h, (6.33)



= y

+ s

h, (6.34)

onde r

e s

s˜ao sequˆencias de n´umeros aleat´orios distribu´ıdos uniformemente no

intervalo [−0.24, 0.24].

109

Figura 6.8: Malha distorcida aleatoriamente

Resultados de convergˆencia para as aproxima¸c˜oes de Galerkin, GLS e QOPG,

comparadas com a interpolante, s˜ao apresentados na Figura 6.9.

E interessante

notar a robustez das aproxima¸c˜oes por elementos ﬁnitos. As aproxima¸c˜oes de Ga-

lerkin e GLS s˜ao melhores do que as mesmas em malhas uniformes. A aproxima¸c˜ao

do QOPG tem praticamente a mesma precis˜ao obtida com malhas uniformes.

Devemos notar que o parˆametro de estabiliza¸c˜ao do GLS ´e obtido por meio

de an´alise espectral numa malha uniforme. Por constru¸c˜ao, a aproxima¸c˜ao de GLS

de uma onda plana ´e indˆentica `a interpolante quando θ = π/8 e 3π/8, como pode-

mos ver nas Figuras 6.2 e 6.4. No entanto, o mesmo n˜ao ocorre necessariamente

110

2 2.2 2.4 2.6

−log(h)

-2.5

-2

-1.5

-1

-0.5

0.5

log (u − u

/u)

Interpolant

Galerkin

GLS

QOPG

2 2.2 2.4 2.6

−log(h)

-1.5

-1

-0.5

0.5

log (|u − u

|/|u|)

Interpolant

Galerkin

GLS

QOPG

Figura 6.9: Estudo de convergˆencia para condi¸c˜oes de Robin, malha n˜ao uniforme,

k = 100 e θ = π/4.

com malhas n˜ao uniformes. Para checar isso, apresentamos na Figura 6.10 o erro

relativo em L

para as aproxima¸c˜oes de Galerkin, GLS e QOPG para ondas planas

com v´arias dire¸c˜oes em θ ∈ [0, π/2], considerenado malhas com 100×100 elementos

distorcidos. O erro relativo m´ınimo n˜ao ocorre em θ = π/8 e θ = 3π/8, mas em

valores de θ pr´oximos de 3π/16 e 5π/16. Notamos tamb´em que o erro m´aximo

n˜ao ocorre em θ = π/4, mas em θ = 0 and θ = π/2. Assim, com essas malhas

n˜ao uniformes, o efeito de polui¸c˜ao deve ser mais intenso para ondas planas em

θ = 0 ou θ = π/2. Para ilustrar isso, apresentamos resultados de convergˆencia

para θ = 0 na Figura 6.11, para k = 100. A polui¸c˜ao nas solu¸c˜oes aproximandas

de Galerkin e GLS se torna mais intensa nesse caso.

6.6.3 Malhas n˜ao uniformes. Transforma¸c˜ao biquadr´atica

Conseideramos agora malhas n˜ao uniformes para o dom´ınio (0, 1) × (0, 1)

geradas transformando os n´os de uma malha uniformes por uma transforma¸c˜ao

biquadr´atica. A Figura 6.12 mostra uma malha com 20 × 20 elementos gerada

dessa forma. Podemos notar a ocorrˆencia de elementos muito distorcidos.

Novamente estudos de convergˆencia s˜ao feitos para a equa¸c˜ao de Helmholtz

111

3π

-1

-0.5

log (u − u

/u)

Interpolant

Galerkin

GLS

QOPG

Figura 6.10: Erros relativos em L

para ondas planas em diferentes dire¸c˜oes θ ∈

[0, π/2] com k = 100, Condi¸c˜oes de contorno de Robin em uma malha n˜ao uniforme

com 100 × 100 nelementos.

2 2.2 2.4 2.6

−log(h)

-2.5

-2

-1.5

-1

-0.5

0.5

log (u − u

/u)

Interpolant

Galerkin

GLS

QOPG

2 2.2 2.4 2.6

−log(h)

-1.5

-1

-0.5

0.5

log (|u − u

|/|u|)

Interpolant

Galerkin

GLS

QOPG

Figura 6.11: Estudo de convergˆencia com condi¸c˜oes de contorno de Robin. Apro-

xima¸c˜oes de Galerkin, GLS and QOPG comparadas com a interpolante, malhas

n˜ao uniformes, k = 100 e θ = 0

com k = 100, considerando condi¸c˜oes de contorno de Robin apropriadas para que

a solu¸c˜ao seja uma onda plana na dire¸c˜ao θ = π/4.

Resultados de convergˆencia s˜ao apresentados nas ﬁguras Figuras 6.13 e 6.14.

Novamente podemos observar a robusteza do m´etodo QOPG, que produz solu¸c˜oes

com a mesma precis˜ao da interpolante na seminorma do H

e precis˜ao maior que

a interpolante na norma L

112

Figura 6.12: Malha n˜ao uniforme gerada aplicando uma uma transforma¸c˜ao biqua-

dr´atica a uma malha uniforme. Uma camada extra de elementos foi adicionada no

contorno para tratar condi¸c˜oes de Robin.

3π

-1

-0.5

log (u − u

/u)

Interpolant

Galerkin

GLS

QOPG

Figura 6.13: Erros relativos na norma L

para ondas planas em diferentes dire¸c˜oes

θ ∈ [0, π/2] com k = 100, condi¸c˜oes de contorno de Robin e uma malha distorcida

por um mapeamento biquadr´atico com 100 × 100 elementos.

6.6.4 Problema de Helmholtz n˜ao homogˆeneo

O QOPG ´e formulado com foco em solu¸c˜oes que s˜ao ondas planas. Al´em

disso, a matriz global resultante do QOPG n˜ao ´e sim´etrica em geral. Assim,

n˜ao podemos esperar taxas de convergˆencia ´otimas em L

para aproxima¸c˜oes com

termos de fonte.

O objetivo desta se¸c˜ao ´e checar o desempenho da formula¸c˜ao QOPG para

113

2 2.2 2.4 2.6

−log(h)

-2.5

-2

-1.5

-1

-0.5

0.5

log (u − u

/u)

Interpolant

Galerkin

GLS

QOPG

2 2.2 2.4 2.6

−log(h)

-1.5

-1

-0.5

0.5

log (|u − u

|/|u|)

Interpolant

Galerkin

GLS

QOPG

Figura 6.14: Estudo de convergˆencia para condi¸c˜oes de contorno de Robin. Erros

relativos na norma L

(direita) e seminorma H

(esquerda). Onda plana com

k = 100 e θ = π/4.

problema com termos de fonte suaves. Resultados s˜ao apresentados para k = 10 e

k = 100 considerando condi¸c˜oes de contorno de Robin e um termo de fonte f(x, y)

tal que a solu¸c˜ao exata em um quadrado unit´ario seja

u(x, y) = sin(πx) sin(πy) . (6.35)

Para estudos de convergˆencia usando sequˆencias de malhas uniformes, sempre

obtivemos taxas de convergˆencia ´otimas para aproxima¸c˜oes com o QOPG na norma

do L

e seminorma do H

. Para malhas muito distorcidas foram obtidas taxas n˜ao

´otimas na norma do L

. Para malhas moderadamente distorcidas foram tamb´em

obtidas taxas ´otimas em L

, como podemos observar na Figura 6.15 para k = 10 e

na Figura 6.16 para k = 100. Nesses estudos consideramos reﬁnamentos de malhas

n˜ao uniformes, geradas pelo mesmo mapeamento biquadr´atico como ilustrado na

Figura 6.12.

114

1 1.5 2

−log(h)

-5

-4.5

-4

-3.5

-3

-2.5

-2

-1.5

-1

log (u − u

/u)

Interpolant

Galerkin

GLS

QOPG

1 1.5 2

−log(h)

-2.2

-2

-1.8

-1.6

-1.4

-1.2

-1

-0.8

log (|u − u

|/|u|)

Interpolant

Galerkin

GLS

QOPG

Figura 6.15: Estudo de convergˆencia para condi¸c˜oes de contorno de Robin. Erros

relativos na norma L

(direita) e seminorma H

(esquerda). k = 10 com termo de

fonte.

1.8 2 2.2 2.4 2.6 2.8

−log(h)

-5.5

-5

-4.5

-4

-3.5

log (u − u

/u)

Interpolant

Galerkin

GLS

QOPG

1.8 2 2.2 2.4 2.6 2.8

−log(h)

-2.8

-2.6

-2.4

-2.2

-2

-1.8

log (|u − u

|/|u|)

Interpolant

Galerkin

GLS

QOPG

Figura 6.16: Estudo de convergˆencia para condi¸c˜oes de contorno de Robin. Erros

relativos na norma L

(direita) e seminorma H

(esquerda). k = 100 com termo

de fonte.

6.6.5 Ondas n˜ao direcionais

Nesta se¸c˜ao consideramos a equa¸c˜ao de Helmholtz num dom´ılio anelar com

raio interior R

= 0.5 e raio exterior R

= 1, k = 100 e f(x, y) = 0, com condi¸c˜oes

de contorno de Dirichlet em r = R

e de Robin em R = R

, tais que a solu¸c˜ao em

115

coordenadas polares seja dada pela solu¸c˜ao descrita na Se¸c˜ao 2.3.3:

(r, θ) = cos(nθ)H

(1)

(kr), (6.36)

As malhas usadas neste estudo s˜ao constru´ıdas da mesma forma descrita na

Se¸c˜ao 4.4.4 (Figura 4.13).

Estudos de convergˆencia para o erro relativo na norma L

e na seminorma

s˜ao apresentados para n = 0 na Figura 6.17, n = 3 na Figura 6.18 e n = 7 na

Figura 6.19. As aproxima¸c˜oes de Galerkin e GLS s˜ao novamente muito afetadas

pela polui¸c˜ao, enquanto a aproxima¸c˜ao pelo QOPG s˜ao t˜ao precisas quanto a

interpolante.

1.8 2 2.2 2.4 2.6

−log(h)

-2.5

-2

-1.5

-1

-0.5

0.5

log (u − u

/u)

Interpolant

Galerkin

GLS

QOPG

1.8 2 2.2 2.4 2.6

−log(h)

-1.5

-1

-0.5

0.5

log (|u − u

|/|u|)

Interpolant

Galerkin

GLS

QOPG

Figura 6.17: Convergˆencia para uma onda n˜ao direcional (eq. (6.36) com n=0).

116

1.8 2 2.2 2.4 2.6

−log(h)

-2.5

-2

-1.5

-1

-0.5

0.5

log (u − u

/u)

Interpolant

Galerkin

GLS

QOPG

1.8 2 2.2 2.4 2.6

−log(h)

-1.5

-1

-0.5

0.5

log (|u − u

|/|u|)

Interpolant

Galerkin

GLS

QOPG

Figura 6.18: Convergˆencia para uma onda n˜ao direcional (eq. (6.36) com n=3).

1.8 2 2.2 2.4 2.6

−log(h)

-2.5

-2

-1.5

-1

-0.5

0.5

log (u − u

/u)

Interpolant

Galerkin

GLS

QOPG

1.8 2 2.2 2.4 2.6

−log(h)

-1.5

-1

-0.5

0.5

log (|u − u

|/|u|)

Interpolant

Galerkin

GLS

QOPG

Figura 6.19: Convergˆencia para uma onda n˜ao direcional (eq. (6.36) com n=7).

117

Cap´ıtulo 7

Conclus˜oes e futuros desenvolvimentos

Iniciamos este trabalho com o objetivo de desenvolver m´etodos de elementos

ﬁnitos e de diferen¸cas ﬁnitas para o problema de Helmholtz com n´umeros de onda

elevados, com boas propriedades de estabilidade, precis˜ao e rubustez a distor¸c˜oes

da malha. Apesar de n˜ao termos apresentado uma an´alise matem´atica e num´erica

completa dos m´etodos propostos, o grande n´umero de estudos de convergˆencia

realizados, comparativamente com outros m´etodos, mostra que o objetivo inicial

foi alcan¸cado em grande parte.

Inicialmente ﬁzemos uma revis˜ao da literatura, particularmente sobre m´eto-

dos de elementos ﬁnitos baseados na formula¸c˜ao cl´assica de Galerkin e em t´ecnicas

de estabiliza¸c˜ao residuais do tipo Galerkin M´ınimos Quadrados. Neste contexto

propomos dois m´etodos estabilizados: o primeiro combina res´ıduos de Galerkin,

m´ınimos quadrados e de volumes ﬁnitos e o segundo utiliza res´ıduos de Galerkin

e m´ınimos quadrados ponderados. Ambos s˜ao capazes de minimizar os efeitos

de polui¸c˜ao do erro, no sentido deﬁnido por Babuˇska, quando restritos a malhas

uniformes de elementos quadrados. No entanto, estes m´etodos n˜ao apresentam a

robustez a distor¸c˜oes da malha que buscamos.

Em continua¸c˜ao desenvolvemos novos m´etodos de diferen¸cas ﬁnitas gene-

ralizadas aplic´aveis a malhas n˜ao uniformes e n˜ao estruturadas, bem como no-

vos m´etodos de elementos ﬁnitos apoiados em formula¸c˜oes conformes de Petrov-

Galerkin que se mostraram bastante precisos e robustos a distor¸c˜oes da malha.

118

As an´alises das ordens de aproxima¸c˜ao e de dispers˜ao em malhas uniformes e uma

s´erie bastante ampla de experimentos num´ericos realizados suportam as conclus˜oes

apresentadas a seguir.

Sobre diferen¸cas ﬁnitas

(1) O m´etodo de diferen¸cas ﬁnitas QOFD (Quasi Optimal ﬁnite Diﬀerence) foi

desenvolvido de forma diferente da usual, no sentido de que foi formulado

de modo a ser aplic´avel a qualquer distribui¸c˜ao dos pontos onde a solu¸c˜ao

do problema ´e aproximada, para qualquer n´umero n ≥ 1 de dimens˜oes.

(2) O QOFD foi constru´ıdo de forma que o erro de truncamento do operador

de Helmholtz aproximado fosse minimizado para ondas planas, levando

em conta igualmente todas as dire¸c˜oes poss´ıveis. O resultado ´e equivalente

em termos de ordem de aproxima¸c˜ao ao QSFEM em malhas uniformes de

elementos quadrados, mas, diferentemente do QSFEM, o QOFD ´e aplic´avel

a malhas n˜ao uniformes e n˜ao estruturadas.

(3) Resultados num´ericos sobre uma s´erie de problemas em 2D mostram taxas

de convergˆencia ´otimas na seminorma de H

e na norma de L

, equivalentes

`a interpolante em malhas n˜ao uniformes mas estruturadas, com estˆenceis

de 9 pontos.

(4) Estudos de convergˆencia com malhas n˜ao estruturadas mostram igual-

mente taxas de convergˆencia ´otimas nestas normas.

(5) Tamb´em foram apresentadas generaliza¸c˜oes de aproxima¸c˜oes com estˆenceis

envolvendo um n´umero qualquer de pontos. Nesse caso, especial aten¸c˜ao

deve ser dada ao problema num´erico, dado o mal condicionamento das

matrizes dos problemas locais. Quando esses problemas num´ericos s˜ao

contornados, solu¸c˜oes com grande precis˜ao podem ser obtidas aumentando

o n´umero de pontos dos estˆenceis.

119

Sobre elementos ﬁnitos

(1) Foram formulados m´etodos de elementos ﬁnitos de Petrov-Galerkin para

espa¸cos de aproxima¸c˜ao lagrangianos C

lineares ou bilineares em dom´ınios

bidimensionais, que resultam em matrizes globais com a mesma esparsi-

dade da matriz global do m´etodo de Galerkin. Em cada m´etodo, o espa¸co

das fun¸c˜oes peso foi modiﬁcado usando bolhas C

polinomiais por partes

em macroelementos, visando obter estabilidade, precis˜ao e robustez.

(2) O RPPG (Reduced Pollution Petrov-Galerkin), o mais simples dos m´etodos

de elementos ﬁnitos desenvolvido neste trabalho, mostrou-se bastante eﬁ-

ciente em reduzir o efeito de polui¸c˜ao em malhas com poucas distor¸c˜oes.

(3) O RPPG n˜ao tem parˆametros de estabiliza¸c˜ao dependentes do n´umero de

onda ou dos tamanhos dos elementos como ´e comum nos m´etodos estabi-

lizados para o problema de Helmholtz.

(4) Visando uma maior robustez a distor¸c˜oes da malha, o QSPG Quasi Stabi-

lized Petrov-Galerkin foi desenvolvido de forma semelhante ao RPPG mas

usando parˆametros de estabiliza¸c˜ao dependentes de k e da geometria dos

elementos. Com isso conseguimos um m´etodo equivalente ao QSFEM nas

situa¸c˜oes em que esses m´etodos podem ser comparados diretamente.

(5) O QSPG ´e aplic´avel naturalmente a elementos distorcidos de forma razo-

avelmente robusta, especialmente para elementos muito alongados, uma

situa¸c˜ao que geralmente causa problemas para muitos m´etodos estabiliza-

dos propostos na literatura.

(6) Para algumas distor¸c˜oes mais extremas, o QOFD ´e bem mais robusto do

que o RPPG e o QSPG. Com isso em mente, o QOPG (Quasi Optimal

Petrov-Galerkin) foi desenvolvido, tendo como objetivo aliar a robustez

resultante do crit´erio usado para formular o QOFD `as vantagens de um

120

m´etodo variacionalmente consistente. Os estudos num´ericos indicam que

esse objetivo foi alcan¸cado em grande parte.

(7) Foram feitos estudos de convergˆencia para ondas planas, ondas direcio-

nais e ondas evanescentes, bem como para solu¸c˜oes n˜ao homogˆeneas do

problema de Helmholtz com condi¸c˜oes de contorno de Diriclet, Robin e

mistas. Em todos estes estudos foram observadas taxas de convergˆencia

´otimas e precis˜oes equivalentes `a interpolante na seminorma de H

. Para

condi¸c˜oes de contorno de Robin foram observadas precis˜oes melhores do

que a da interpolante na norma do L

(8) Dentre todas as formula¸c˜oes aqui estudadas, no contexto dos m´etodos

de elementos ﬁnitos lagrangianos lineares ou bilineares, a formula¸c˜ao de

Petrov-Galerkin QOPG ´e a que apresenta melhores propriedades de esta-

bilidade, precis˜ao, generalidade e robustez a distor¸c˜oes de malhas.

Futuros desenvolvimentos

Conforme resumimos acima, conclu´ımos este trabalho com alguns resultados

que julgamos signiﬁcativos. Infelizmente, ou felizmente, algumas quest˜oes conti-

nuam em aberto, como a pr´opria an´alise num´erica dos m´etodos propostos. Surgem

outras que podem ser relevantes, como uma maior conex˜ao entre os m´etodos de

elementos ﬁnitos e de diferen¸cas ﬁnitas generalizados, assim como uma nova es-

trat´egia para determina¸c˜ao dos parˆametros de estabiliza¸c˜ao de m´etodos de elemen-

tos ﬁnitos ou de diferen¸cas ﬁnitas atrav´es da minimiza¸c˜ao do erro de truncamento

da interpolante de solu¸c˜oes fundamentais do problema de Helmholtz homogˆeneo.

Esta estrat´egia n˜ao est´a limitada ao Problema de Helmholtz, e pode vir a ser bem

explorada em outros problemas que tamb´em demandam estabiliza¸c˜ao, como, por

exemplo, problemas de convec¸c˜ao difus˜ao ou de rea¸c˜ao difus˜ao predominantemente

convectivos ou reativos.

Por ﬁm, listamos algumas destas quest˜oes que ﬁcaram em aberto, bem como

121

algumas poss´ıveis extens˜oes dos resultados obtidos:

• A an´alise num´erica de todos os m´etodos propostos. Quest˜ao ainda em

aberto e dif´ıcil, visto que mesmo o m´etodo de Galerkin carece de uma

an´alise completa para o problema de Helmholtz em mais de uma dimens˜ao.

• Extens˜oes do RPPG e do QSPG para problemas tridimensionais.

• Extens˜ao do QOPG para problemas n˜ao homogˆeneo com termos de fonte

singulares.

• Tratamento de condi¸c˜oes de contorno n˜ao essenciais no QOFD.

• Extens˜oes dos m´etodos propostos para o problema de Helmholtz vetorial

(eletromagnetismo).

• Adaptatividade.

E interessante observar que o funcional de m´ınimos qua-

drados usado na determina¸c˜ao dos parˆametros de estabiliza¸c˜ao do QOPG,

bem como dos coeﬁcientes do estˆencil do QOFD, ´e um indicador a priori

de qualidade da aproxima¸c˜ao da solu¸c˜ao do problema homogˆeneo.

122

Referˆencias Bibliogr´aﬁcas

Gustavo Benitez Alvarez, Abimael Fernando Dourado Loula, Eduardo Gomes Du-

tra do Carmo, e Fernando Alves Rochinha. A discontinuous ﬁnite element for-

mulation for Helmholtz equation. Computer Methods in Applied Mechanics and

Engineering, 195(33-36):4018 – 4035, 2006.

Ivo M. Babuˇska e Stefan A. Sauter. Is the pollution eﬀect of the FEM avoidable

for the Helmholtz equation considering high wave numbers? SIAM Journal on

Numerical Analysis, 34(6):2392–2423, 1997.

Ivo M. Babuˇska, F. Ihlenburg, Stefan A. Sauter, e E. T. Paik. A generalized

ﬁnite element method for solving the Helmholtz equation in two dimensions

with minimal pollution. Comp. Meth. Appl. Mech. Eng., 128:325–359, 1995.

Eduardo Gomes Dutra do Carmo, Gustavo Benitez Alvarez, Abimael Fer-

nando Dourado Loula, e Fernando Alves Rochinha. A nearly optimal Galer-

kin projected residual ﬁnite element method for Helmholtz problem. Computer

Methods in Applied Mechanics and Engineering, 197(13-16):1362 – 1375, 2008.

ISSN 0045-7825.

Jim Douglas, Juan E. Santos, Dongwoo Sheen, e Lynn Schreyer Bennethum. Fre-

quency domain treatment of one-dimensional scalar waves. Mathematical Models

and Methods in Applied Sciences, 3(2):171– 194, 1993.

Charbel Farhat, Isaac Harari, e Ulrich Hetmaniuk. A discontinuous Galerkin

method with lagrange multipliers for the solution of Helmholtz problems in the

123

mid-frequency regime. Computer Methods in Applied Mechanics and Enginee-

ring, 192(11-12):1389 – 1419, 2003.

Daniel Thomes Fernandes e Abimael Loula. M´etodos de elementos ﬁnitos estabi-

lizados para o problema de Helmholtz. XXVII CILAMCE, 2006.

Daniel Thomes Fernandes e Abimael Loula. A ﬁnite diﬀerence method for

Helmholtz equation for unstructured meshes. 8th. World Congress on Com-

putational Mechanics, 2008.

Leopoldo P. Franca, Charbel Farhat, Antonini P. Macedo, e Michel Lesoinne.

Residual-free bubbles for the Helmholtz equation. International Journal for

Numerical Methods in Engineering, 40(21):4003–4009, 1997.

Murthy N. Guddati e Bin Yue. Modiﬁed integration rules for reducing dispersion

error in ﬁnite element methods. Computer Methods in Applied Mechanics and

Engineering, 193(3-5):275–287, 2004.

Isaac Harari e Kirill Gosteev. Bubble-based stabilization for the Helmholtz equa-

tion. International Journal for Numerical Methods in Engineering, 70(10):1241–

1260, 2007.

Isaac Harari e Thomas J. R. Hughes. Finite element methods for the Helmholtz

equation in an exterior domain: Model problems. Computer Methods in Applied

Mechanics and Engineering, 87:59–96, 1991.

F. Ihlenburg e Ivo M. Babuˇska. Finite element solution of the Helmholtz equa-

tion with high wave number part i: The h-version of the fem. Computers &

Mathematics with Applications, (30, issue 9):9–37, 1995a.

F. Ihlenburg e Ivo M. Babuˇska. Finite element solution of the Helmholtz equation

with high wave number part ii: The hp-version of the fem. SIAM Journal on

Numerical Analysis, (34, issue 1):315–358, 1997.

124

Frank Ihlenburg. Finite Element Analysis of Acoustic Scattering. Springer, New

York, 1 edition, 1998.

Frank Ihlenburg e Ivo Babuˇska. Dispersion analysis and error estimation of Ga-

lerkin ﬁnite element methods for the Helmholtz equation. International Journal

for Numerical Methods in Engineering, 38(22):3745–3774, 1995b.

Abimael Loula e Daniel Thomes Fernandes. A quasi optimal Petrov-Galerkin

method for Helmholtz problem. International Journal for Numerical Methods

in Engineering, submetido.

Abimael Loula e Daniel Thomes Fernandes. Quasi optimal Petrov-Galerkin

methods for Helmholtz problem. 8th. World Congress on Computational Me-

chanics, 2008.

Abimael F.D. Loula, Gustavo B. Alvarez, Eduardo G.D. do Carmo, e Fernando A.

Rochinha. A discontinuous ﬁnite element method at element level for Helmholtz

equation. Computer Methods in Applied Mechanics and Engineering, 196(4-6):

867 – 878, 2007.

Jens Markus Melenk. On generalized ﬁnite element methods. PhD thesis, Univer-

sity of Maryland, 1995.

P. Monk e Da-Qing Wang. A least-squares method for the Helmholtz equation.

Computer Methods in Applied Mechanics and Engineering, 175(1-2):121 – 136,

1999.

Assad A. Oberai e Peter M. Pinsky. A residual-based ﬁnite element method for

the Helmholtz equation. International Journal for Numerical Methods in Engi-

neering, 49(3):399–419, 2000.

Jonathan Richard Shewchuk. Triangle: Engineering a 2D Quality Mesh Generator

and Delaunay Triangulator. In Ming C. Lin e Dinesh Manocha, editors, Ap-

plied Computational Geometry: Towards Geometric Engineering, volume 1148

125

of Lecture Notes in Computer Science, pages 203–222. Springer-Verlag, May

1996. From the First ACM Workshop on Applied Computational Geometry.

I. Singer e E. Turkel. High-order ﬁnite diﬀerence methods for the Helmholtz equa-

tion. Computer Methods in Applied Mechanics and Engineering, 163(1-4):343 –

358, 1998.

Godehard Sutmann. Compact ﬁnite diﬀerence schemes of sixth order for the

Helmholtz equation. Journal of Computational and Applied Mathematics, 203

(1):15 – 31, 2007.

Lonny L. Thompson e Peter M. Pinsky. A Galerkin least-squares ﬁnite element

method for the two-dimensional Helmholtz equation. International Journal for

Numerical Methods in Engineering, 38(3):371–397, 1995.

Wolfram Research. Mathematica. Wolfram Research, Inc., Champaign, Illinois,

version 5.2 edition, 2005.

O. C. Zienkiewicz. Achievements and some unsolved problems of the ﬁnite element

method. International Journal for Numerical Methods in Engineering, 47(1-3):

9–28, 2000.

126

Livros Grátis
( http://www.livrosgratis.com.br )
 
Milhares de Livros para Download:
 
Baixar livros de Administração
Baixar livros de Agronomia
Baixar livros de Arquitetura
Baixar livros de Artes
Baixar livros de Astronomia
Baixar livros de Biologia Geral
Baixar livros de Ciência da Computação
Baixar livros de Ciência da Informação
Baixar livros de Ciência Política
Baixar livros de Ciências da Saúde
Baixar livros de Comunicação
Baixar livros do Conselho Nacional de Educação - CNE
Baixar livros de Defesa civil
Baixar livros de Direito
Baixar livros de Direitos humanos
Baixar livros de Economia
Baixar livros de Economia Doméstica
Baixar livros de Educação
Baixar livros de Educação - Trânsito
Baixar livros de Educação Física
Baixar livros de Engenharia Aeroespacial
Baixar livros de Farmácia
Baixar livros de Filosofia
Baixar livros de Física
Baixar livros de Geociências
Baixar livros de Geografia
Baixar livros de História
Baixar livros de Línguas

Baixar livros de Literatura
Baixar livros de Literatura de Cordel
Baixar livros de Literatura Infantil
Baixar livros de Matemática
Baixar livros de Medicina
Baixar livros de Medicina Veterinária
Baixar livros de Meio Ambiente
Baixar livros de Meteorologia
Baixar Monografias e TCC
Baixar livros Multidisciplinar
Baixar livros de Música
Baixar livros de Psicologia
Baixar livros de Química
Baixar livros de Saúde Coletiva
Baixar livros de Serviço Social
Baixar livros de Sociologia
Baixar livros de Teologia
Baixar livros de Trabalho
Baixar livros de Turismo
 
 