CENTRO FEDERAL DE EDUCAÇÃO TECNOLÓGICA
DE MINAS GERAIS
Diretoria de Pesquisa e Pós-Graduação
Programa de Pós-Graduação em Modelagem Matemática e
Computacional
Modelagem e Simulação Numérica da
Radiação Sonora de um Cilindro Infinito
Pulsante
Aluno: Emerson de Sousa Costa
Orientadora: Profª. Dra. Ester Naves Machado Borges
Co-Orientador:
Prof. Dr. Márcio Matias Afonso
Belo Horizonte, 29 de agosto de 2008.
CENTRO FEDERAL DE EDUCAÇÃO
TECNOLÓGICA DE MINAS GERAIS
Diretoria de Pesquisa e Pós-Graduação
Programa de Pós-Graduação em Modelagem
Matemática e Computacional
Modelagem e Simulação Numérica da
Radiação Sonora de um Cilindro Infinito
Pulsante
Dissertação de Mestrado, submetida ao Programa
de Pós-Graduação em Modelagem Matemática e
Computacional, como parte dos requisitos exigidos
para a obtenção do título de Mestre em Modelagem
Matemática e Computacional.
Aluno: Emerson de Sousa Costa
Orientadora: Profª. Dra. Ester Naves Machado Borges
Co-Orientador:
Prof. Dr. Márcio Matias Afonso
Belo Horizonte, 29 de agosto de 2008.
C837m
Costa, Emerson de Sousa
Modelagem e simulação numérica da radiação sonora de um
cilindro infinito pulsante. – 2008.
72 f.
Orientadora: Ester Naves Machado Borges
Dissertação (mestrado) – Centro Federal de Educação Tecnológica
de Minas Gerais.
1. Método matemático – Teses. 2. Métodos de simulação
I. Borges, Ester Naves Machado. II. Centro Federal de Educação
Tecnológica de Minas Gerais. III. Título.
CDD 519.6
Elaboração da ficha catalográfica por Biblioteca-Campus II / CEFET-MG
Dedicatória
Para meu pai (in memorian), todos meus
familiares e amigos.
i
Agradecimentos
Agradeço, primeiramente, a Deus, por ter me dado força, paciência, persistência,
determinação, conforto e paz necessárias à conclusão de mais uma etapa da minha formação
acadêmica.
A professora Ester Naves Machado Borges pela orientação inestimável, pela amizade
e confiança depositada em mim, além da paciência em todos os momentos de dúvidas.
Ao professor Márcio Matias Afonso pela orientação sempre presente, críticas
construtivas e valiosas contribuições ao meu trabalho.
Aos demais professores pelo apoio na realização deste trabalho e aos colegas de
mestrado pelo companheirismo e amizade.
À minha mãe, minhas irmãs e irmão pelo carinho e pelas palavras de incentivo que
não me deixaram desanimar nunca.
A todos os colegas das escolas, que incentivaram e ajudaram no desenvolvimento
deste trabalho. Em especial à Nedina, Marcílio e Donizete que tiveram tolerância e me
apoiaram em todos os momentos.
Por fim, aos não menos importantes, demais amigos e familiares pelo apoio.
Agradeço de coração a todos vocês.
ii
“O conhecimento é o processo de acumular dados;
a sabedoria reside na sua simplificação."
Martin H. Fischer
iii
Resumo
Neste trabalho é obtida a solução da equação diferencial de onda acústica radiada por um
cilindro infinito pulsante em um meio livre e homogêneo. A solução analítica é encontrada
utilizando-se o Método de Separação de Variáveis e os resultados obtidos são comparados
com os da literatura. A formulação e implementação da equação de onda que é regida pelo
operador diferencial de Helmholtz, são obtidas utilizando-se o Método dos Elementos de
Contorno (MEC). A transformação da equação de Helmholtz em Equação Integral de
Contorno, bem como sua solução, é apresentada de forma detalhada no texto. O Método de
Elementos de Contorno apresenta grande aplicação na solução de determinados problemas
acústicos em campo aberto, em relação aos métodos diferenciais. Este método reduz a
dimensão do problema, simplificando com isso os dados de entrada a serem trabalhados e
reduzindo o tempo computacional utilizado.
Palavras-chave: Radiação; Acústica; Elementos de Contorno; Modelagem.
iv
Abstract
This paper obtained the solution of differential equation of acoustic wave radiated by an
infinite cylinder in a pulsating half free and homogeneous. The analytical solution is found
using the method of separation of variables and results are compared with those of literature.
The formulation and implementation of the wave equation that is governed by the spread of
Helmholtz, are obtained using the Boundary Element method (BEM). The transformation of
the Helmholtz equation in the Boundary Integral equation and its solution is presented in
detail in the text. The Boundary Element Method has a great application in the solution of the
problem of acoustic radiation in open fields, when compared with the differential methods.
The BEM reduces the size of the problem because it simplifies the input data to be worked
and reduces the computational time used.
Keywords: Radiation; acoustics; Boundary Element; Modeling.
v
Conteúdo
INTRODUÇÃO
1
1.1 − Considerações Gerais
1
1.2 − Objetivos
3
1.2.1 - Objetivo Geral
3
1.2.2 - Objetivos Específicos
3
1.3 − Organização do Trabalho
4
EQUAÇÃO DE HELMHOLTZ PARA MEIOS ACÚSTICOS
2.1 - Equação Linear da Onda Acústica
5
5
2.1.1 - Equação de Estado
6
2.1.2 - Equação de Continuidade
8
2.1.3 - Equação de Euler
8
2.1.4 - Equação Linear da Onda
9
2.2 – Equação de Helmholtz
11
2.3 – Solução Analítica da Equação de Helmholtz
12
2.4 – Sumário
22
TÉCNICAS NUMÉRICAS
23
3.1– Técnicas Numéricas Diferenciais
23
3.2 – Técnicas Numéricas Integrais
25
3.3 – Sumário
27
A EQUAÇÃO INTEGRAL DE CONTORNO
28
4.1– Introdução
28
4.2– Equação Integral de Contorno
28
4.3 – Equação Integral no Contorno
32
4.4 – O Coeficiente c, termo livre do sinal de integração
35
4.5 – Tipos de Elementos de Contorno
36
vi
4.6– Discretização e Colocação
40
v
n
46
4.8 – Cálculo dos Elementos não Singulares das Matrizes G e H
48
4.9 – Cálculo dos Elementos Singulares das Matrizes G e H
49
4.10 – Sumário
51
4.7 – Derivada da solução fundamental em relação a normal
RESULTADOS NUMÉRICOS E COMPARAÇÕES
52
5.1 – Resultado Analítico e Comparações
52
5.2 – Resultado Numérico e Comparações
54
5.3 – Sumário
61
CONCLUSÃO GERAL E PROPOSTA PARA TRABALHOS FUTUROS
62
6.1 − Conclusão Geral
62
6.2 − Propostas para trabalhos futuros
63
REFERÊNCIAS BIBLIOGRÁFICAS
64
ANEXO A
66
ANEXO B
69
ANEXO C
70
ANEXO D
71
Função Delta de Dirac e Solução Fundamental
ANEXO E
71
72
vii
Lista de Tabelas
Tabela 5.1 – Comparação de resultados
53
Tabela 5.2 – Resultados para discretização com 8 elementos
55
Tabela 5.3 – Resultados para discretização com 16 elementos
56
Tabela 5.4 – Resultados para discretização com 32 elementos
57
Tabela 5.5 – Resultados para discretização com 64 elementos
58
Tabela 5.6 – Resultados para discretização com 128 elementos
59
Tabela 5.7 – Resultados para discretização com 256 elementos
60
Tabela B1 – Valores de xi e wi para integração gaussiana, utilizando-se 10
pontos de Gauss
69
Tabela E.1 – Nós e respectivos elementos
72
viii
Lista de Figuras
Figura 2.1 – Cilindro infinito de raio r = a
13
Figura 2.2 – As coordenadas cilíndricas (r, ψ, z)
14
Figura 4.1 – Representação do domínio e dos contornos
30
Figura 4.2 – Representação do domínio sobre um contorno 2D
33
Figura 4.3 – Malhas compostas de diferentes tipos de elementos
38
Figura 4.4 – Malhas com o mesmo número de elementos
38
Figura 4.5 – Aproximação realizada por elementos contínuos e descontínuos
39
Figura 4.6 (a) – Elementos de contorno: elementos constantes
39
Figura 4.6 (b) – Elementos de contorno: elementos lineares
39
Figura 4.6 (c) – Elementos de contorno: elementos quadráticos
40
Figura 4.7 – Contorno discretizado
41
Figura 4.8 – Elemento Linear
45
Figura 4.9 – Sistema de coordenadas no elemento
45
Figura 5.1 – Comparação de resultados implementados
53
Figura 5.2 – Ilustração bidimensional da discretização do cilindro em 12
elementos
Figura 5.3 – Representação gráfica dos resultados obtidos com a discretização
com 8 elementos
Figura 5.4 – Representação gráfica dos resultados obtidos com a discretização
com 16 elementos
Figura 5.5 – Representação gráfica dos resultados obtidos com a discretização
com 32 elementos
Figura 5.6 – Representação gráfica dos resultados obtidos com a discretização
com 64 elementos
Figura 5.7 – Representação gráfica dos resultados obtidos com a discretização
com 128 elementos
Figura 5.8 – Representação gráfica dos resultados obtidos com a discretização
com 256 elementos
ix
54
55
56
57
58
59
60
Simbologia
CARACTERES LATINOS
Variável
Descrição
Unidade
s
adensamento em um ponto
R
constante dos gases
F
força
β
módulo adiabático (coeficiente de expansão volumétrica do fluido)
r
r
posição de equilíbrio de uma partícula de fluido em (x, y, z)
[m]
p
pressão acústica em um ponto
[Pa]
P0
pressão de equilíbrio no fluido
[Pa]
P
pressão instantânea em um ponto
[Pa]
r
raio do cilindro
[m]
representam vetores unitários nas direções x, y e z, respectivamente.
[m]
T
temperatura em graus Celsius
[°C]
TK
temperatura em Kelvin
[K]
t
tempo
[s]
r
u
velocidade da partícula de fluido
[m/s]
c
velocidade de fase da onda
[m/s]
xˆ, yˆ , e zˆ
[J/kgK]
[N]
x
-
CARACTERES GREGOS
Variável
Descrição
Dimensão
ρ0
densidade de equilíbrio do fluido
[kg/m3]
ρ
densidade instantânea em um ponto
[kg/m3]
r
deslocamento da partícula de fluido em relação à
ξ
posição de equilíbrio
[m]
ω
freqüência angular
[rad/s]
φ
potencial de velocidade
[m/s]
OPERADORES
Variável
∂ / ∂i
Descrição
derivada parcial
∇ ⋅ (*) = div(*)
divergente
∇(*) = grad (*)
gradiente
∇ 2 (*)
*
Laplaciano
norma de um vetor ou matriz
xi
RELAÇÕES ENTRE CARACTERES
Variável
s=
ρ − ρ0
ρ0
r
ξ = ξ x xˆ + ξ y yˆ + ξ z zˆ
r
r = xxˆ + yyˆ + zzˆ
Descrição
adensamento em um ponto
deslocamento da partícula de fluido em relação à posição de
equilíbrio
posição de equilíbrio de uma partícula de fluido em (x, y, z)
p = P − P0
pressão acústica em um ponto
T + 273.15 = Tk
temperatura em graus Celsius
r
r ∂ξ
r
u=
= u x xˆ + u y yˆ + u z z
∂t
velocidade de partícula
LISTA DE ABREVIATURAS
Abreviação
EI
Descrição
Equação Integral
EIC
Equações Integrais de Contorno
MDF
Método das Diferenças Finitas
MEC
Método dos Elementos de Contorno
MEF
Método dos Elementos Finitos
xii
Capítulo 1
Introdução
_________________________
1.1 − Considerações Gerais
O som pode ser definido como uma variação da pressão ambiente detectável pelo
sistema auditivo e o ruído como um som sem harmonia. Um mecanismo bastante comum para
gerar sons consiste em fazer vibrar uma estrutura. Estruturas vibrantes movimentam
ciclicamente as moléculas do ar ao seu redor, gerando localmente regiões de concentração e
de rarefação destas, o que provoca variações de pressão.
A propagação sonora ao ar livre é normalmente estudada em termos de três
componentes: a fonte sonora, a trajetória de transmissão e o receptor. Primeiramente, a fonte
emite uma certa potência sonora, gerando um nível sonoro que pode ser medido nas
proximidades da fonte. A partir daí, o nível sonoro é atenuado à medida que o som se
propaga, entre a fonte e o receptor, ao longo de determinada trajetória.
A modelagem da radiação acústica é de fundamental importância para se compreender
a propagação das ondas acústicas e, conseqüentemente, desenvolver mecanismos para
atenuação de ruídos acústicos. Para estimativas de níveis de pressão sonora, em certas
ocasiões, é preciso conhecer os níveis de potência sonora das fontes em questão. É este o
caso, por exemplo, quando se deseja determinar o nível de pressão sonora gerado pelo
maquinário de um ambiente industrial e o nível de pressão sonora devido ao tráfego de uma
rodovia, entre outros.
1
Ondas acústicas em fluidos não viscosos, como ar e água, são ondas longitudinais, isto
é, as moléculas do fluido se movem para frente e para trás na direção de propagação da onda,
produzindo regiões adjacentes de compressão e de rarefação.
Na análise que se segue, são feitas algumas suposições necessárias para o estudo da
propagação do som em fluidos na sua forma mais simples. Dentre elas, tem-se que os efeitos
das forças gravitacionais são considerados desprezíveis; o fluido é homogêneo, isotrópico e
perfeitamente elástico e as ondas são de amplitude relativamente pequena, de modo que as
variações na densidade do meio são pequenas quando comparadas com seu valor de
equilíbrio.
O termo partícula de fluido é utilizado para identificar um elemento de volume, grande
o suficiente, para conter um número muito grande de moléculas, de forma que este possa ser
considerado como um meio contínuo, contudo pequeno o bastante para que todas as variáveis
acústicas possam ser consideradas constantes ao longo deste elemento de volume.
Neste trabalho é mostrada a solução analítica da equação de onda para um cilindro
infinito que está vibrando (expandindo e contraindo) uniformemente na direção radial com
amplitude constante. Essa solução encontrada é então comparada com a literatura. Uma
formulação numérica para este problema também é apresentada, utilizando-se a forma direta
do Método dos Elementos de Contorno.
Atualmente, o Método de Elementos de Contorno é um dos métodos mais avançados e
é utilizado especialmente para o tratamento de problemas envolvendo meios infinitos e semiinfinitos. Uma das grandes vantagens desse método é que ele permite a redução da dimensão
do problema, diminuindo o número de equações utilizadas, permitindo a solução apenas no
contorno, sem necessidade de se analisar todo o domínio.
2
1.2 − Objetivos
1.2.1 - Objetivo Geral
O objetivo deste trabalho é fazer uma modelagem analítica e numérica de um cilindro
infinito pulsante que está vibrando uniformemente na direção radial com amplitude constante.
A modelagem numérica deste problema é feita a partir da formulação direta do
Método dos Elementos de Contorno.
1.2.2 - Objetivos Específicos
● Realizar um estudo teórico e numérico sobre a equação de onda radiada por um cilindro
infinito pulsante.
● Comparar os resultados obtidos na implementação, utilizando-se o Método de Elementos de
Contorno, aos resultados de Yoon, (1990) e Papini, (1999).
3
1.3 − Organização do Trabalho
A estrutura do trabalho apresenta, inicialmente, no Capítulo 2, a formulação analítica
da radiação acústica de um cilindro infinito pulsante. Logo após é descrita a implementação
da formulação analítica.
No Capítulo 3 é apresentada uma breve introdução dos métodos numéricos que são
comumente utilizados.
No Capítulo 4 são desenvolvidos os conceitos básicos de Equações Integrais e do
Método de Elementos de Contorno, empregados neste trabalho.
No Capítulo 5 são apresentados os resultados analíticos e numéricos e estes são
comparados com os resultados obtidos na literatura.
No Capítulo 6 são apresentadas a conclusão geral e perspectivas de futuros trabalhos.
4
Capítulo 2
Equação de Helmholtz para Meios
Acústicos
______________________________
Neste capítulo é apresentada a dedução da Equação de Helmholtz para o fenômeno de
radiação acústica. Esta equação será obtida a partir da Equação Linear da Onda.
2.1 - Equação Linear da Onda Acústica
A equação que descreve o fenômeno da radiação acústica é obtida utilizando-se as
equações de estado, de conservação da massa e da conservação da quantidade de movimento,
na presença de um sinal acústico, segundo Ziomek, (1995). Para obter a formulação analítica
da radiação acústica, devem ser feitas algumas suposições:
•
Os efeitos das forças gravitacionais são desprezíveis;
•
O fluido é homogêneo, isotrópico e perfeitamente elástico;
•
Efeitos dissipativos devido à viscosidade ou condução de calor não estão presentes;
•
A amplitude de onda sonora é relativamente pequena, de modo que as variações na
densidade do meio são pequenas quando comparadas com seu valor de equilíbrio.
Estas suposições são necessárias para se estudar a propagação de som em fluidos na
sua forma mais simples. Experimentos mostram que uma análise deste tipo descreve
adequadamente os fenômenos acústicos mais comuns.
5
2.1.1 - Equação de Estado
Para meios fluidos, a equação de estado relaciona grandezas físicas que descrevem o
comportamento termodinâmico de um fluido. Por exemplo, para um gás perfeito tem-se
P = ρ r TK ,
(2.1)
que relaciona a pressão instantânea P (em Pascal-Pa), a densidade instantânea ρ (kg/m3) e a
temperatura absoluta TK (Kelvin). A grandeza r =  R  é a constante específica do gás e
M 
depende da constante universal dos gases R e do seu peso molecular M. Para o ar, R ≈ 287
J/kg.K.
Processos acústicos são aproximadamente isentrópicos (adiabáticos e reversíveis). A
condutividade térmica do fluido e os gradientes de temperatura da perturbação (onda acústica)
são pequenos de forma que a troca de energia térmica entre camadas adjacentes do fluido
pode ser desprezada. Nessas condições, a entropia do fluido permanece aproximadamente
constante. O comportamento acústico de um gás sob essas condições é descrito pela equação
adiabática de estado descrita por
γ
P  ρ 
=  ,
P0  ρ 0 
(2.2)
onde γ é a razão entre os valores específicos à pressão e volume constantes. Para que as
perturbações acústicas do fluido sejam consideradas adiabáticas, não pode haver troca de
energia térmica entre elementos adjacentes do fluido. Isto significa que a condutividade
térmica do fluido e que os gradientes de temperatura da perturbação devem ser pequenos o
bastante para que não ocorra fluxo térmico significante durante o tempo da perturbação. Isso
ocorre nas freqüências e amplitudes de interesse em acústica, Kinsler et al., (1982).
Se o fluido não se comporta como um gás perfeito, a equação que descreve seu
comportamento é mais complicada. Neste caso, a relação isentrópica entre pressão e
6
densidade de flutuação é determinada experimentalmente. Esta relação pode ser representada
utilizando-se uma expansão em série de Taylor
 ∂P 
1  ∂2P 
2
P = P0 +   (ρ − ρ 0 ) +  2  (ρ − ρ 0 ) + ... .
2  ∂ρ  ρ
 ∂ρ  ρ 0
0
(2.3)
Na expressão acima, as derivadas parciais são constantes determinadas para
compressão adiabática e para a expansão do fluido em torno de sua densidade de equilíbrio.
Se tais flutuações são pequenas, somente o termo de ordem mais baixa em (ρ − ρ 0 ) é
considerado. Assim, obtém-se uma relação linear entre a flutuação de pressão e variação de
densidade
P − P0 = β
(ρ − ρ 0 )
ρ0
,
(2.4)
onde β é o módulo adiabático ou coeficiente de expansão volumétrica do fluido, dado por
 ∂P 
 .
 ∂ρ  ρ 0
β = ρ 0 
(2.5)
Pode-se definir o adensamento s em um ponto, como a razão entre a variação de
densidade e seu valor de equilíbrio, s =
ρ − ρ0
e a pressão acústica p como a variação da
ρ0
pressão em relação a seu valor de equilíbrio, p = P − P0 . Dessa maneira pode-se reescrever a
equação (2.4) em termos da pressão acústica p e do adensamento s,
p≈β s.
(2.6)
A restrição essencial é que o adensamento s deve ser muito pequena, s << 1 , Kinsler
et al., (1982).
7
2.1.2 - Equação de Continuidade
Para relacionar o movimento do fluido com sua compressão ou expansão, precisa-se
r
de uma função que relacione a velocidade u da partícula do fluido com sua densidade
instantânea ρ .
Considera-se um fenômeno de transporte de massa através de um elemento de volume
infinitesimal de fluido, fixo no espaço. A equação de continuidade relaciona a taxa de
crescimento ou decrescimento de massa no interior desse elemento de volume com o fluxo de
massa através da superfície fechada que o envolve e tem a seguinte expressão
r
∂ρ
+ ∇.(ρ u ) = 0 .
∂t
Como a densidade instantânea
(2.7)
ρ
pode ser expressa em função do adensamento
ρ = ρ 0 (1 + s) , pode-se usar o fato que ρ 0 é constante no espaço e no tempo, e que s é muito
pequena, isto é, s << 1 , para expressar a equação acima da seguinte maneira:
r
∂s
+ ∇⋅u = 0,
∂t
(2.8)
que é a equação da continuidade linearizada.
2.1.3 - Equação de Euler
Para fluidos reais, a existência de viscosidade e o fato de que processos acústicos não
são perfeitamente adiabáticos introduzem termos dissipativos. Entretanto, uma vez que os
efeitos da condutividade térmica na equação de estado foram considerados desprezíveis,
pode-se também ignorar os efeitos da viscosidade e considerar o fluido como sendo não
viscoso.
8
A equação da conservação da taxa de variação de quantidade de movimento relaciona
r
a pressão acústica p com a velocidade u instantânea da partícula, para um fluido adiabático e
não viscoso, isto é, os efeitos da viscosidade do fluido podem ser desprezados.
Dessa maneira obtém-se a equação de Euler (equação de força) que é a equação de
movimento para fluidos invíscidos:
r
∂u
= −∇p .
ρ0
∂t
(2.9)
2.1.4 - Equação Linear da Onda
Aplicando-se o operador divergente em ambos os lados da equação (2.9), obtém-se
r
 ∂u 
2
 = −∇ p ,
 ∂t 
ρ0 ∇ ⋅ 
(2.10)
onde ∇2 é o operador Laplaciano.
Derivando-se a equação (2.8) em relação ao tempo e usando o fato de que
v
v
∂(∇ ⋅ u )
 ∂u 
= ∇ ⋅  ,
∂t
 ∂t 
(2.10-a)
obtém-se
r
∂2s
 ∂u 
+ ∇ ⋅  = 0.
∂t 2
 ∂t 
(2.11)
As equações (2.10) e (2.11) podem ser combinadas numa única equação:
∇2 p = ρ0
∂2s
.
∂t 2
(2.12)
9
Usando-se a equação de estado (2.6) para eliminar o adensamento s , obtém-se
∇2 p =
1 ∂2 p
,
c 2 ∂t 2
(2.13)
onde a constante c , definida como
c=
β
,
ρ0
(2.14)
é denominada velocidade de fase (propagação) da onda acústica no meio.
A equação (2.13) é a equação linear de onda para a propagação sonora em meios
homogêneos e sem perdas.
r
Para fluidos não viscosos, a velocidade da partícula é irrotacional, ∇ × u = 0 . Isto
significa que ela pode ser expressa como o gradiente de uma função escalar φ , denominada
potencial de velocidade,
r
u = ∇φ .
(2.15)
O significado físico deste importante resultado é que a excitação acústica de um fluido
invíscido não envolve fluxo rotacional: efeitos como tensões cisalhamento ou turbulência não
estão presentes. Em fluidos reais, para os quais a viscosidade é finita, a velocidade de
partícula não é irrotacional em todos os pontos do fluido. Na maioria dos processos acústicos,
a presença de pequenos efeitos rotacionais limita-se à vizinhança ao redor dos contornos e
exerce pouca influência sobre a propagação do som.
Substituindo-se a equação (2.15) na equação (2.9), obtém-se
 ∂φ

∇  ρ0
+ p = 0 .
∂t


(2.16)
10
Verifica-se que a expressão entre parênteses na equação (2.16) pode ser escolhida
nula, caso não haja excitação acústica, Kinsler et al., (1982). Dessa forma, tem-se que
p = −ρ 0
∂φ
.
∂t
(2.17)
Substituindo-se a equação (2.17) na equação (2.10), obtém-se a equação de onda
linearizada, expressa em termos do potencial de velocidade de onda acústica:
∇ 2φ =
1 ∂ 2φ
.
c 2 ∂t 2
(2.18)
2.2 – Equação de Helmholtz
Para se obter a solução da equação (2.18), supõe-se que o potencial de velocidade
v
φ (t , r ) tem dependência harmônica no tempo e pode ser escrito da seguinte forma
r
r
φ (t , r ) = φ f (r ) eiωt .
(2.19)
Na equação (2.19), φ f representa a parte espacial do potencial de velocidade e ω
representa a freqüência angular da perturbação. Substituindo-se a equação (2.19) na equação
(2.18) obtém-se a Equação de Helmholtz linear, tridimensional, homogênea, em coordenadas
cilíndricas, independente do tempo para um meio sem perdas, expressa em termos do
potencial de velocidade da onda acústica
()
()
∇ 2φ f r + κ 2φ f r = 0 ,
(2.20)
onde
κ=
2πf 2π
=
,
c
λ
(2.21)
é o número de onda e
11
c = fλ .
(2.22)
Para meios homogêneos, o valor de c (velocidade de fase) é constante.
2.3 – Solução Analítica da Equação de Helmholtz
Neste trabalho, considera-se um cilindro de raio r = a cuja superfície está vibrando
uniformemente na direção radial com uma amplitude de Vo metros por segundo, numa
freqüência de f Hertz, representado na Figura 2.1.
Figura 2.1 – Cilindro infinito de raio r = a
Considera-se ainda que o problema obedece a condição de contorno de Neumann
∂φ
v
= V0 = u (r ), em r = a .
∂r
(2.23)
Na Equação de Helmholtz (2.20), utiliza-se o operador Laplaciano ∇2 em coordenadas
r
cilíndricas (r, ψ, z), devido à simetria do problema e dessa forma, φ f (r ) é o potencial de
velocidade no ponto (r ,ψ , z ) .
12
Fig. 4 – As coordenadas cilíndricas (r, ψ, z)
Figura 2.2 – As coordenadas cilíndricas (r, ψ, z)
A solução da equação (2.20) é obtida utilizando-se o método de separação de
variáveis, isto é, supondo-se que a solução pode ser escrita da seguinte forma:
v
φ f (r ) = φ f (r ,ψ , z ) = R(r )Ψ (ψ )Z ( z ) .
(2.24)
Substituindo a equação (2.24) na equação (2.20) obtém-se
1 d2
1 d
1
d2
1 d2
2
R
(
r
)
+
R
(
r
)
+
Ψ
(
ψ
)
+
k
=
−
Z (z ) .
R (r ) dr 2
rR (r ) dr
r 2 Ψ (ψ ) dψ 2
Z ( z ) dz 2
(2.25)
Como o lado esquerdo da equação (2.25) é função de r e ψ, e o lado direito é função
de z, esta igualdade só se verifica se ambos os lados forem iguais a uma constante, isto é,
−
1 d2
Z ( z ) = k z2 ,
Z ( z ) dz 2
(2.26)
onde k z2 é chamada de constante de separação. A equação (2.26) pode ser reescrita como
d2
Z ( z ) + k z2 Z (z ) = 0 ,
dz 2
(2.27)
13
que é uma equação diferencial ordinária de segunda ordem, homogênea, que tem como
solução exata
Z ( z ) = AZ e −ik z z + BZ e ik z z ,
(2.28)
onde AZ e BZ são constantes em geral complexas, cujos valores são determinados a partir das
condições de contorno. Se k z é positivo o primeiro termo na equação acima representa uma
onda plana viajando no sentido positivo de z, enquanto o segundo termo representa uma onda
plana viajando no sentido negativo de z.
Retornando à equação (2.25), observamos que o lado esquerdo desta deve também ser
igual a k z2 . Assim,
1 d2
1 d
1
d2
(
)
(
)
R
r
+
R
r
+
Ψ (ψ ) + k 2 = k z2 .
R(r ) dr 2
rR(r ) dr
r 2 Ψ (ψ ) dψ 2
(2.29)
Multiplicando-se por r 2 ambos os lados da equação acima, obtém-se
r2 d 2
r d
1
d2
(
)
(
)
R
r
+
R
r
+
Ψ (ψ ) + k 2 r 2 = r 2 k z2 ,
R(r ) dr 2
R (r ) dr
Ψ (ψ ) dψ 2
(2.30)
r2 d 2
r d
1
d2
2
2 2
R
(
r
)
+
R
(
r
)
+
k
−
k
r
=
−
Ψ (ψ ) .
z
Ψ (ψ ) dψ 2
R(r ) dr 2
R (r ) dr
(2.31)
ou
(
)
Fazendo-se (k 2 − k z2 ) = k r2 encontra-se
r2 d 2
r d
1
d2
2 2
R
(
r
)
+
R
(
r
)
+
k
r
=
−
Ψ (ψ ) .
r
R (r ) dr 2
R(r ) dr
Ψ (ψ ) dψ 2
(2.32)
14
Como o lado esquerdo da equação acima é função somente de r e o lado direito é
função de ψ, esta igualdade só se verifica se ambos os lados forem iguais a uma constante,
isto é, se
−
d2
Ψ (ψ ) = n 2 .
Ψ (ψ ) dψ 2
1
(2.33)
Na equação acima n2 é chamada de constante de separação. Multiplicando-se ambos
os lados da equação (2.33) por Ψ (ψ ) obtém-se
d2
Ψ (ψ ) + n 2 Ψ (ψ ) = 0 ,
dψ 2
(2.34)
que é uma equação diferencial ordinária de segunda ordem, homogênea, que tem como
solução exata
Ψ (ψ ) = Aψ e − inψ + Bψ e inψ ,
(2.35)
ou, utilizando-se a fórmula de Euler,
Ψ (ψ ) = Aψ cos (nψ ) + Bψ sen (nψ ) .
Na equação (2.36)
(2.36)
Aψ e Bψ são constantes em geral complexas, cujos valores são
determinados impondo-se as condições de contorno.
Finalmente, o lado esquerdo da equação (2.32), deve ser igual a n2. Com isso,
r2 d 2
r d
R (r ) +
R (r ) + kr2 r 2 = n 2 .
2
R (r ) dr
R (r ) dr
(2.37)
Multiplicando ambos os lados da equação acima por R(2r ) obtém-se
r
15
 2 n2 
d2
1 d
(
)
(
)
R
r
+
R
r
+
kr − 2  R(r ) = 0 ,
dr 2
r dr
r 

(
(2.38)
)
onde k r2 = k 2 − k z2 .
O próximo passo é transformar a equação (2.38) em uma equação diferencial de
Bessel, que tem solução exata conhecida.
Seja R(r) uma função g (k r r ) tal que
R(r ) = g (k r r ) .
(2.39)
Desta forma,
dg (k r r )
d
.
R(r ) =
dr
dr
(2.40)
Fazendo-se k r r = u tem-se que d = d du . Assim,
dr
du dr
d
dg (u ) du
,
R (r ) =
dr
du dr
(2.41)
dg (k r r ) d
dg (k r r ) .
d
(k r r ) = k r
R (r ) =
(
)
dr
d k r r dr
d (k r r )
(2.42)
dg (k r r )
d
,
R(r ) = k r
dr
d (k r r )
(2.43)
ou
Isto é,
e também
16
2
d2
d d
d  dg (u )
d  dg 
d  dg 
2 d g (k r r )
R
(
r
)
=
R
(
r
)
=
k
=
k
=
k
k
=
k
,
r
r
r
r
r
2
dr dr
dr 
du 
du  dr 
du  du 
dr 2
d 2 (k r r )
(2.44)
2
d2
2 d g (k r r )
.
R
r
k
(
)
=
r
2
dr 2
d (k r r )
(2.45)
ou
Substituindo-se na equação (2.38) as equações (2.43), (2.44) e (2.45) na equação (2.38)
obtém-se
k r2
d 2g k r dg  2 n 2 
+
+  k r − 2  g = 0 .
r 
du 2 r du 
(2.46)
Dividindo-se a equação acima por k r2 obtém-se
d 2g
1 dg 
n2 

g = 0 .
+
+
1
−
du 2 k r r du 
k r2 r 2 
(2.47)
Como u = k r r verifica-se que a equação (2.47) é uma equação diferencial de Bessel
d 2g 1 dg 
n2

+
+ 1−
du 2 u du  u 2

 g = 0 .

(2.48)
A solução exata da equação (2.48), para valores arbitrários de n é dada por:
g (u ) = Ar J n (u ) + Br Yn (u ) ,
(2.49)
g (u ) = Ar H n(1) (u ) + B r H n( 2 ) (u ) ,
(2.50)
ou
17
onde J n (u ) e Yn (u ) são as funções de Bessel da 1ª e 2ª espécies e de n-ésima ordem, e H n(1) (u ) e
H n( 2 ) (u ) são as funções de Hankel da 1ª e 2ª espécies e de n-ésima ordem, respectivamente,
também conhecidas como funções de Bessel da 3ª espécie, e são dadas por
H n(1) (u ) = J n (u ) + iYn (u ) ,
(2.51)
H n( 2) (u ) = J n (u ) − iYn (u ) .
(2.52)
Utilizando-se as equações (2.46) e (2.50), pode-se finalmente escrever a solução exata
da equação (2.38) para valores arbitrários de n,
R (r ) = Ar J n (k r r ) + B r Yn (k r r ) ,
(2.53)
ou
R (r ) = Ar H
(1)
n
(kr r ) + Br H ( 2n) (kr r ) .
(2.54)
As constantes Ar e Br , em geral complexas, são determinadas pelas condições de
contorno. Substituindo-se as equações (2.28), (2.36) e (2.54) na equação (2.24) obtém-se
φ f (r ,ψ , z ) = [Ar H n(1) (k r r ) + B r H n( 2 ) (k r r )][Aψ cos(nψ ) + Bψ sen(nψ )] [AZ e −ik z + B Z e ik z ]. (2.55)
z
z
r
Na equação (2.55), k r e k z são as componentes do vetor de propagação k e satisfazem a
relação kr2 + k z2 = k 2 .
Substituindo-se a equação (2.55) na equação (2.19) obtém-se a solução harmônica no
tempo da equação de onda (2.18).
Se o fluido envolve todo o eixo z, isto é, se o fluido envolve todo o cilindro, então
Ψ (ψ )
deve ser uma função periódica com período igual a 2π, e tem-se
18
φ f (r ,ψ , z ) = φ f (r ,ψ + 2π , z ) ,
(2.56)
e os valores de n devem ser inteiros, isto é, n = 0 , ± 1, ± 2, .... Se n = 0, utilizando-se a
equação (2.36) tem-se
Ψ (ψ ) = Aψ ,
n = 0.
(2.57)
Como resultado, o potencial de velocidade φ f (r ,ψ , z ) dado pela equação (2.55) não é função
do ângulo azimutal ψ . Este caso é chamado de axissimétrico.
Para pontos no campo distante, ou seja, pontos que satisfazem a relação
k r r >>
n2
2
,
(2.58)
podemos escrever as seguintes aproximações para as funções de Hankel
H n(1) (k r r ) ≈
 
2
nπ π  
exp + i k r r −
− ,
2
4 
π kr r
 
(2.59)
e
H n( 2 ) (kr r ) ≈
 
2
nπ π 
exp − i kr r −
−  ,
2
4 
π kr r
 
(2.60)
para valores arbitrários de n.
Neste trabalho foi feita a suposição de que o potencial de velocidade é expresso na
forma descrita pela equação (2.19) e dessa forma, as equações (2.59) e (2.60) representam
ondas que se propagam no sentido de decrescimento e crescimento da coordenada cilíndrica r,
respectivamente.
19
Como estamos buscando a solução para uma superfície cilíndrica que está irradiando
no espaço livre, não existem ondas refletidas e portanto devemos escolher Ar = 0 na equação
(2.55).
Dessa maneira, para um cilindro pulsante, o potencial de velocidade é dado por
φ f ,n (r,ψ , z ) = H n(2) (kr r ) [An cos(nψ ) + Bn sen(nψ )] [AZ e −ik z + BZ eik z ].
z
(2.61)
z
O subscrito n foi adicionado a φ f (r ,ψ , z ) para enfatizar que o potencial de velocidade é
função do número inteiro n.
Como a soma das possíveis soluções também é uma solução, temos
∞
φ f (r ,ψ , z ) = ∑ φ f , n (r ,ψ , z ) ,
(2.62)
n=0
isto é,
∞
[
]
φ f (r ,ψ , z ) = ∑ H ( 2n) (k r r )[An cos(nψ ) + Bn sen(nψ )] AZ e −ik z + BZ e ik z ,
z
z
(2.63)
n =0
onde k r2 + k z2 = k 2 .
A equação (2.63) representa a parte espacial do potencial de velocidade do campo
acústico irradiado e os valores das constantes An , Bn , Az e Bz são determinados pela
aplicação das condições de contorno, equação 2.23, na superfície do cilindro vibrante.
Fazendo-se a suposição de que a parte espacial do vetor velocidade da superfície do
cilindro pulsante tem apenas componente radial e que sua dependência no tempo é harmônica,
pode-se escrever
v r (t , a ,ψ , z ) = v f ,r (a,ψ , z ) exp (+ i 2π f t ) ,
(2.64)
20
onde v f ,r (a,ψ , z ) é uma função conhecida de ψ e z . Isto é equivalente a conhecer a
componente radial da velocidade do fluido,
u r , na superfície do cilindro, ou seja, em r = a.
Assim,
ur (t , a,ψ , z ) = vr (t , a,ψ , z ) =
∂
φ (t , r ,ψ , z ) .
∂r
r =a
(2.65)
No presente trabalho, assume-se que o raio do cilindro é unitário e que a velocidade do
cilindro vibrante pode ser escrita da seguinte forma
vr (t , a,ψ , z ) = V0 exp(+ i 2π f t ) ,
(2.66)
onde V0 = 1 .
Dessa forma, a condição de contorno descrita na equação (2.23) pode ser reescrita:
∂φ
= V0 = 1 ,
∂r
(2.67)
em r = a = 1.
É importante lembrar que a velocidade acústica do fluido pode ser expressa da
seguinte forma
r r
r
u (t , r ) = ∇φ (t , r ) ,
(2.68)
e tem-se ainda que
r
r r
r
u r (t , r ) = rˆ . u (t , r ) = rˆ . ∇φ (t , r ) .
(2.69)
Portanto em r = a = 1, tem-se
21
u r (t ,1,ψ , z ) = vr (t ,1,ψ , z ) =
∂
φ (t , r ,ψ , z ) .
∂r
r =1
(2.70)
A equação (2.66) descreve o caso especial de radiação uniforme, isto é, a superfície do
cilindro está vibrando uniformemente na direção radial com uma amplitude de V0 metros por
segundo, numa freqüência de f Hertz. Este caso é também conhecido como modo monopolo
de vibração. Substituindo-se as equações (2.19), (2.55) e (2.64) na equação (2.70), obtém-se
φ (r ,ψ , z ) = −
V0 H 0( 2) (kr )
,
k H 1( 2 ) (ka )
r ≥ a,
(2.71)
que é a parte espacial da solução da equação de Helmholtz (2.20).
É importante ressaltar que essa solução também pode ser aplicada com boa
aproximação para o campo acústico gerado por um cilindro de comprimento finito, sempre
que os efeitos da radiação de suas extremidades possam ser desprezados. Este é, por exemplo,
o caso da radiação de um transdutor eletroacústico de forma cilíndrica e comprimento finito,
Ziomek, (1994).
2.4 – Sumário
Neste capítulo foi obtida a equação que governa o fenômeno da radiação acústica,
utilizando-se as equações de estado, de conservação da massa, e da conservação da
quantidade de movimento. Foram feitas algumas suposições necessárias para se estudar a
propagação de ondas sonoras em fluidos na sua forma mais simples.
A partir da equação de onda (2.18) foi obtida a equação de Helmholtz (2.20),
independente do tempo para um meio sem perdas. Para a solução desta equação, utilizou-se o
método de separação de variáveis. A seguir é apresentada uma breve descrição das principais
técnicas numéricas que podem ser utilizadas a fim de determinar a solução numérica para essa
equação de Helmholtz.
22
Capítulo 3
Técnicas Numéricas
_________________________
Este capítulo descreve algumas técnicas numéricas utilizadas para solucionar
problemas físicos.
3.1– Técnicas Numéricas Diferenciais
As técnicas numéricas diferenciais, também conhecidas como métodos de domínio,
são eficazes na solução de problemas de contorno em domínios fechados, preenchidos por
materiais heterogêneos, não-lineares ou anisotrópicos. Os Métodos de Diferenças Finitas e
Elementos Finitos constituem os maiores representantes desta classe, Afonso,( 2003).
Método de Diferenças Finitas
O Método de Diferenças Finitas - MDF - é uma técnica de solução aproximada para
solucionar equações diferenciais parciais. O Método de Diferenças Finitas é, talvez a mais
antiga técnica de solução numérica aplicada na resolução de equações diferenciais parciais, e
devido à sua simplicidade, ainda é muito utilizado atualmente. Baseia-se na substituição da
operação de diferenciação por uma simples operação de subtração de pontos considerados.
Dessa forma, é possível substituir uma equação contínua, com infinitos graus de liberdade,1
por uma equação discretizada, com número finito e regular de nós. Utilizando-se este
processo, a equação diferencial parcial original é transformada em um conjunto de equações
algébricas, e a solução simultânea desse sistema de equações fornece a solução aproximada da
equação original do problema de contorno, Afonso, (2003).
1
Incógnitas.
23
O Método de Diferenças Finitas é de simples implementação computacional. Além
disso, é capaz de tratar problemas não lineares e anisotrópicos. Entretanto, este método possui
duas limitações muito sérias:
i.
a obrigatoriedade de uma malha regular, o que não permite uma boa modelagem dos
campos cujos gradientes são intensos ou a modelagem correta de problemas que
possuem superfícies curvas;
ii.
dificuldade em representar os campos na interface entre meios diferentes. Devido a
esses inconvenientes, o método das diferenças finitas não foi escolhido para solucionar
os problemas abordados neste trabalho.
Método de Elementos Finitos
O Método de Elementos Finitos –MEF– é uma técnica numérica de solução
aproximada para equações diferenciais. Este método começou a ser utilizado para solucionar
problemas da mecânica de estruturas e expandiu-se rapidamente para outras áreas. O princípio
do Método de Elementos Finitos consiste em dividir o domínio do problema em pequenos
subdomínios, com forma e comprimentos arbitrários, designados por elementos. Esse
procedimento permite que o MEF modele, precisamente, objetos cuja forma geométrica seja
complexa. Além disso, a densidade de elementos pode ser ajustada de acordo com a
necessidade de cada problema. No interior de cada elemento, a incógnita é aproximada por
uma função de interpolação e, utilizando-se o método dos resíduos ponderados ou método
variacional, a equação diferencial parcial é transformada em um sistema algébrico de
equações, no qual a matriz dos coeficientes é esparsa e, em alguns casos também é simétrica.
A grande vantagem do MEF reside na sua flexibilidade. Em particular, destacam-se sua
capacidade em modelar problemas com geometrias complexas e cujo domínio esteja
preenchido por diferentes materiais, Afonso, (2003).
24
3.2 – Técnicas Numéricas Integrais
Método de Elementos de Contorno
A redução das dimensões do problema analisado constitui uma das razões pela qual o
Método de Elementos de Contorno é tão atrativo: em problemas bidimensionais, apenas o
contorno unidimensional do domínio necessita ser discretizado em elementos; em problemas
tridimensionais, apenas as superfícies do contorno necessitam ser discretizadas.
As soluções nos pontos internos são calculadas após as incógnitas de contorno terem
sido calculadas, de maneira semelhante a um pós-processamento. A densidade, distribuição e
localização dos pontos internos não interferem na malha de contorno, tampouco nos valores
das incógnitas de contorno.
Como outra característica do Método de Elementos de Contorno tem-se que as
condições
de
contorno
para
domínios
infinitos
e
semi-infinitos
são
satisfeitas
automaticamente, o que elimina a necessidade de discretização numérica em contornos
remotos.
As origens do método em discussão estão diretamente ligadas às equações integrais
que durante os quinze primeiros anos da sua aplicação moderna, a partir de cerca de 1960 até
1975, utilizou-se a nomenclatura "Equação Integral" quase exclusivamente, não importando
qual forma utilizada, direta ou indireta.
No caso do método indireto, as grandezas matemáticas utilizadas para a formulação do
problema não são as variáveis físicas do problema original. Por sua vez, no método direto, as
variáveis originais presentes no contorno do problema físico são utilizadas na formulação das
EIC, Ciskowski, (1991).
Considerando ainda o período de 1960 a 1975, o Método de Equações Integrais de
Contorno foi intensamente aplicado em duas grandes aplicações da mecânica, quais sejam, a
teoria do potencial utilizada principalmente em pesquisas de fluidos perfeitos e pesquisas em
acústica, especialmente em acústica submarina. É importante constatar que essas duas áreas
apresentam dificuldades comuns, pois ambas desenvolvem análises em domínios infinitos, o
25
que para os padrões da época era incompatível com as técnicas numéricas de elementos
finitos e diferenças finitas.
Em várias áreas de pesquisa nos ramos das engenharias, existe a demanda de
conhecimentos sobre o comportamento de propagação de ondas em diferentes meios.
Fenômenos importantes de propagação de ondas ocorrem em domínio limitados e também em
domínios ilimitados. As análises de problemas realistas, que incluem domínios com
geometrias complexas e condições de contorno arbitrárias, não podem ser realizadas com uso
do ferramental analítico que fornece repostas matematicamente exatas. Estes problemas mais
complexos e realistas exigem a aplicação de ferramentas numéricas de aproximação de sua
solução. As ferramentas numéricas ganharam grande impulso na segunda metade do século
passado com o advento do computador digital.
Segundo Ciskowski, (1991) os problemas que envolvem domínios ilimitados foram
intensamente investigados no período de 1960 a 1975, tendo em vista duas grandes
aplicações, quais sejam, a teoria potencial aplicada a fluidos perfeitos e os problemas
acústicos lineares. O Método dos Elementos de Contorno teve sua origem em pesquisas de
problemas de acústica, nos trabalhos de Chen e Schweikert, (1963), Chertok, (1964), Shaw e
Fiedman, (1962), Banaugh e Goldsmith, (1963), e Shenck, (1968). Estes trabalhos tratavam de
problemas de radiação e espalhamento de ondas acústicas fazendo uso de Equações Integrais
de Contorno (EIC).
Embora o Método de Equações Integrais de Contorno tenha sido parte integrante no
desenvolvimento do Método dos Elementos de Contorno (MEC), como pode ser visto nos
artigos já citados, no período de 1960 até 1975 usava-se quase que exclusivamente a
nomenclatura Equação Integral (EI), sem distinguir a forma direta da indireta. Na verdade, a
primeira conferência de inovação desses métodos, organizado, no Rensselaer Polytecnhic
Institute em junho de 1975 por Tom Cruse e Frank Rizzo, (1975), sob o patrocínio da ASME,
foi denominada “Boundary Integral Equation Method: Computational in Applied Mechanics”.
Dois anos mais tarde foi publicado o texto de M.A. Jawson e G.T. Symm, (1977)
denominado “Integral Equation Methods in Potential Theory and Elastostatics”. Nenhum
desses trabalhos usava a frase “Métodos de Elementos de Contorno”. Nesse mesmo ano, outra
conferência tratando de métodos de inovações numéricas, foi realizada em Paris. Nesta
26
conferência uma pequena parte dos trabalhos, tratavam dos Métodos de Equações Integrais de
Contorno e a maior parte desses trabalhos, tratavam de análises numéricas baseadas no
Método dos Elementos Finitos e no Método das Diferenças Finitas, Ciskowski, (1991).
Passado um ano, o nome Método de Elementos de Contorno foi usado na segunda conferência
dedicada a estes métodos, organizado por C.A. Brebbia, (1978) em Southamptom, U.K. O
primeiro “texto introdutório”, como oposição à pesquisa monográfica, foi também publicado
por C.A. Brebbia em 1978. Desde aquele tempo o método tornou-se conhecido quase que
exclusivamente como Método de Elementos de Contorno (MEC).
A síntese das Equações Integrais de Contorno para problemas de elastodinâmica teve
o mérito de aproximar a comunidade que atuava em mecânica dos sólidos, daquela que
trabalhava com mecânica dos fluidos e problemas potenciais. A partir desse ponto houve um
rápido aumento nas aplicações em estruturas mecânicas e sólidos, baseados nestas
formulações de contorno, mas agora com o título “Métodos de Elementos de Contorno”.
As equações integrais de contorno (EICs), por sua vez foram utilizadas para solucionar
uma série de problemas de radiação e espalhamento de ondas acústicas. Alguns dos trabalhos
mais significativos estão citados a seguir. Chen e Schweikert, (1963) estudaram o problema
de radiação harmônica no tempo, seguido por Chertock, (1964) ainda tratando de problemas
de radiação do som. O problema de radiação e espalhamento de ondas acústicas com pequeno
comprimento de onda foi tratado por Yoon, (1990), estudado e implementado por Papini,
(1999). Também no trabalho de Debain, et al., (2004), é apresentada uma formulação para o
espalhamento de ondas.
3.3 – Sumário
Neste capítulo, foi visto que o Método de Elementos de Contorno destaca-se entre as
diferentes técnicas numéricas para solução numérica de problemas de campo aberto na área de
Acústica. Foram apresentadas de maneira sucinta as principais técnicas numéricas utilizadas
para solução de diversos fenômenos físicos em geral.
Uma breve história das Equações Integrais de Contorno aplicadas a problemas de
Acústica também foi apresentada. A fim de solucionar a equação de onda (2.18), apresentada
no capítulo 2, será desenvolvida, no próximo capítulo, a formulação bidimensional para o
problema de radiação acústica.
27
Capítulo 4
A Equação Integral de Contorno
_________________________
Este capítulo descreve a transformação da equação diferencial de Helmholtz,
deduzida no capítulo 2, em uma equação integral de contorno.
4.1– Introdução
Problemas de engenharia são freqüentemente descritos por leis físicas as quais são
comumente expressas por equações diferenciais parciais. Em muitos casos, uma representação
matemática alternativa e equivalente do problema é encontrada em termos de Equações
Integrais de Contorno. Com o avanço nas técnicas de modelagem numérica e o incremento na
capacidade de processamento dos computadores, métodos de modelagem baseados em
equações integrais de contorno podem ser agora usados para a simulação de muitos problemas
práticos de engenharia.
4.2– Equação Integral de Contorno
A equação de Helmholtz (2.20), equação diferencial clássica que descreve o problema
de radiação acústica de um cilindro infinito pulsante foi determinada no capítulo 2. Para que
seja encontrada a solução numérica dessa equação, faz-se necessário encontrar a Equação
Integral no Contorno para a mesma.
Seja B um corpo bidimensional imerso em um domínio infinito Ω, como representado
na figura 4.1 a seguir.
28
Figura 4.1 – Representação do domínio e dos contornos
A equação de Helmholtz (2.20) é válida em todos os pontos do domínio Ω. Para
transformá-la em uma equação integral, faz-se a suposição que ela seja satisfeita em um
conjunto limitado de pontos do domínio. Em alguns desses conjuntos, a equação pode ser não
nula, e para resolver essa questão, considera-se que a equação de Helmholtz não é mais nula
em todo o domínio, gerando assim um resíduo r. Com isso ela pode ser escrita da seguinte
forma:
r
r
r
∇2φ (r ) + k 2φ (r ) = r, ∀r ∈ Ω .
(4.1)
O resíduo da equação (4.1) é avaliado em cada ponto do domínio. Para solução da
mesma, utiliza-se o Método dos Resíduos Ponderados. Este método é baseado na inserção de
uma de ponderação u * , que faz com que a soma dos resíduos em todo o domínio torne-se
nula. Dessa maneira
∫r u
*
dΩ = 0 .
(4.2)
Ω
Substituindo-se a equação (4.1) na equação (4.2) obtém-se:
∫ (∇ φ + k φ ) u
2
2
*
dΩ = 0 .
(4.3)
Ω
Para obter-se a solução da equação acima, tem-se que fazer uso de algumas
identidades vetoriais, mostradas no anexo C.
29
Fazendo-se algumas manipulações algébricas e usando algumas propriedades do
gradiente verifica-se que
(
)
u* ∇2φ = ∇ ⋅ u* ∇ φ − ∇u* ⋅ ∇φ .
(4.4)
Retornando a equação (4.3) e aplicando a propriedade distributiva, obtém-se
∫ (u ∇ φ + uk φ ) dΩ = 0 .
*
2
2
(4.5)
Ω
Separando a equação (4.5) como uma soma de integrais, tem-se
∫ u ∇ φ dΩ + ∫ uk
*
2
Ω
2
φ dΩ = 0 .
(4.6)
Ω
Aplicando novamente as identidades vetoriais na equação (4.6) chega-se à seguinte equação
∫ [∇ ⋅ (u
*
]
)
∇ φ − ∇u* ⋅ ∇φ dΩ + ∫ u*k 2φ dΩ = 0 .
Ω
(4.7)
Ω
Aplicando algumas manipulações algébricas, temos
∫ ∇ ⋅ (u
*
)
(
)
∇ φ dΩ + ∫ u*k 2φ − ∇u* ⋅ ∇φ dΩ = 0 .
Ω
(4.8)
Ω
Aplicando o teorema da divergência na equação (4.8) obtém-se:
∫u
*
Γ
∂φ
dΓ + ∫ u * k 2φ − ∇u * ⋅ ∇φ dΩ = 0 .
∂n
Ω
(
)
(4.9)
Fazendo uso da identidade vetorial C.5 na equação (4.9), tem-se
∫u
Γ
*
∂φ
dΓ + ∫ u * k 2φ + φ ∇ 2 u * − ∇ ⋅ φ ∇u * dΩ = 0 .
∂n
Ω
[
(
)]
(4.10)
30
Separando a segunda integral em duas partes obtém-se
∫u
*
Γ
∂φ
dΓ + ∫ u * k 2 φ + φ ∇ 2 u * dΩ − ∫ ∇ ⋅ φ ∇ u * dΩ = 0 .
∂n
Ω
Ω
(
)
(
)
(4.11)
Aplicando-se o teorema da divergência na equação (4.11) encontra-se
*
∫u
Γ
∂φ
∂u *
dΓ − ∫ φ
dΓ + ∫ φ ∇ 2 u * + k 2 u * dΩ = 0 .
∂n
∂n
Γ
Ω
(
)
(4.12)
Utilizando-se a propriedade da função delta de Dirac (anexo D) na equação (4.12) obtém-se:
*
∫u
Γ
r
∂φ
∂u *
dΓ − ∫ φ
dΓ − φ (r ) = 0 .
∂n
∂n
Γ
(4.13)
r
A equação acima foi obtida para pontos de colocação r pertencentes ao domínio.
r
v
Considerando Y = r ' e X = r , a partir da equação (4.13) obtém-se
*
∫ u ( X,Y )
Γ
∂φ ( X )
∂u * ( X,Y )
dΓ ( X ) − ∫ φ ( X )
dΓ( X ) = φ (Y ) .
∂n( X )
∂n( X )
Γ
(4.14)
A equação acima mostra que é possível solucionar a equação integral (4.13) em
v
qualquer ponto Y = r ' interior ao domínio, utilizando apenas integrais ao longo do contorno.
Para isso, basta conhecer as funções φ e
∂φ
∂u *
no contorno, pois u* e
são funções
∂n
∂n
analíticas conhecidas (funções de Green). Neste trabalho, tem-se que
∂φ
= 1 e com isso,
∂n
precisa-se determinar somente φ .
A solução fundamental u* que representa a função de Green, para a Equação de
Helmholtz é apresentada para o corpo bidimensional, Ciskowski et al., (1991), como
31
u* =
i 1 r
H 0 kR ,
4
( )
(4.15)
e sua derivada por
r
∂u *
ik 1 r ∂R
,
q =
= − H 1 kR
∂n
4
∂n
*
( )
(4.16)
v
onde R é a distância entre o ponto X e o ponto de aplicação Y no domínio Ω.
Para resolver o problema, é necessário discretizar o contorno, dividindo em elementos
que podem ser lineares, quadráticos, cúbicos, etc.
4.3 – Equação Integral no Contorno
A equação (4.14) é válida para qualquer ponto Y pertencente ao domínio Ω. No Método
dos Elementos de Contorno, esta equação é aplicada no contorno. Com isso, faz-se necessária
a análise do que acontece com a equação, quando o ponto Y torna-se um ponto pertencente ao
r
contorno. Quando o ponto de aplicação Y coincidir com o ponto de colocação X, o valor de R
será zero, causando um problema de singularidade nas equações (4.15) e (4.16). Para resolver
esta singularidade, pode-se envolver o ponto Y por um pequeno círculo de raio ε, como
mostrado na figura 4.2, e examinar a solução no limite quando o raio ε → 0.
Figura 4.2 – Representação do domínio sobre um contorno 2D
32
O ponto Y é escolhido como o centro deste círculo, sendo que no limite, o raio ε → 0.
Dessa maneira, o ponto Y torna-se um ponto de contorno Γ , o que resulta numa expressão
que é uma particularização da equação (4.14). Como este processo limite depende apenas da
ordem da singularidade do potencial de velocidade
φ , que é a mesma para os operadores de
Helmholtz e Laplace, segundo Papini, (1999), é realizado um estudo das integrais de contorno
da equação (4.14) utilizando a solução fundamental para a equação de Laplace no domínio Ω.
A Equação de Laplace é dada por
∇ 2φ = 0 .
(4.17)
A solução fundamental para o caso isotrópico bidimensional, segundo Brebbia et al,
(1980), é
u* =
1 1
ln r  .
2π  R 
(4.18)
Suponha que o contorno é do tipo Γ2 e considere cada integral de contorno na equação
(4.14) decomposta em duas partes, sendo a primeira uma integral ao longo de Γ2−ε e a
segunda uma integral ao longo de Γε , isto é,
∫u
*
Γ
∂φ
∂φ
∂φ
dΓ = ∫ u *
dΓ + ∫ u *
dΓ ,
∂n
∂n
∂n
Γ2 − ε
Γε
(4.19)
∂u *
∂u *
∂u *
dΓ = ∫ φ
dΓ + ∫ φ
dΓ .
∂n
∂n
∂n
Γ2 − ε
Γε
(4.20)
e
∫
φ
Γ
Considerando a segunda parcela do segundo termo das equações (4.19) e (4.20) e substituindo
*
u * e ∂u pelo valor de suas soluções fundamentais, obtém-se
∂n
∫u
Γε
*
∂φ
1  1  ∂φ
dΓ = ∫
ln 
dΓ ,
∂n
2
 ε  ∂n
Γε π
(4.21)
33
e
−
∫
Γε
φ
1
∂u *
dΓ = − ∫ φ
dΓ .
2πε
∂n
Γε
(4.22)
Aplicando o limite nas equações (4.21) 3 (4.22) tem-se


 ∂φ
∂φ
lim  ∫ u *
dΓ  = lim  ∫
ε →0
∂n
Γε
 ε → 0 Γε ∂n

 ∂φ πε 2  1  
 1  1 
ln
d
Γ
=
lim
ln    = 0 ,


 ε →0 
 2π

 ε 


 ∂n 2π  ε  
(4.23)
e




∂u *
1
1
πε 

lim − ∫ φ
dΓ  = lim − ∫ φ
dΓ  = lim − φ
 = − φ (Y ) .
ε →0
ε
→
0
ε
→
0
2πε 
2

 Γε ∂n

 Γε 2πε

Note que como
(4.24)
ε vai pra zero no limite, o contorno Γ ε novamente pode ser
2−
representado como Γ2 e quando se integra sobre a singularidade obtém-se − 1 φ (Y ) . Com
2
isso, quando o ponto Y está sobre o contorno, aplicam-se as equações (4.21) e (4.22) na
equação (4.35), obtendo a seguinte expressão
*
∫u
Γ
∂φ
∂u *
1
dΓ − ∫ φ
dΓ = φ (Y ) .
∂n
∂n
2
Γ
(4.25)
Esta equação é conhecida como Equação Integral no Contorno. Uma maneira mais
geral de representá-la, na qual Y pode estar localizado no domínio, no contorno ou fora do
domínio, pode ser formulada utilizando-se um termo livre c(Y ) relacionado à posição de Y,
*
∫u
Γ
∂φ
∂u *
dΓ − ∫ φ
dΓ = c(Y )φ (Y ) .
∂n
∂n
Γ
(4.26)
34
4.4 – O Coeficiente c, termo livre do sinal de integração
O coeficiente c que aparece na equação (4.26), tem seu valor determinado quando é
conhecida a posição do ponto fixo Y em relação à (Ω ∩ Γ). Assim o valor determinado pela
relação
 θe
c = 1 −
 2π

 ,

(4.27)
depende somente do ângulo com origem no ponto Y, Menoni, (2004).
Se o ponto Y não pertencer a Ω ∩ Γ, sempre vai existir ε > 0, que determina um
círculo de raio ε e centro Y, tal que nenhum dos pontos deste círculo pertencem a Ω ∩ Γ, isto
é, a região angular é o ângulo externo θ e = 2π e que substituído na relação (4.27) determina
c = 0.
Por outro lado, se o ponto Y pertencer ao interior de Ω, sempre vai existir ε > 0, de
modo que o círculo de centro Y e raio ε esteja contido em Ω, determinando um ângulo interno
θ i = 2π , consequentemente θ e = 0 , da relação (4.27) resulta c = 1.
Consideremos agora, o ponto Y sobre o contorno Γ:
a) se é possível traçar uma reta tangente ao contorno Γ, passando pelo ponto Y, então o
ângulo externo ao domínio Ω é θ e = π , e portanto c =
1
. Nesse caso dizemos que o
2
ponto Y pertence à parte suave do contorno Γ.
b) O ponto Y está sobre a parte não suave do contorno Γ, isto é, o ponto Y está sobre um
canto no contorno. Nesse caso, o ângulo externo deve ser determinado.
Pode-se resumir essas possibilidades descritas acima da seguinte maneira:
35
0, se o ponto não pertence ao domínio ;
1
 , se o ponto Y, pertence à parte suave do contorno;
2
c=
1, se o ponto Y é um ponto interno ao domínio;
 θe
, se o ponto Y pertence a parte não suave do contorno.
 1 2π
4.5 – Tipos de Elementos de Contorno
O Método dos Elementos de Contorno pode ser formulado utilizando-se elementos
contínuos ou descontínuos para discretizar o contorno de um objeto. Os elementos
descontínuos têm nós que estão localizados entre as extremidades dos elementos. Logo, os
elementos não dividem nós com elementos adjacentes. Os elementos contínuos, por outro
lado, têm nós em suas extremidades e dividem estes nós com elementos adjacentes como na
Figura 4.3. Parte da razão pela preferência de usuários e pesquisadores quanto ao uso de
elementos contínuos justifica-se quando os dois tipos de elementos são comparados em
problemas que usam o mesmo tipo de malha. Tais comparações mostram, como na Figura 4.4,
que o modelo que utiliza elementos descontínuos contém 24 nós, já o modelo com elementos
contínuos tem apenas 12 nós. Como era de se esperar, o modelo com elementos descontínuos
requer mais tempo computacional no procedimento para a formação das matrizes. A partir da
Figura 4.5 uma resposta hipotética (como temperatura) é comparada com o tipo de resposta
que pode ser modelada ao utilizar-se os dois tipos de elemento. Esta figura justifica o uso do
nome descontínuo. Existem saltos entre as respostas dadas pelos nós das extremidades desses
elementos. Os elementos descontínuos são indicados para tratar hipersingularidades, segundo
Papini, (1999).
Desta forma, utilizando-se elementos contínuos o contorno do objeto pode ser
discretizado por elementos isoparamétricos constantes, lineares ou de ordem superior, como
os cúbicos e quadráticos, segundo a Figura 4.6. A definição de elementos isoparamétricos é
dada para aqueles elementos cuja geometria é interpolada usando-se as mesmas funções que
são usadas para interpolar as respostas destes elementos.
36
No caso particular de elementos constantes os valores de φ e
∂φ
são considerados
∂n
constantes ao longo de cada elemento e são iguais aos valores do ponto nodal, o qual tem
posição convencionada no centro do elemento. Os pontos nas extremidades do elemento são
usados apenas para definir a geometria do contorno de um objeto, que está inserido no
domínio de interesse. Assim, para este tipo de elemento o contorno é sempre plano (ou liso),
conseqüentemente o coeficiente de c(Y) é sempre 1 2 no contorno, segundo Ciskowski et al.,
(1991); Brebbia et al., (1998).
Figura 4.3 – Malhas compostas de diferentes tipos de elementos, Papini (1999)
37
Figura 4.4 – Malhas com o mesmo número de elementos, Papini (1999)
Figura 4.5 – Aproximação realizada por elementos contínuos e descontínuos, Papini (1999)
38
Figura 4.6 (a) – Elementos de contorno: elementos constantes, Ciskowski et al., (1991)
Figura 4.6 (b)– Elementos de contorno: elementos lineares, Ciskowski et al., (1991)
Figura 4.6 (c) – Elementos de contorno: elementos quadráticos, Ciskowski et al., (1991)
39
4.6– Discretização e Colocação
A discretização do contorno Γ consiste na divisão do mesmo, em n partes. O
segmento de reta que une os extremos consecutivos de cada uma destas partes, é denominado
de elemento de contorno. Desse modo, se Γ é o contorno de um domínio Ω ⊂ ℜ 2 , e Γ é
discretizado em n elementos, temos que:
n
Γ = UΓj ,
(4.28)
j =1
em que cada Γ j , é um elemento de contorno Γ , conforme pode ser visto na figura 4.7.
Figura 4.7 – Contorno discretizado
A avaliação do contorno por meio da Equação Integral de Contorno, segundo a
equação diferencial que governa o problema no meio, é determinada pela contribuição que
cada elemento exerce sobre esse contorno. Assim, faz-se necessário a caracterização de cada
um dos elementos no contorno discretizado.
Seja Γ j um elemento de Γ , (x j , y j ) e (x j +1 , y j +1 ) as coordenadas de seus extremos e L
se comprimento, definido por:
40
L=
(x
− x j ) + ( y j +1 − y j ) .
2
j +1
2
(4.29)
v
O vetor normal n , para o elemento Γ j , pode ser calculado considerando-se um vetor
v
v 2 , normalizado e unitário, construído utilizando-se as coordenadas das extremidades do
v
v
elemento, e outro vetor v1 = (0,0,1) ortogonal a v 2 . Assim, o vetor normal ao elemento é
v
v
calculado pelo produto vetorial entre v1 e v 2 .
v v v
n = v1 × v 2 =
i
(x j +1 − x j )
L
0
j
k
(y j +1 − y j )
0,
L
0
1
x j +1 − x j

v  y j +1 − y j
n = 
i ,−
j ,0k  .
L
L


(4.30)
(4.31)
Caracterizado cada um dos elementos Γ j , considere x j = ( x, y ) , um nó do elemento
Γ j e y i = ( xi , y i ) um ponto de Ω ∪ Γ . É comum na literatura denominar xi como ponto de
v
colocação e y j como ponto de campo. A distância entre os pontos xi e y j , denotada por R é
definida como sendo a distância Euclidiana dada por:
v
R = y j − xi = ( x − xi )i + ( y − y i ) j ,
(4.32)
v
R ( x, y ) =
(4.33)
e
( x − x i )2 + ( y − y i )2 .
Para a discretização das variáveis físicas e geométricas do problema, assume-se uma
distribuição constante das variáveis u * e
∂u *
ao longo dos elementos em que o contorno foi
∂n
discretizado.
41
Assim, da equação (4.25) pode-se escrever
rr
n
n ∂φ
r
r
r ∂u * (r ,r ')
1
j
φ j (r ') + ∑ φ j (r ') ∫
dΓ(r ) = ∑
2
∂n
j =1
j =1 ∂n
Γj
As integrais
∂u *
∫ ∂n dΓ
e
Γj
∫ u dΓ
*
rr
r
∫ u (r ,r ') dΓ(r ) .
*
(4.34)
Γj
são chamadas de coeficientes de influência, pois relacionam a
Γj
influência da solução no ponto P, quando a solução fundamental é integrada sobre o elemento
Q. Renomeando-se as integrais, tem-se
Gij = ∫ u * dΓ ,
(4.35)
Γj
e
)
∂u *
H ij = ∫
dΓ ,
∂n
Γj
(4.36)
onde i representa o ponto de colocação e j o elemento em consideração a ser integrado. Assim
a equação (4.34) pode ser escrita da seguinte maneira
n
n ∂φ
)
1
v
j
φ j (r ') + ∑ φ j H ij = ∑
Gij .
2
∂
n
j =1
j =1
(4.37)
Denominando Hij por
)
 H ij

H ij =  )
1
 H ij + 2
para i ≠ j
para i = j
(4.38)
r
pois φ (r ') é zero para os elementos que não contém a singularidade e um para os que contém.
Assim, chega-se a seguinte equação
42
n
∑φ
n
j
H ij = ∑ Gij q j ,
j =1
onde
(4.39)
j =1
∂φ
= 1 (condição de Neumann) conforme definido no capítulo 2.
∂n
Uma vez conhecida a Equação Integral de Contorno para o operador de Helmholtz em
sua forma discretizada, deve-se descrever sobre cada elemento a variação da geometria e das
variáveis envolvidas no problema. Estas variações podem ser constantes, lineares,
quadráticas, cúbicas ou de ordem superior. A forma mais conveniente de descrever o
comportamento de um elemento Γ j é através da utilização das chamadas funções de forma,
que utilizam um número de pontos (nós) sobre cada elemento onde o valor da variável é
conhecido. Assim, uma função de forma quadrática exigirá a presença de três nós sobre cada
elemento, ou seja, um no ponto médio e os outros nos extremos do elemento, Ciskowski
(1991).
A geometria do elemento Γ j , pode ser descrita pelas coordenadas dos nós e pela
função de forma correspondente, isto é:
m
x(ξ ) = ∑ N s (ξ )x s ,
(4.40)
s =1
e
m
y (ξ ) = ∑ N s (ξ ) y s ,
(4.41)
s =1
em que ξ corresponde a um novo sistema de coordenadas que é local, para aquele elemento
considerado, N s são as funções de forma satisfazendo as seguintes condições:
s
1, no nó
N s (ξ ) = 
,
0, nos demais nós
(4.42)
43
e m corresponde ao número de nós sobre cada elemento e está relacionado com o tipo de
função de forma a ser utilizada, Menoni, (2004).
Considere um elemento Γ j , para o qual se deseja aproximar sua geometria, a partir de
uma função de forma linear. Para isto define-se um novo sistema de coordenadas que é local a
Γ j , utilizando a variável ξ com sua origem no ponto médio do elemento e assumindo valores
no intervalo [-1, 1] (figura 4.8 e 4.9).
Figura 4.8 – Elemento Linear
Figura 4.9 – Sistema de coordenadas no elemento, Ciskowski et al., (1991).
A geometria de Γ j pode ser aproximada pelas coordenadas dos dois nós e pela função
de forma, como segue:
2
x(ξ ) = ∑ N i (ξ )xi = N 1 (ξ )x1 + N 2 (ξ )x 2 ,
i =1
44
e
(4.43)
2
y (ξ ) = ∑ N i (ξ ) y i = N 1 (ξ ) y1 + N 2 (ξ ) y 2 ,
i =1
em que neste caso, N i (ξ ) assume a seguinte forma em cada nó:
N1 (ξ ) = a1ξ + a 2 ,
e
(4.44)
N 2 (ξ ) = a3ξ + a 4 .
Os coeficientes a1 , a 2 , a3 , a 4 , são constantes que podem ser determinadas utilizando as
condições de contorno, como segue:
N 1 (− 1) = 1, N 1 (1) = 0 ,
e
(4.45)
N 2 (− 1) = 0, N 2 (1) = 1 ,
e substituindo as equações (4.45) nas equações (4.44), tem-se:
N1 (ξ ) =
1
(1 − ξ ) ,
2
e
(4.46)
N 2 (ξ ) =
1
(ξ + 1) .
2
Este procedimento é aplicável, para qualquer que seja a ordem da função de forma
utilizada. Estas funções de forma podem ser usadas para aproximar as variáveis do problema.
r ∂φ
Assim, se φ representa o potencial e q =
o fluxo, elas podem ser descritas como:
∂n
m
φ (ξ ) = ∑ N s (ξ )φ s ,
s =1
e
(4.47)
45
q (ξ ) =
∂u (ξ ) m
 ∂φ 
= ∑ N s (ξ )  .
∂n
 ∂n  s
s =1
Para o comprimento do elemento Γ j dado pela equação (4.29), pode-se considerar o
seu infinitésimo dΓs , dado como:
2
dΓs =
(dx )
2
+ (dy )
2
2
 dx   dy 
=   +   dξ ,
 dξ   dξ 
(4.48)
em que:
m
dN s
dx
,
= ∑ xs
dξ s =1
dξ
e
(4.49)
m
dN s
dy
.
= ∑ ys
dξ s =1
dξ
Fazendo-se o jacobiano da transformação
2
2
 dx   dy 
J (ξ ) =   +   ,
 dξ   dξ 
(4.50)
dΓs = J (ξ )dξ .
(4.51)
tem-se
v
4.7 – Derivada da solução fundamental em relação a normal n
Para calcular o termo
∂ u*
, é necessário determinar o valor da derivada direcional de
∂n
v
R , pois:
46
∂ u * ∂u * ∂ R
.
=
∂ n ∂R ∂ n
(4.52)
Pela definição de derivada direcional, tem-se:
r
∂R
=∇R⋅ n ,
∂n
(4.53)
onde,
r  y − y1
x −x

n = 2
i,− 2 1 j , 0k  ,
L
 L

(4.54)
 ∂R ∂R 
 ,
∇R = 
,
∂ x ∂ y
(4.55)
e
v
representam o vetor unitário normal ao elemento e o gradiente de R , respectivamente. Assim,
tem-se que:
∂ R x-xc
,
=
∂x R
(4.56)
∂ R y-yc
.
=
∂y R
(4.57)
e
v
Realizando o produto escalar entre o gradiente de R e o vetor unitário normal, tem-se:
 ( x − xc )( y 2 − y1 ) + ( y − yc )(x1 − x2 ) 
=
.
RL


r
Derivando a solução fundamental u * em relação a n ,tem-se:
r
∂R
=∇R⋅ n
∂n
(4.58)
47
∂u *
∂R
ik
= − H 11 (kR )
,
∂n
4
∂n
(4.59)
ou
∂u *
ik
 ( x − xc )( y 2 − y1 ) + ( y − yc )( x1 − x2 ) 
= − H 11 (kR ) 
.
∂n
4
RL


(4.60)
4.8 – Cálculo dos Elementos não Singulares das Matrizes G e H
Sendo X1 e X2 as coordenadas das extremidades do elemento, tem-se a partir das
equações (4.15) e (4.35), que:
X2
Gij =
i
∫ 4 H (kR )dΓ .
1
0
(4.61)
X1
Integrando essa equação pela regra da quadratura gaussiana (desenvolvida no
Apêndice B), tem-se:
nptg
i 1
L
i
L
H 0 (kR(ξ )) dξ = ∑ H 01 (kR(ξ m ))wm .
−1 4
2
2
m =1 4
Gij = ∫
+1
(4.62)
Expandindo a função de Hankel utilizando-se as funções de Bessel (mostradas no
Apêndice A), pode-se rescrever a equação (4.62) como:
nptg
L 1
L
i
Gij = ∑  J 0 (kR(ξ m ))wm − Yo (kR(ξ m ))wm  .
2 4
2
m =1  4
(4.63)
Para calcular os elementos da matriz H, tem-se a partir das equações (4.16) e (4.36)
que:
48
X2
H ij = ∫ −
X1
ik 1
∂R
H 1 (kR ) r dΓ .
4
∂n
(4.64)
Integrando essa equação pela regra da quadratura, tem-se:
X2
H ij = ∫ −
X1
nptg
ik 1
∂R L
ik
∂R
L
H 1 (kR(ξ )) r dξ = ∑ − H 11 (kR(ξ )) r wm .
4
∂n 2
4
∂n
2
m =1
(4.65)
Expandindo a função de Hankel a partir das funções de Bessel, tem-se:
nptg
H ij = ∑ −
m =1
ik
∂R
L k
∂R
L
J 1 (kR(ξ )) r wm + Y1 (kR(ξ )) r wm .
4
∂n
2 4
∂n
2
(4.66)
4.9 – Cálculo dos Elementos Singulares das Matrizes G e H
Os elementos da diagonal das matrizes G e H referem-se aos elementos singulares, ou
seja, é o caso onde i = j. Assim, a presença da singularidade nesses elementos, devido à
r
solução fundamental quando R → 0 , requer uma integração com fórmulas especiais para se
obter uma maior acuidade. No caso dos elementos de contorno constantes as fórmulas são
calculadas analiticamente.
Para o calculo do elemento H ii , sabe-se que o vetor unitário normal à superfície de
contorno e a coordenada desse elemento são perpendiculares, e que o raio de integração está
sobre o próprio elemento. Desta forma, tem-se:
r
∂R
=∇R ⋅ n =0 .
∂n
(4.67)
)
Automaticamente, o integrando de H ii , de acordo com as equações (4.36) e (4.52),
serão igual a zero, consequentemente H ii =
1
. Os elementos Gii da diagonal requerem
2
manuseio especial. Assim, expandindo a equação de Hankel presente na equação (4.61), temse:
49
Gij = ∫
X2
X1
X
X
i 1
1 2
i 2
H 0 (kR )dΓ = − ∫ Y0 (kR )dΓ + ∫ J 0 (kR )dΓ .
4
4 X1
4 X1
(4.68)
A singularidade apresentada na equação (4.68) foi resolvida, neste trabalho, com o
aumento do número de pontos de Gauss, para 10 pontos, na integração feita utilizando-se a
quadratura gaussiana.
Assim, retornando a equação (4.39), a solução fundamental é aplicada em cada nó, o
que possibilita verificar a influência de todos os outros elementos no nó da singularidade e
dele nele mesmo, obtendo-se um sistema de equações expresso na forma matricial, para cada
ponto do contorno, como segue
r
r
Hφ = Gq ,
(4.69)
r
r
onde H e G são duas matrizes N x N, e φ e q são vetores de comprimento N.
Inserindo-se todos os elementos das matrizes H e G correspondentes às condições de
contorno desconhecidas do lado esquerdo e aqueles correspondentes às condições de contorno
conhecidas do lado direito, e multiplicando as matrizes do lado direito, forma-se o seguinte
sistema de equações:
r r
Ay = b ,
(4.70)
r r
r
r
onde y é o vetor dos valores de contorno desconhecidos de φ e q . O vetor b é encontrado
r
r
multiplicando as colunas correspondentes de H ou de G pelos valores conhecidos de φ e q .
É interessante notar que as variáveis podem ser uma mistura dos potenciais e suas derivadas.
Após a resolução no contorno é possível calcular qualquer valor interno do potencial ou de
sua derivada.
50
4.10 – Sumário
Neste capítulo, foi desenvolvida a formulação direta do MEC para a Equação de
Helmholtz (2.20). Após o desenvolvimento das integrais para obtenção dos elementos das
matrizes G e H, seguiu-se a implementação computacional, utilizando-se o software Matlab.
Foi verificada a existência de singularidades em algumas integrais e apresentado o método
para solução das mesmas. O próximo capítulo apresenta os resultados obtidos e os compara
com os resultados obtidos por Papini, (1999).
51
Capítulo 5
Resultados Numéricos e
Comparações
_________________________
Neste capítulo são apresentados os resultados analíticos e numéricos encontrados
durante este trabalho sendo comparados com os resultados encontrados por Papini, (1999) e
Yoon, (1990).
5.1 – Resultado Analítico e Comparações
A solução analítica da equação diferencial de Helmholtz (eq 2.20) foi implementada
no software MatLab e os resultados obtidos foram comparados com os obtidos por Papini,
(1999), mostrados na Tabela 5.1 e exibidos na figura 5.1.
Na análise dos dados, foram calculados o erro absoluto, e e o erro relativo, er
calculados da seguinte forma
e = val1 − val ,
(5.1)
e
er =
e
,
val1
(5.2)
52
onde e é o erro absoluto, val1 é o valor do potencial de velocidade obtido por Papini, (1999),
val é o valor do potencial de velocidade obtido neste trabalho e er é o erro relativo ao
resultado obtido por Papini, (1999).
Freqüência
Analítico
Analítico
(Hz)
Papini
(Implementada)
2
(m /s)
2
(m /s)
62,5
0,770946776
125
Erro absoluto
Erro relativo
0,770946769
-7,000000E-09
-9,079745E-09
0,419396006
0,419396013
7,000000E-09
1,669067E-08
250
0,215903228
0,215903243
1,500000E-08
6,947557E-08
500
0,108860012
0,108860013
1,000000E-09
9,186110E-09
Tabela 5.1 – Comparação de resultados Implementados
Analítica Papini
Analítica Implementada
0,9
2
Módulo do Potencial de Velocidade (m /s)
0,8
0,7
0,6
0,5
0,4
0,3
0,2
0,1
0
62,5
125
Freqüência (Hz)
250
500
Figura 5.1 – Comparação de resultados implementados
Os resultados apresentados na Tabela 5.1 e mostrados na figura 5.1 são bastante
satisfatórios quando comparados aos resultados analíticos. Para facilitar a visualização os
pontos foram unidos por linhas contínuas.
53
5.2 – Resultado Numérico e Comparações
Para a solução numérica do potencial de velocidade radiado (φ ) , cuja condição de
contorno de Neumann ao longo de todo o seu contorno é
∂φ
= 1 , foram utilizados 6 pontos de
∂n
Gauss para a integração nos elementos que não singulares e 10 pontos de Gauss para os
elementos singulares. Os resultados são comparados com os da literatura e apresentados por
meio de figuras e tabelas.
Na figura 5.2 a seguir está representado o cilindro (bidimensional) discretizado em 12
elementos. Os pontos, que representam o início e o fim dos elementos, estão indicados pelos
asteriscos. A linha, localizada no centro do elemento, representa a reta normal a cada
elemento.
Figura 5.2 – Ilustração bidimensional da discretização do cilindro em 12 elementos
As discretizações no contorno do cilindro foram consideradas para 8, 16, 32, 64, 128 e
256 elementos lineares (respectivamente 8, 16, 32, 64, 128, 256 nós), testando as freqüências
centrais das bandas de uma oitava (de 62,5 Hz até 500 Hz) e analisando os valores
54
encontrados para o módulo do potencial de velocidade radiado. A seguir são apresentados os
resultados encontrados e suas respectivas representações gráficas.
8 Elementos
Freqüência
Analítica
Numérica
(Hz)
(m2/s)
(m2/s)
Erro absoluto Erro relativo
62,5
0,770946776 0,789957763
0,019010987
0,024659273
125
0,419396006 0,439270497
0,019874491
0,047388365
250
0,215903228 0,233108076
0,017204848
0,079687776
500
0,108860012 0,112451904
0,003591892
0,032995514
Tabela 5.2 – Resultados para discretização com 8 elementos
Analítica
Numérica
0,9
2
Módulo do Potencial de Velocidade (m /s)
0,8
0,7
0,6
0,5
0,4
0,3
0,2
0,1
0
62,5
125
Freqüência (Hz)
250
500
Figura 5.3 – Representação gráfica dos resultados obtidos com a discretização com 8 elementos
Na implementação com 8 elementos verifica-se que a solução encontrada é
satisfatória. Mesmo com um baixo número de elementos, observa-se o erro absoluto é de
muito pequeno.
55
16 Elementos
Freqüência
Analítica
Numérica
(Hz)
2
(m2/s)
(m /s)
Erro absoluto
Erro relativo
62,5
0,770946776 0,778077803
0,007131027
0,009249701
125
0,419396006 0,441867827
0,022471821
0,053581390
250
0,215903228 0,221561809
0,005658581
0,026208876
500
0,108860012 0,101519375
-0,007340637
-0,067431896
Tabela 5.3 – Resultados para discretização com 16 elementos
Analítica
Numérica
0,9
2
Módulo do Potencial de Velocidade (m /s)
0,8
0,7
0,6
0,5
0,4
0,3
0,2
0,1
0
62,5
125
Freqüência (Hz)
250
500
Figura 5.4 – Representação gráfica dos resultados obtidos com a discretização com 16 elementos
Na implementação com 16 elementos, verifica-se uma melhor aproximação para as
freqüências acima de 250 Hz. Mas no geral, a solução ainda pode ser melhorada.
56
32 Elementos
Freqüência
Analítica
Numérica
(Hz)
2
(m2/s)
(m /s)
Erro absoluto
Erro relativo
62,5
0,770946776 0,774036535
0,003089759
0,004007746
125
0,419396006 0,437708957
0,018312951
0,043665058
250
0,215903228 0,218662812
0,002759584
0,012781578
500
0,108860012 0,104656653
-0,004203359
-0,038612516
Tabela 5.4 – Resultados para discretização com 32 elementos
Analítica
Numérica
0,9
2
Módulo do Potencial de Velocidade (m /s)
0,8
0,7
0,6
0,5
0,4
0,3
0,2
0,1
0
62,5
125
Freqüência (Hz)
250
500
Figura 5.5 – Representação gráfica dos resultados obtidos com a discretização com 32 elementos
Com a implementação para 32 elementos, observa-se que a solução para freqüências
abaixo de 125 Hz, e acima de 250 Hz aproximam-se mais rapidamente da solução procurada.
57
64 Elementos
Freqüência
Analítica
Numérica
(Hz)
2
(m2/s)
(m /s)
Erro absoluto
Erro relativo
62,5
0,770946776 0,772387266
0,001440490
0,001868469
125
0,419396006 0,430581709
0,011185703
0,026670981
250
0,215903228 0,217315829
0,001412601
0,006542751
500
0,108860012 0,106478591
-0,002381421
-0,021875994
Tabela 5.5 – Resultados para discretização com 64 elementos
Analítica
Numérica
0,9
2
Módulo do Potencial de Velocidade (m /s)
0,8
0,7
0,6
0,5
0,4
0,3
0,2
0,1
0
62,5
125
Freqüência (Hz)
250
500
Figura 5.6 – Representação gráfica dos resultados obtidos com a discretização com 64 elementos
Com 64 elementos, a solução implementada se aproxima, ainda mais da solução da
literatura, mostrando que quanto maior o número de elementos, maior a precisão da solução
encontrada.
58
128 Elementos
Freqüência
Analítica
Numérica
(Hz)
2
(m2/s)
(m /s)
Erro absoluto
Erro relativo
62,5
0,770946776 0,771642618
0,000695842
0,000902581
125
0,419396006 0,425471932
0,006075926
0,014487324
250
0,215903228 0,216622257
0,000719029
0,003330330
500
0,108860012 0,107585143 -0,001274869
-0,011711086
Tabela 5.6 – Resultados para discretização com 128 elementos
Analítica
Numérica
0,9
2
Módulo do Potencial de Velocidade (m /s)
0,8
0,7
0,6
0,5
0,4
0,3
0,2
0,1
0
62,5
125
Freqüência (Hz)
250
500
Figura 5.7 – Representação gráfica dos resultados obtidos com a discretização com 128 elementos
Com a implementação utilizando-se 128 elementos, tem-se uma boa aproximação da
solução numérica. Faz-se ainda uma última implementação com 256 elementos, mostradas a
seguir.
59
256 Elementos
Freqüência
Analítica
Numérica
(Hz)
2
(m2/s)
(m /s)
Erro absoluto
Erro relativo
62,5
0,770946776 0,771288794
0,000342018
0,000443634
125
0,419396006 0,422547304
0,003151298
0,007513896
250
0,215903228 0,216266385
0,000363157
0,001682036
500
0,108860012 0,108199459
-0,000660553
-0,006067912
Tabela 5.7 – Resultados para discretização com 256 elementos
Analítica
Numérica
0,9
2
Módulo do Potencial de Velocidade (m /s)
0,8
0,7
0,6
0,5
0,4
0,3
0,2
0,1
0
62,5
125
Freqüência (Hz)
250
500
Figura 5.8 – Representação gráfica dos resultados obtidos com a discretização com 256 elementos
Verifica-se que para a implementação com 256 elementos, a aproximação da solução
numérica está bem próxima da solução de referência.
O erro está na terceira casa decimal, mostrando que essa é uma boa aproximação. Um
melhor resultado deverá ser obtido se for feito o tratamento das singularidades apresentadas.
60
Pode-se concluir, pelas análises feitas anteriormente, que ao aumentar o número de
elementos, o valor do módulo do potencial de velocidade aproxima-se do valor obtido por
Papini, (1999) e Yoon, (1990).
5.3 – Sumário
Neste capítulo foi solucionado o problema da radiação acústica de um cilindro infinito
pulsante. Aplicou-se para a solução numérica, o Método dos Elementos de Contorno. Os
resultados obtidos a partir deste método mostraram-se me bem próximos dos resultados
previamente publicados na literatura. Os erros encontrados devem-se ao fato do não
tratamento adequado das singularidades encontradas no cálculo da matriz G. O artifício
utilizado para contornar o problema da singularidade, foi aumentar o número de pontos de
Gauss, para 10 pontos na realização da integração Gaussiana, nos elementos da diagonal
principal da matriz G.
61
Capítulo 6
Conclusão Geral e Proposta para
Trabalhos Futuros
________________________
6.1 − Conclusão Geral
O presente trabalho dedicou-se a formulação analítica do problema de radiação
acústica de um cilindro infinito pulsante, onde foi obtida a solução para o caso em que a
superfície do cilindro está vibrando uniformemente na direção radial com uma amplitude de
Vo metros por segundo, numa freqüência de f Hertz. Esta solução analítica foi implementada
e comparada com os valores exatos, observando-se que os resultados foram bastante
satisfatórios.
Foi realizado o estudo da formulação numérica do mesmo problema, utilizando-se o
Método dos Elementos de Contorno (MEC), a partir da forma direta. A discretização do
cilindro foi escolhida como elementos de contorno constantes, para a solução da Equação de
Helmholtz.
Nessa formulação do MEC, foi obtida a equação integral de contorno para a equação
que descreve o problema da radiação acústica de um cilindro infinito pulsante. O
desenvolvimento matemático das integrais de contorno, a partir do MEC, às vezes, não é tão
fácil. As dificuldades estão relacionadas com as soluções fundamentais requeridas nessa
formulação que envolvem as Funções de Hankel e Bessel. Além disso, não se pode esquecer
dos problemas de singularidades de algumas integrais e que a matriz do sistema final é cheia e
não simétrica.
62
A solução numérica, utilizando-se o método direto do MEC, foi implementada e foram
realizados testes, alterando-se o número de elementos na discretização do cilindro. Verificouse que mesmo para um número pequeno de elementos, a solução é bem próxima da solução
exata. Isso mostra a robustez do Método de Elementos de Contorno para solução deste
problema.
A aproximação dos resultados apresentou-se satisfatória, mesmo sem o tratamento das
singularidades apresentadas por algumas integrais. O devido tratamento destas singularidades
deve melhorar ainda mais a precisão da solução encontrada, mesmo com poucos elementos
6.2 − Propostas para trabalhos futuros
O presente trabalho apresenta uma pequena aplicação do Método de Elementos de
Contorno a problemas de Acústica. Algumas sugestões de continuidade da análise são
formuladas na seqüência:
a) Implementação de elementos polinomiais de ordem mais alta, por exemplo,
quadrática e cúbica, visando analisar ganhos de precisão e eficiência computacional
associados a cada tipo de discretização.
b) Investigação das singularidades apresentadas pelos elementos da diagonal da matriz
G e analisar melhores soluções que melhor aproximem os resultados encontrados.
c) Analisar o caso de espalhamento de ondas por um cilindro, utilizando-se o Método
de Elementos de Contorno.
d) Estudar o comportamento de barreiras acústicas, utilizando o Método de Elementos
de Contorno.
63
Referências Bibliográficas
________________________
Afonso, M.M., Métodos Híbridos na solução de problemas de espalhamento eletromagnético.
Tese de Doutorado, UFMG, Belo Horizonte, 2003.
Bistafa, S.R., Acústica aplicada ao controle do Ruído, Ed. Edgard Blücher, 2006.
Boyce, W.E. e Diprima, R.C., Equações diferenciais elementares e problemas de valor de
contorno, Ed. Guanabara Koogan, 1990, Cap.6, pp 236 – 239.
Brebbia, C. A., Dominguez J., Boundary Elements: An Introductory Course. WitPress,
Boston, 1998.
Brebbia, C. A., Walker, S., Boundary Element Techniques in Engineering. London, 1980.
Chapman, S.J., Programação em Matlab para engenheiros, Editora Thompson Learning, 2006.
Chen, L.H., and Schweikert, D.G, Sound radiation from an arbitrary body, The Journal of the
Acoustical Society of America, Volume 35, pp. 1626-1632, 1963.
Chertock, George, Sound Radiation from Vibrating Surfaces, The Journal of the Acoustical
Society of America, Volume 36, Número 7, Julho 1964, Páginas 1305 – 1313.
Ciskowski, R.D., Brebbia, C.A. Boundary Element Methods in Acoustics. SouthamptonBoston: Computational Mechanics Publications, 1991.
Debain, E.P., Trevelyan, J., Bettess, P., Wave Boudary Elements: a theoretical overview
presenting applications in scattering of short waves. Engineering Analysis with Boundary
Elements, Volume 28, Issue 2, February 2004, Pages 131-141
64
Greco, M., Análise do problema harmônico de radiação e difusão acústica, usando o Método
dos Elementos de Contorno. Dissertação de mestrado, USP, São Carlos, 2000.
Juhl, P.M., The Boundary Element Method for sound field calculations. Tese de Doutorado,
Technical University of Denmark, Lyngby, 1993.
Kinsler, L.E., Frey, A.R., Coppens, A.B., Sanders, J.V. Fundamentals of Acoustics. New
York: John Wiley & Sons, 1982.
Menoni J.A., Formulação e implementação da versão direta do Método dos Elementos de
Contorno para tratamento de problemas acústicos estacionários bidimensionais diretos e
inversos. Tese de Doutorado, Unicamp, Campinas, 2004.
Nepomuceno, L.X., Acústica. Editora Edigard Blücher Ltda, 1977.
Papini, G. S., Estudo Numérico de Barreiras Acústicas. Dissertação de Mestrado, UFMG,
Belo Horizonte, 1999.
Schenk, Harry A., Improved Integral Formulation for acoustic radiation problems, J. Acous.
Soc. Amer., vol. 44, 1968, pp 41 – 58.
Yoon, W.S., Park, J.M., Eversman,W. Two- Dimensional radiation and scattering at short
wave lenght, Journal of Vibration and Acoustics, V. 112, 1990.
Ziomek, L.J. Fundamentals of Acoustic Field Theory and Space-Time Signal Processing.
Boca Raton: CRC Press, 1995.
65
Anexo A
________________________
Funções de Bessel utilizadas
Funções de Hankel de 1ª classe
A solução fundamental apresentada para o caso bidimensional e sua derivada em
relação à reta normal são constituídas de funções de Hankel de 1ª classe de ordem 0 e 1,
respectivamente.
H 0(1) = J 0 ( x) + i ⋅ Y0 ( x)
,
(a)
e
(A.1)
H1(1) = J1 ( x) + i ⋅ Y1 ( x ) .
(b)
J n ( x) é uma função de Bessel de primeiro tipo, representada por:
2⋅ k + n
 x
(−1) ⋅  
∞
2
J n ( x) = ∑
,
k = 0 k!⋅Γ (k + n + 1)
k
(A.2)
onde:
n = 0,1,2,...
Função Gama: Γ(m + 1) = m!
66
Yn (x) é uma função de Bessel do segundo tipo representada por:
Yn ( x) =

2   x
1 n −1
x
⋅ ln  + γ  ⋅ J n ( x) − ⋅ ∑ (n − k − 1)!⋅ 
π  2 
π k =0
2
2⋅ k + n
2⋅ k + n
x
 
∞
1
2
k
− ⋅ ∑ (− 1) ⋅ [Φ(k ) + Φ (n + k )] ⋅  
π k =0
k!⋅(n + k )!
(A.3)
onde:
n = 0, 1, 2, ...
Constante de Euler: γ ≅ 05772156649015329 , Greco, (2000)
p
1
Φ( p ) = ∑  , Φ (0 ) = 0 .
m=1  m 
Das equações (A.2) e (A.3), são retiradas as seguintes funções de Bessel:
(− 1)m (x )2m
2m
(m!)2
m=1 2
∞
J 0 ( x) = 1 + ∑
(A.4)
e
Y0 ( x) =
∞
(− 1)m+1 hm ( x) 2m 
2 
x
γ
+
ln
J
(
x
)
+
 0


∑
2m
π 
2
(m!)2
m=1 2

(A.5)
onde:
m
1
hm = ∑ .
h =1 h
O comportamento assintótico das funções de Bessel é caracterizado conforme mostra
as equações a seguir (Boyce, 1990):
Lim J 0 ( x ) = 1
x →0
(a)
67
e
(A.6)
Lim Y0 ( x) =
x →0
2
π
ln( x)
(b)
E por outro lado se x cresce ilimitadamente,
π

cos x − 
4

Lim J 0 ( x) =
x →∞
πx
(a)
e
(A.7)
π

sen x − 
4

Lim Y0 ( x) =
x →∞
πx
(b)
68
Anexo B
________________________
Integração Numérica
Quadratura Gaussiana
Em problemas envolvendo o método de elementos de contorno, a integração
Gaussiana é encontrada com maior freqüência, dentre as técnicas de integração conhecidas na
literatura. Isso se deve ao fato talvez, de que esta técnica, avalia a função interpolante de
modo otimizado, ou seja: a partir da escolha dos pontos x1, x2, x3, ..., xnpg, no intervalo de
integração, seus respectivos pesos w1, w2, w3, ..., wnpg, (npg: número de pontos de Gauss) são
determinados, de forma a minimizar o erro durante a aproximação da integral Menoni (2004).
Assim temos:
b
∫
a
npg
f ( x)dx ≈ ∑ wi f (xi )
(B.1)
i
em que f ( x ) é uma função qualquer , wi são conhecidos como os pesos e npg o número de
pontos de Gauss. Os pesos wi são tabelados e assumem valores no intervalo [-1, 1].
± xi
wi
0,1488743389816312108
0,43339539412924719080
0,67940956829902440623
0,86506336668898451073
0,97390652851717172008
0,29552422471475287017
0,26926671930999635509
0,21908636251598204400
0,14945134915058059315
0,06667134430868813759
Tabela B1 – Valores de xi e wi para integração gaussiana, utilizando-se 10 pontos de Gauss
69
Anexo C
________________________
Identidades Vetoriais
Sejam a, b e c campos vetoriais e f e g campos escalares, então as seguintes
identidades vetoriais são válidas, Afonso, (2003):
a ⋅ (b × c ) = b ⋅ (c × a ) = c ⋅ (a × b );
(C.1)
a × b × c = (a ⋅ c )b − (a ⋅ b )c;
(C.2)
∇ ⋅ (a × b ) = b ⋅ ∇ × a − a ⋅ ∇ × b;
(C.3)
∇ ⋅ ( ga ) = ∇g ⋅ a + g∇ ⋅ a.
(C.4)
∇ ⋅ (b ∇ a) = ∇b ⋅ ∇a + b∇2 a
(C.5)
70
Anexo D
________________________
Função Delta de Dirac e Solução Fundamental
r r
A função delta de Dirac, δ (r , r ') representa a aplicação de quantidades (cargas,
potenciais, etc.) unitárias concentradas em um ponto
ponto fonte. O vetor
r
r'
r
r , chamado ponto de colocação ou
é um ponto do domínio, onde se realiza a integração. Essa função
pode ser definida como
v r
 0 para r ≠ r '
r r
δ i (r , r ') = 
v r
∞ para r = r '
(D.1)
Ela possui a possui a seguinte propriedade
r r
r
r
∫ δ (r , r ') φ (r ) dΩ = φ (r ) .
(D.2)
Ω
Dessa maneira, tem-se que se a função ponderadora escolhida u * utilizada na equação (4.2)
satisfaz a equação de Helmholtz, que governa o problema, e representa a solução para o
potencial de velocidade no ponto de colocação
fonte unitária aplicada no ponto
r
r , no problema de domínio infinito com uma
r
r ' , então u* é chamada de solução fundamental.
71
Anexo E
________________________
Nós e Elementos
A figura 5.2 apresenta a discretização do cilindro em 12 elementos. A tabela E.1 a
seguir, mostra os 12 elementos e os respectivos nós que os formam.
Nós
Elemento
1e2
1
2e3
2
3e4
3
4e5
4
5e6
5
6e7
6
7e8
7
8e9
8
9 e 10
9
10 e 11
10
11 e 12
11
12 e 1
12
Tabela E.1 – Nós e respectivos elementos
Para a implementação computacional, utilizou-se o número de elementos igual ao
número de nós.
72