ACOPLAMENTO MEC-MEF PARA PROBLEMAS ACÚSTICO

ACOPLAMENTO MEC-MEF PARA PROBLEMAS ACÚSTICO-

ELASTODINÂMICOS AXISSIMÉTRICOS NO DOMÍNIO DO TEMPO

Arnaldo Warszawski

DISSERTAÇÃO SUBMETIDA AO CORPO DOCENTE DA COORDENAÇÃO

DOS PROGRAMAS DE PÓS-GRADUAÇÃO DE ENGENHARIA DA

UNIVERSIDADE FEDERAL DO RIO DE JANEIRO COMO PARTE DOS

REQUISITOS NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE MESTRE

EM CIÊNCIAS EM ENGENHARIA CIVIL.

Aprovada por:

Prof. Webe João Mansur, Ph.D.

Prof. Delfim Soares Júnior, D.Sc.

Prof. José Antonio Marques Carrer, D.Sc.

Prof. José Antonio Fontes Santiago, D.Sc.

Prof. Ney Augusto Dumont, Dr.-Ing.

RIO DE JANEIRO, RJ - BRASIL

DEZEMBRO DE 2005

ii

WARSZAWSKI, ARNALDO

Acoplamento MEC-MEF para problemas

acústico-elastodinâmicos axissimétricos no

domínio do tempo [Rio de Janeiro] 2005

XIII, 145 p. 29,7 cm (COPPE/UFRJ, M.Sc.,

Engenharia Civil, 2005)

Dissertação – Universidade Federal do Rio

de Janeiro, COPPE

1. Propagação de ondas

2. Método dos Elementos de Contorno

3. Método dos Elementos Finitos

4. Interação fluido-estrutura

5. Axissimetria

I. COPPE/UFRJ II. Título (série)

iii

Agradecimentos

À minha família pelo apoio e incentivo em todos os momentos.

Ao prof. Webe João Mansur pelos conselhos e sugestões sempre úteis e pelos

ensinamentos passados.

Ao prof. Delfim Soares Júnior pela colaboração dada a este trabalho.

Aos colegas de mestrado pelo companheirismo.

Aos funcionários do LAMEC pela prontidão em ajudar.

Ao CNPq e à FAPERJ pelo auxílio financeiro.

iv

Resumo da Dissertação apresentada à COPPE/UFRJ como parte dos requisitos

necessários para a obtenção do grau de Mestre em Ciências (M. Sc.)

ACOPLAMENTO MEC-MEF PARA PROBLEMAS ACÚSTICO-

ELASTODINÂMICOS AXISSIMÉTRICOS NO DOMÍNIO DO TEMPO

Arnaldo Warszawski

Dezembro/2005

Orientadores: Webe João Mansur

Delfim Soares Júnior

Programa: Engenharia Civil

Neste trabalho é desenvolvido um código computacional para análise de

propagação de ondas em corpos axissimétricos que envolvem acoplamento acústico-

elástico, como é o caso, por exemplo, de problemas de interação fluido-estrutura. Esta

análise é feita sempre no domínio do tempo.

O meio acústico é modelado através do Método dos Elementos de Contorno

(MEC), cujas integrais no tempo são efetuadas analiticamente, utilizando o conceito

de parte finita de integral. Todas as singularidades presentes para a integração no

contorno são tratadas adequadamente. O meio elástico é modelado através do Método

dos Elementos Finitos (MEF), empregando-se o método de Newmark para avanço no

tempo.

O acoplamento entre os dois meios é feito na interface, através de

procedimento iterativo, preservando as características de resolução de cada meio.

São apresentadas algumas aplicações a fim de comprovar a validade das

expressões analíticas geradas para o MEC, bem como do algoritmo de acoplamento

implementado.

v

Abstract of Dissertation presented to COPPE/UFRJ as a partial fulfillment of the

requirements for the degree of Master of Science (M. Sc.)

BEM-FEM COUPLING FOR ACOUSTIC-ELASTODYNAMIC

AXISYMMETRIC PROBLEMS IN THE TIME DOMAIN

Arnaldo Warszawski

December/2005

Advisors: Webe João Mansur

Delfim Soares Júnior

Department: Civil Engineering

In this work, a computer code to model wave propagation in axisymmetric

bodies with acoustic-elastic coupling, as, for example, is the case of fluid-structure

interaction problems, is developed. This analysis is always processed in time domain.

The acoustic medium is modeled by the Boundary Element Method (BEM),

whose time integrals are evaluated analytically, using the concept of finite part

integrals. All the singularities for space integration are treated adequately. The elastic

medium is modeled by the Finite Element Method (FEM), employing the Newmark

method for time marching.

The coupling between the two media is carried out on interfaces, through an

iterative procedure, preserving the features of the solution of each medium.

Some applications are presented in order to prove the validity of the

expressions generated for the BEM, as well as of the implemented coupling

algorithm.

vi

Índice Introdução ......................................................................................................................1

I.1. A propagação de ondas na engenharia ..............................................................1

I.2. Conceitos básicos e breve revisão bibliográfica................................................2

I.3. Objetivos e resumo dos capítulos ......................................................................6

1. Formulação do MEC para problemas acústicos axissimétricos.................................8

1.1. Equação integral e solução fundamental 3-D ...................................................9

1.2. Solução fundamental axissimétrica.................................................................13

2. Implementação numérica do MEC ..........................................................................16

2.1. Discretização e aproximações.........................................................................17

2.1.1. Funções de interpolação.........................................................................18

2.1.2. Montagem do sistema de equações........................................................20

2.1.3. Propriedade de translação ......................................................................22

2.2. Integração analítica no tempo .........................................................................24

2.2.1. Ponto fonte fora do eixo de axissimetria................................................24

2.2.2. Ponto fonte pertencente ao eixo de axissimetria....................................33

2.3. Integração no espaço e singularidades............................................................35

2.3.1. Singularidades em r = 0 .........................................................................36

2.3.2. Singularidades nas frentes de onda da solução fundamental .................42

2.3.3. Singularidades para ponto campo no eixo de axissimetria ....................45

2.3.4. Transformada polinomial de 3o grau para a integração .........................46

2.4. Termos relativos a fontes pontuais .................................................................48

2.5. Esquema Teta..................................................................................................51

3. Acoplamento MEC-MEF.........................................................................................53

3.1. Expressões obtidas para o MEF......................................................................54

3.1.1. Meio acústico .........................................................................................54

3.1.2. Meio elástico..........................................................................................59

3.2. Equações governantes do acoplamento ..........................................................64

3.2.1. Acoplamento acústico-acústico MEC-MEF ..........................................64

3.2.2. Acoplamento acústico-elástico MEC-MEF ...........................................65

3.3. Procedimentos numéricos para o acoplamento iterativo ................................67

3.3.1. Acoplamento acústico-acústico MEC-MEF ..........................................67

3.3.2. Acoplamento acústico-elástico MEC-MEF ...........................................70

vii

4. Aplicações................................................................................................................72

4.1. Meio discretizado exclusivamente pelo MEC ................................................72

4.1.1. Barra prismática circular submetida a fluxo não-nulo em uma

extremidade......................................................................................................73

4.1.2. Cavidade em meio acústico infinito.......................................................87

4.1.3. Membrana submetida a carga impulsiva................................................91

4.1.4. Barra prismática circular submetida a pressão não-nula em uma

extremidade......................................................................................................96

4.2. Acoplamento acústico-acústico ......................................................................98

4.3. Acoplamento acústico-elástico .....................................................................101

4.3.1. Barra prismática circular vazada..........................................................101

4.3.2. Fonte emissora no interior de um duto metálico..................................107

4.3.3. Parede circundada por fluido submetida a carga impulsiva.................112

Conclusões .................................................................................................................116

Bibliografia ................................................................................................................120

A. Expressões analíticas para as integrais no tempo..................................................124

A.1. Expressões para g.........................................................................................124

A.2. Expressões para hI e hF ................................................................................125

B. Integrais elípticas...................................................................................................130

B.1. Tipos de integrais elípticas...........................................................................130

B.2. Cálculo das integrais ....................................................................................131

C. Parte finita de integrais..........................................................................................133

D. Expressões para integrais das singularidades na frente de onda...........................136

D.1. Singularidade em r .......................................................................................136

D.2. Singularidade em d.......................................................................................137

E. Expressões analíticas para os termos relativos a fontes ........................................140

F. Soluções analíticas de interesse .............................................................................142

F.1. Barra engastada em uma extremidade e livre na outra, submetida a

carregamento axial de impacto ............................................................................142

F.2. Cavidade esférica em meio infinito ..............................................................143

F.3. Membrana circular submetida a um campo de velocidades iniciais.............143

F.4. Barra acústica fechada em uma extremidade e livre na outra, submetida a

pressão de impacto...............................................................................................145

viii

Índice de figuras

Figura 1.1 – variáveis necessárias para solução axissimétrica ....................................13

Figura 2.1 – (a) contorno original; (b) contorno discretizado por elementos lineares.17

Figura 2.2 – funções de interpolação no tempo φ m(t) e θ m(t) .....................................19

Figura 2.3 – funções de interpolação no espaço η j(X) e ν

j(X) ao longo dos elementos

adjacentes ao nó j .........................................................................................................20

Figura 2.4 – exemplificação da propriedade de translação..........................................22

Figura 2.5 – região onde é necessário efetuar as integrais de convolução em cada

passo de tempo, com suas respectivas funções de interpolação ..................................23

Figura 2.6 – representação dos seis casos de limites de integração.............................27

Figura 2.7 – gráfico da diferença entre as integrais elípticas do primeiro tipo

incompleta e completa, para ui = 0,6; ξ1 = 1,0; elemento a 45o em relação à horizontal

(0 ≤ r ≤ ui – trecho cujos pontos se enquadram no caso 2 – expressão (A.1)).............36

Figura 2.8 – produto da função de interpolação r/l pela integral elíptica incompleta de

primeiro tipo, com os mesmos dados numéricos da figura 2.7....................................37

Figura 2.9 – gráfico da diferença entre os dois membros da equação (2.47), para os

mesmos dados numéricos do gráfico da figura 2.7......................................................38

Figura 2.10 – integrando da primeira parcela do segundo membro de (2.51) .............39

Figura 2.11 – gráfico da diferença entre integrais elípticas existente em (A.3), (A.9) e

(A.17), para os mesmos dados numéricos das figuras 2.7 a 2.9, e uf = 0,2, traçado para

0 ≤ r ≤ uf e d ≥ ui – trecho cujos pontos se enquadram no caso 4................................41

Figura 3.1 – extrapolação da pressão obtida pelo MEF...............................................68

Figura 3.2 – interpolação do fluxo obtido pelo MEC ..................................................69

Figura 4.1 – (a) esquema do contorno axissimétrico; (b) seção transversal, para

situações 1, 2 e 3; (c) malha empregada na análise .....................................................74

Figura 4.2 – comparação entre o resultado da situação 1 (em preto) e o analítico (em

vermelho) .....................................................................................................................77

Figura 4.3 – comparação entre o resultado das situações 2 (em preto, linha cheia) e 3

(em preto, linha pontilhada) e o analítico (em vermelho)............................................79

Figura 4.4 – pressão no ponto A para análise estendida ..............................................82

Figura 4.5 – pressão no ponto A para análise estendida, situação 2 – comparação entre

amortecimentos numéricos ..........................................................................................83

ix

Figura 4.6 – variação da resposta para pressão no ponto A com a mudança de θ.......84

Figura 4.7 – resposta para pressão no ponto A para integração no tempo numérica e

analítica ........................................................................................................................86

Figura 4.8 – (a) esquema da cavidade esférica em meio infinito; (b) malha de oito

elementos empregada...................................................................................................87

Figura 4.9 – resultados para 8 elementos e θ = 1,0 – pontos do contorno...................88

Figura 4.10 – resultados para o ponto 1, variando-se o valor de θ ..............................89

Figura 4.11 – comparação entre malhas ......................................................................89

Figura 4.12 – pontos internos (em preto), comparação com resultado analítico (em

vermelho) .....................................................................................................................90

Figura 4.13 – esquema do modelo da membrana: (a) em perfil; (b) em planta. ..........92

Figura 4.14 – gráfico no tempo da carga aplicada .......................................................92

Figura 4.15 – comparação de espessuras: (a) e = 0,1m; 0,2 m e teórico; (b) e = 1,0 m

......................................................................................................................................93

Figura 4.16 – comparação entre resultado numérico e analítico – situação 1 .............94

Figura 4.17 – derivada normal ao contorno direito – situação 1 .................................95

Figura 4.18 – derivada normal ao contorno direito – situação 2 .................................95

Figura 4.19 – pressão em um ponto a meia barra ........................................................96

Figura 4.20 – fluxo na extremidade da barra ...............................................................97

Figura 4.21 – modelo para acoplamento acústico-acústico .........................................98

Figura 4.22 – pressão no ponto A – acoplamento x MDF...........................................99

Figura 4.23 – pressão no ponto B – acoplamento x MDF .........................................100

Figura 4.24 – esquema da aplicação do item 4.3.1: (a) situação 1; (b) situação 2; (c)

malha empregada; cotas em metros. ..........................................................................102

Figura 4.25 – deslocamentos uz em A (a) e B (b) para tolerâncias de 10-3 e 10-6......103

Figura 4.26 – deslocamentos uz em A (a) e B (b) variando passo de tempo do MEF103

Figura 4.27 – deslocamentos uz nos pontos A (a) e B (b) variando refinamento da

malha do MEF............................................................................................................104

Figura 4.28 – comparação entre resultados analíticos e numéricos - situação 1 .......105

Figura 4.29 – comparação entre resultados analíticos e numéricos - situação 2 .......106

Figura 4.30 – modelo de um meio afetado por uma fonte geradora de pressão; cotas

em metros...................................................................................................................107

Figura 4.31 – pressão no ponto A – comparação entre métodos ...............................108

x

Figura 4.32 – pressão no ponto B – comparação entre métodos ...............................109

Figura 4.33 – deslocamento horizontal no ponto C – comparação entre métodos ....109

Figura 4.34 – pressão em A – comparação com e sem o duto...................................110

Figura 4.35 – pressão em B – comparação com e sem o duto ...................................110

Figura 4.36 – esquema do modelo para aplicação do item 4.3.3; cotas em metros...112

Figura 4.37 – deslocamento horizontal em A – comparação entre métodos .............113

Figura 4.38 – deslocamento horizontal em B – comparação entre métodos .............113

Figura 4.39 – deslocamento horizontal em A – comparação com e sem fluido ........114

Figura 4.40 – deslocamento horizontal em B – comparação com e sem fluido ........114

Figura F.1 – barra engastada submetida a carga normal............................................142

Figura F.2 – cavidade esférica em meio infinito........................................................143

Figura F.3 – membrana circular.................................................................................144

xi

Lista de símbolos

Símbolos romanos

a coeficiente da transformação de coordenadas.

A função auxiliar para integração de partes finitas.

A matriz de coeficientes do sistema do MEC.

A matriz efetiva do MEF.

b coeficiente da transformação de coordenadas.

B função auxiliar para integração de partes finitas.

B matriz de formação do vetor de condições de contorno do MEC; matriz

de derivadas das funções de interpolação do MEF.

B vetor efetivo do MEF.

c velocidade de propagação da onda primária; coeficiente geométrico do

MEC; coeficiente da transformação de coordenadas.

C1, C2 coeficientes auxiliares para integração.

C matriz de amortecimento do MEF.

d maior distância entre o ponto campo e o anel de fontes; coeficiente da

transformação de coordenadas.

D matriz constitutiva.

E integral elíptica do 2o tipo completa ou incompleta; módulo de

elasticidade.

f função descrita pela emissão de uma fonte pontual.

F integral elíptica incompleta do 1o tipo.

F vetor de cargas nodais do MEF.

g expressão obtida da integração analítica no tempo para montagem de

G.

G jacobiano da transformação de coordenadas.

G matriz de convolução.

h expressão obtida da integração analítica no tempo para montagem de

H; menor distância entre o ponto fonte e a reta que contém o elemento onde está o

ponto campo.

H função Heaviside.

H matriz de convolução.

xii

K integral elíptica completa do 1o tipo; módulo de compressibilidade.

K matriz de rigidez do MEF.

l comprimento.

L matriz de transformação deslocamento-deformação.

M matriz de massa do MEF.

N função de interpolação do MEF.

n vetor normal ao contorno.

p pressão; coeficiente da transformação de coordenadas.

p* pressão da solução fundamental.

p vetor de pressões em certo instante de tempo.

q fluxo; coeficiente da transformação de coordenadas.

q* fluxo da solução fundamental.

q vetor de fluxos em certo instante de tempo.

Q~ vetor que contém fluxos nodais.

r distância entre o ponto fonte no plano de axissimetria e o ponto campo.

R distância entre o ponto fonte e o ponto campo.

S vetor com a contribuição das fontes pontuais.

t tempo.

u variável auxiliar para integração; deslocamento.

U vetor de deslocamentos.

v derivada da pressão em relação ao tempo.

w expressão obtida da integração analítica no tempo do termo relativo a

fontes pontuais; função de ponderação.

X ponto campo.

x vetor de incógnitas do MEC.

y vetor com condições de contorno do MEC.

Símbolos gregos

α parâmetro de relaxação.

β parâmetro de refinamento do MEC; parâmetro de Newmark.

γ densidade de fontes; parâmetro de Newmark; coordenada natural dos

pontos de Gauss.

xiii

Γ contorno do meio.

δ função Delta de Dirac.

ε parâmetro de convergência das iterações.

ζ coeficiente de amortecimento viscoso.

η função de interpolação espacial de p; coordenada natural dos pontos de

integração transformados.

θ função de interpolação temporal de q; ângulo definidor de direção do

elemento de contorno; parâmetro para estabilização do MEC.

ν função de interpolação espacial de q; coeficiente de Poisson.

ξ ponto fonte; coordenada natural.

ρ densidade.

τ variável temporal.

Φ função de interpolação temporal de p.

Ω domínio do meio.

Abreviaturas

MEC Método dos Elementos de Contorno.

MEF Método dos Elementos Finitos.

Notações

Yi termo i do vetor Y (Y = (Y1, Y2, Y3) em 3D e Y = (Y1, Y2) em 2D).

Y& derivada temporal de Y ( tYY ∂∂= /& ).

∇Y vetor gradiente de Y (∇Y = (Y,1; Y,2; Y,3) em 3D e ∇Y = (Y,1; Y,2) em

2D).

∇2 operador de Laplace (∇2Y = Y,ii).

Y norma do vetor Y ( ii YYY .= )

Introdução 1

Introdução

I.1. A propagação de ondas na engenharia

O estudo do fenômeno da propagação de ondas vem sendo aprimorado há

séculos, constituindo-se de inúmeras subdivisões, entre as quais podem-se destacar a

ótica (propagação da luz) e a acústica (propagação do som). Estes dois campos são de

extrema importância, pela sua participação na comunicação de seres humanos e

inúmeros animais. Dois dos cinco sentidos do ser humano devem sua existência às

ondas luminosas e sonoras. Além destas, diversos outros tipos de ondas mecânicas e

eletromagnéticas (sísmicas, ultrassonoras, microondas) têm um vasto campo de

aplicação na ciência.

No contexto da engenharia atual, a análise da propagação de ondas ganha

cada vez mais importância, em virtude da diversidade de aplicações dela

provenientes. O desenvolvimento de equipamentos capazes de emitir ondas, geradas

através de sinais, bem como de receber estas ondas e decodificá-las novamente em

sinais, possibilita a ampliação cada vez maior deste campo de aplicações.

Na indústria do petróleo, pode-se destacar primeiramente a aplicação do

estudo de propagação de ondas na realização de levantamentos sísmicos, inclusive em

regiões com lâmina d’água elevada. Ondas acústicas emitidas na superfície da água

(ou do solo), ao encontrarem obstáculos ou mudanças de propriedades entre meios de

propagação, refletem-se, fazendo com que parte da energia emitida retorne à

superfície. As ondas refletidas são medidas em receptores também localizados na

superfície, rebocados por navios e, através do procedimento conhecido como

migração em profundidade, é possível, através das respostas medidas, identificar

camadas, intrusões salinas e, principalmente, regiões contendo hidrocarbonetos. Esse

estudo assume uma grande responsabilidade em vista dos vultosos custos envolvidos

na perfuração de poços de petróleo.

Uma aplicação comum em engenharia civil é a que utiliza GPR (“ground

penetrating radar”), a fim de se investigar a composição geológica do subsolo para

profundidades pequenas. Procedimentos semelhantes aos expostos no parágrafo

anterior permitem obter informações a respeito das camadas presentes no subsolo. O

GPR, apesar de ainda não comumente empregado, pode ser uma alternativa aos

Introdução 2 métodos tradicionais de investigação do subsolo, quando não for vantajoso, por

diversos motivos, que estes sejam utilizados, podendo ser útil também como

ferramenta complementar. Para profundidades pequenas (até 100 metros), como é o

caso em obras civis, ondas eletromagnéticas são em geral utilizadas. Em

profundidades maiores, utilizam-se ondas sísmicas, que viajam distâncias maiores no

subsolo, mas não são capazes de identificar obstáculos pequenos.

Uma última aplicação mencionada aqui é na engenharia estrutural, pelo

desenvolvimento de equipamentos que possibilitam a identificação de danos (fissuras,

alteração de propriedades) em estruturas de concreto, também através da medição da

resposta a uma determinada fonte emitindo ondas na superfície da estrutura. Essa

identificação é de extrema importância em estruturas como barragens, em que a

presença de fissuras pode gerar graves conseqüências, pela possibilidade de

propagação das mesmas. Deste modo, pode-se reduzir a necessidade da retirada de

testemunhos, que explicitariam a existência da fissura.

Existe uma variedade extensa de problemas cujo entendimento requer a

utilização de metodologias baseadas nos fenômenos de propagação de ondas; um

deles é o problema de interação acústico-elástica (fluido-estrutura) estudado nesta

monografia.

I.2. Conceitos básicos e breve revisão bibliográfica

De acordo com o material de que é composto o meio onde ocorre a

propagação de ondas, pode-se fazer a distinção entre meios elásticos e acústicos. A

princípio, POISSON [1] foi o primeiro a identificar que a onda que se propaga em um

meio elástico é composta de uma onda primária, também chamada de irrotacional e

longitudinal, e uma onda secundária, também chamada de transversal, distorcional e

de cisalhamento. Um meio acústico, essencialmente constituído por fluidos, não

resiste ao cisalhamento, por isso não ocorre nele a propagação de ondas secundárias.

Este conceito é de fundamental importância, já que, dependendo do material de que é

composto o meio em estudo, a modelagem física é completamente diferente.

Com o aumento da diversidade de aplicações, o estudo da propagação de

ondas cresceu bastante nas últimas décadas, impulsionado também pelo aumento da

capacidade de processamento e memória dos microcomputadores, o que possibilitou o

Introdução 3 desenvolvimento de métodos mais sofisticados e eficientes visando à modelagem

numérica do problema.

A modelagem de problemas transientes pode ser realizada no domínio do

tempo ou no domínio da freqüência. As análises baseadas no domínio do tempo são

mais familiares, visto que a resolução se processa em passos de tempo,

correspondentes ao tempo real de análise (BATHE [2], COOK [3]). As análises no

domínio da freqüência necessitam de domínios transformados, tendo se desenvolvido

bastante ultimamente, já que em alguns casos constituem-se no modo mais adequado

para análise (THOMSON [4], CLOUGH e PENZIEN [5]). Neste trabalho, todo o

desenvolvimento é feito no domínio do tempo.

O Método dos Elementos Finitos (MEF) é provavelmente o método numérico

mais popular tanto nos meios acadêmicos quanto em empresas produtoras de

programas computacionais comerciais. Este método é mais adequado para problemas

relativos a meios não-homogêneos e anisotrópicos, bem como para lidar com as não-

linearidades passíveis de existência no problema. Mais detalhes sobre o método

podem ser encontrados em ZIENKIEWICZ e TAYLOR [6] e BATHE [2].

O MEF, no entanto, apresenta deficiências quando se lida com domínios

infinitos, já que há a necessidade de limitar a malha. Deste modo, são introduzidos

contornos artificiais. Tem-se tentado desenvolver contornos chamados de

transmissores, ou não-refletores, que numericamente eliminem reflexões, absorvendo

a energia, a fim de representar corretamente a continuidade do meio (GIVOLI e

NETA [7], SARMA et al. [8]), porém sua eficiência ainda é questionável. Com isso,

muitas vezes é necessário que a malha de elementos finitos seja estendida a uma

dimensão muito grande, sendo que este procedimento em alguns casos inviabiliza as

análises, produzindo malhas de tamanho proibitivo.

O Método dos Elementos de Contorno (MEC) teve seu desenvolvimento

iniciado mais recentemente, sendo derivado da discretização numérica de equações

integrais obtidas a partir das equações físicas básicas do problema em questão.

Através do MEC um corpo de domínio infinito pode ser corretamente representado, já

que o método exige apenas a discretização do contorno. Mais detalhes sobre o método

podem ser encontrados em DOMINGUEZ [9], BECKER [10] e MANSUR [11]. Uma

desvantagem deste método é que, dependendo do vulto do problema, há a

necessidade, na análise no domínio do tempo, do armazenamento de grande

Introdução 4 quantidade de matrizes de convolução correspondentes aos instantes de tempo de

análise decorridos.

Em alguns tipos de problema, o domínio em análise é formado por sistemas de

diferentes características interagindo entre si. Neste caso, não é possível a solução de

cada sistema isoladamente, já que os resultados de cada um influenciam a resposta do

outro. Estes sistemas são ditos acoplados e, dependendo das suas características, o

acoplamento pode ser do tipo acústico-acústico, acústico-elástico ou elástico-elástico,

por exemplo.

Segundo SOARES JR. [12], os sistemas acoplados podem interagir apenas em

suas interfaces ou seus domínios podem se sobrepor. No primeiro caso, chamado de

acoplamento de interface, o acoplamento se dá através das condições de contorno na

interface, enquanto no segundo, chamado de acoplamento de domínio, o acoplamento

se dá através das equações governantes dos fenômenos físicos. No presente trabalho,

lida-se apenas com acoplamento de interface. Como exemplo deste tipo de

acoplamento, pode-se citar a interação entre fluido e estrutura, como no caso de água

interagindo com estruturas de barragens, portos e estruturas off-shore, ou de

lubrificantes interagindo com elementos de máquinas; outro caso é o de interação

entre sólidos de propriedades diferentes, como entre o solo e as fundações de uma

estrutura.

Os diferentes sistemas que interagem entre si podem ser modelados por

métodos numéricos distintos. Neste trabalho trata-se de sistemas acoplados modelados

pelo MEF e pelo MEC. O estudo deste tipo de acoplamento leva, entre outros, aos

trabalhos de VON ESTORFF e ANTES [13], BELYTSCHKO e LU [14] e YU et al.

[15]. Nestes trabalhos, a resolução do sistema ocorre através da formação de um

sistema de equações acoplado. Segundo SOARES JR. e VON ESTORFF [16], alguns

problemas advêm deste procedimento, sendo importante citar no contexto do presente

trabalho os seguintes: o sistema acoplado de equações não apresenta as características

de esparsidade e configuração em banda daquele obtido no MEF, inviabilizando o

emprego dos programas de solução já desenvolvidos para os sistemas obtidos neste

método; e o intervalo de tempo da solução deve ser único para todas as malhas. Esse

último é um problema grave, visto que há casos em que a diferença entre as

velocidades de propagação de onda dos diversos meios acoplados é considerável, o

Introdução 5 que exigiria discretizações espaciais e temporais distintas, adequadas a cada meio e a

cada método empregado na modelagem.

Estes problemas são evitados empregando-se o método iterativo para

acoplamento apresentado inicialmente por SOARES et al. [17] para análise de

sistemas não-lineares bidimensionais com acoplamento do tipo sólido-sólido.

Segundo esta abordagem, cada sistema é resolvido separadamente, sendo verificadas

as condições de contorno na interface, através de um processo iterativo, até que estas

sejam atingidas (tal procedimento é descrito detalhadamente no item 3.3 do presente

trabalho, adaptado ao caso axissimétrico). Posteriormente, o conceito de acoplamento

iterativo MEC-MEF foi expandido por SOARES et al. [18] para análise de sistemas

acoplados do tipo sólido-fluido, tendo mais uma vez como foco problemas

bidimensionais.

Este trabalho trata de problemas axissimétricos, ou seja, que possuem simetria

em relação a um eixo, que se configura como um eixo de revolução. A consideração

desta simetria representa uma economia, já que evita a análise do problema

tridimensional. Quando for feita referência a um problema axissimétrico, fica

subentendido a partir de agora que está se tratando tanto da geometria quanto do

carregamento aplicado.

A análise de meios acústicos axissimétricos modelados pelo MEC, um ponto

fundamental neste trabalho, inicialmente concentrou-se em soluções no domínio da

freqüência. GRANNELL et al. [19] deduziram um método de alta precisão para

solução do problema de Neumann para a equação de Helmholtz. TSINOPOULOS et

al. [20] e SOENARKO [21] lidaram com problemas de propagação de ondas

acústicas com condições de contorno não-axissimétricas. No domínio do tempo, o

usual era a integração numérica da solução fundamental tridimensional em torno do

eixo de axissimetria, como em ISRAIL et al. [22] e BESKOS [23]. No entanto, em

CZYGAN e VON ESTORFF [24] é apresentado o procedimento de integração

analítica da solução fundamental em torno do eixo de axissimetria, chegando-se então

à solução fundamental axissimétrica. No presente trabalho, o principal foco consiste

na integração analítica ao longo do tempo desta solução.

Introdução 6

I.3. Objetivos e resumo dos capítulos

No Capítulo 1, apresenta-se a formulação do MEC para problemas de

propagação de ondas acústicas em meios homogêneos. Partindo-se da equação

governante do fenômeno, chega-se à equação integral de contorno através da qual se

obtém a pressão em qualquer ponto do domínio ou do contorno em estudo. É

apresentado resumidamente o processo empregado por MANSUR [11] para se obter

esta equação, o qual se utiliza do método dos resíduos ponderados, considerando

como função de ponderação a solução fundamental do problema. Este procedimento

refere-se ao caso genérico tridimensional. Em seguida, é mostrado brevemente o

procedimento utilizado por CZYGAN e VON ESTORFF [24] para obtenção da

solução fundamental axissimétrica.

O Capítulo 2 contém todo o procedimento aplicado às integrais de contorno a

fim de resolver o problema numericamente. O contorno é discretizado em elementos

de contorno de geometria linear e são apresentadas as funções de interpolação

empregadas no método. A equação integral de contorno é escrita para todos os pontos

do contorno para se chegar ao sistema de equações de interesse. Os coeficientes das

matrizes deste sistema são provenientes de integrais no tempo (relativas à

convolução) e no espaço. O próximo passo, então, é a integração analítica no tempo,

cujo procedimento é detalhado no item 2.2. Os cuidados a serem tomados para a

integração espacial são também aprofundados, no item seguinte, no qual se apresenta

o processo semi-analítico empregado a fim de se eliminar as singularidades presentes

nas expressões. No item 2.4, trata-se da integração analítica no tempo dos termos

correspondentes a fontes pontuais presentes no meio. O “Esquema Teta” utilizado

para estabilização da resposta é brevemente mostrado no fim do capítulo. As

integrações analíticas, que produzem as expressões apresentadas nos Anexos A, D e

F, são o cerne do trabalho, sendo o resultado obtido através deste processo comparado

com aquele obtido através de integração numérica.

No Capítulo 3, inicialmente se apresenta a formulação do MEF, tanto para

problemas acústicos quanto para problemas elásticos, sendo mostrada a adaptação

realizada para problemas axissimétricos. Em seguida é abordado o acoplamento de

sistemas fisicamente distintos (acústico-elástico), ou fisicamente similares (acústico-

acústico), porém modelados por métodos distintos. São relacionadas as condições a

Introdução 7 serem obedecidas na interface para o acoplamento. O método iterativo para

implementação numérica do acoplamento é apresentado no item 3.3.

O Capítulo 4 apresenta exemplos de aplicação dos métodos descritos nos

capítulos anteriores visando à validação do que foi desenvolvido aqui. Os exemplos,

primeiramente, referem-se a problemas de meio exclusivamente acústico, a fim de se

testar as expressões analíticas obtidas no Capítulo 2 e mostradas nos Anexos A, D e

E. Posteriormente, passa-se para exemplos de domínios em que ocorre interação entre

meios acústicos, e entre meios acústicos e elásticos, verificando-se o acoplamento

implementado.

O objetivo principal do trabalho é a obtenção das expressões analíticas

provenientes da integração ao longo do tempo. Com isso, é desenvolvido um código

computacional em Fortran para análise de meios acústicos axissimétricos pelo MEC.

Posteriormente, trabalha-se no acoplamento iterativo para meios

axissimétricos. Utilizando-se de um programa já desenvolvido para solução de

problemas acústicos e elastodinâmicos pelo MEF, promove-se a interação entre os

dois programas através de uma sub-rotina de acoplamento de interface e de execução

das iterações.

Formulação do MEC para problemas acústicos axissimétricos 8

Capítulo 1

Formulação do MEC para problemas acústicos axissimétricos

A equação escalar da onda, explicitada pela expressão (1.1), governa a

propagação de ondas através de meios acústicos, de modo geral, podendo ser também

empregada na análise de diversos outros problemas, tais como vibração longitudinal

de barras, vibração transversal de cordas e de membranas. Em alguns casos, um

problema de meio elástico pode ter suas equações governantes simplificadas, de modo

a recaírem na equação escalar da onda, ampliando o campo de aplicações dos

procedimentos de resolução desta equação.

γ−=∂∂

−∇ 2

2

22 1

tp

cp (1.1)

Em (1.1), p é a pressão em qualquer ponto do domínio do problema, c é a

velocidade de propagação da onda acústica no meio e γ a função que descreve a

distribuição espacial e temporal de fontes no domínio. O operador 2∇ representa o

Laplaciano de uma função. O domínio do problema é denominado Ω, sendo seu

contorno expresso por Γ.

São consideradas as condições iniciais de pressão e de sua derivada em relação

ao tempo expressas em (1.2), onde o ponto X representa um ponto qualquer do

domínio.

( )

( ) )(0,

)(0,

0

0

XvXtp

XpXp

=∂∂

= (1.2)

As condições de contorno do problema são expressas por (1.3), onde Γ1 é a

parte do contorno em que são prescritas pressões e Γ2 a parte em que são prescritas as

derivadas na direção normal ao contorno, que equivalem ao fluxo q.

( )

( ) 2

1

,,

,,

Γ∈=∂∂

Γ∈=

XqtXnp

XptXp (1.3)


Este capítulo tem como objetivo a obtenção, através da manipulação da

equação (1.1), da equação integral a ser empregada para a análise pelo Método dos

Elementos de Contorno (MEC). Os procedimentos são aqui expostos de forma

simplificada, já que não é este o objetivo do trabalho. A dedução mais rigorosa pode

ser encontrada em MANSUR [11].

Obtida a equação integral, necessita-se da solução fundamental, que, como

será visto, é uma função de Green do problema. Por fim, a solução fundamental geral,

obtida para o caso tridimensional, é particularizada para o caso em que se lida com a

axissimetria do problema, conforme executado em CZYGAN e VON ESTORFF [24].

Os resultados aqui expostos são a base teórica para os procedimentos

implementados no capítulo seguinte.

1.1. Equação integral e solução fundamental 3-D

Através da utilização do método dos resíduos ponderados, a partir da equação

de interesse (1.1), é escrita a expressão (1.4), na qual τ é a variável ao longo do tempo.

( ) ( )∫∫∫∫∫∫+++

ΓΓΩ

Γ−−Γ−=Ω⎟⎟⎠

⎞⎜⎜⎝

⎛+

∂∂

−∇

ttt

ddqppddpqqddppc

p0

*

0

*

0

*2

2

22

12

1 τττγτ

(1.4)

A função de ponderação p* é chamada de solução fundamental do problema,

apresentada mais adiante. Esta função deve ser solução da equação (1.1), em uma

região descrita pelo domínio Ω*, de contorno Γ*, de modo que Ω* contenha o domínio

e o contorno originais. A função q* representa a derivada da solução fundamental na

direção normal ao contorno, ou seja, o fluxo no contorno relacionado à solução

fundamental. A integração é efetuada até um instante de tempo t+, que corresponde ao

instante t acrescido de um parâmetro infinitesimal. Este procedimento é necessário

para que o limite de integração não coincida com o pico de uma função delta de

Dirac.

O próximo passo consiste em aplicar a segunda identidade de Green,

apresentada em (1.5), ao termo de (1.4) que contém o Laplaciano. Esta segunda

Formulação do MEC para problemas acústicos axissimétricos 10 identidade de Green pode ser obtida através de manipulações do teorema da

divergência [11].

( )∫∫∫ΓΩΩ

Γ−+Ω∇=Ω∇ dpqqpdppdpp ***2*2 (1.5)

Ao mesmo tempo, realiza-se, duas vezes a integração por partes no termo que

contém a derivada segunda da pressão em relação ao tempo, expressa em (1.6).

∫∫+

++

∂∂

+⎥⎦

⎤⎢⎣

⎡∂∂

−∂∂

=∂∂

ttt

dppppppdpp

0

2

*2

0

**

0

*2

2

ττττ

ττ

(1.6)

Com isso, chega-se à expressão (1.7), na qual se considera que a pressão e o

fluxo são perfeitamente conhecidos, respectivamente, em Γ1 e Γ2.

( )

01

1

0

**

2

0

*

0

2

*2

2*2

0

**

=Ω⎥⎦

⎤⎢⎣

⎡∂∂

−∂∂

−Γ+

+Ω⎟⎟⎠

⎞⎜⎜⎝

⎛∂

∂−∇+Γ−

∫∫∫

∫∫∫∫

ΩΩ

ΩΓ

++

++

dppppc

ddp

ddppc

pddpqqp

tt

tt

τττγ

ττ

τ

(1.7)

Sabendo-se que, devido à propriedade da causalidade, é válida a expressão

(1.8), chega-se a (1.9).

0**

=⎥⎦⎤

⎢⎣⎡∂∂

=⎥⎦

⎤⎢⎣

⎡∂∂

++ == tt

pppp

ττ ττ (1.8)

( )

( ) 01

1

*000

*02

0

*

0

2

*2

2*2

0

**

=Ω−−Γ+

+Ω⎟⎟⎠

⎞⎜⎜⎝

⎛∂

∂−∇+Γ−

∫∫∫

∫∫∫∫

ΩΩ

ΩΓ+

++

dpvpvc

ddp

ddppc

pddpqqp

t

tt

τγ

ττ

τ

(1.9)


A solução fundamental p* aqui utilizada é a função de Green para um meio

infinito submetido a uma fonte concentrada dada pela expressão (1.10). O ponto ξ é

chamado de ponto fonte, no qual está localizada a fonte geradora de pressão da

solução fundamental, sendo o ponto X denominado ponto campo, que corresponde a

um ponto qualquer no qual se está medindo a resposta devida à excitação provocada

pela fonte. A fonte corresponde a uma função Delta de Dirac também no tempo,

sendo sua aplicação dada no tempo τ, enquanto a resposta é medida em um tempo

qualquer t.

( ) ( )τδξδγ −−= tX* (1.10)

Substituindo-se a expressão para a fonte em (1.1), pode-se chegar à solução

fundamental dada por (1.11), segundo MORSE e FESHBACH [25].

( )[ ]RtcR

cp −−= τδπ4

* (1.11)

( )[ ] ( )[ ]⎟⎠⎞

⎜⎝⎛ −−+−−−

∂∂

=∂∂

= RtcRtcRc

Rnr

npq τδτδ

π&

41*

* (1.12)

Em (1.11), R é a distância entre o ponto fonte ξ e o ponto campo X.

Substituindo-se o termo entre parênteses da segunda parcela de (1.9) pela expressão

(1.10) e se utilizando das propriedades da integração de uma função Delta de Dirac,

chega-se finalmente à expressão (1.13).

( ) ( ) ( ) ( ) ( )

( ) ( ) ( ) ( )

( ) ( )∫ ∫

∫∫

∫ ∫∫ ∫

+

++

Ω

ΩΩ

ΓΓ

Γ+

+Ω+Ω−

−Γ−Γ=

t

tt

ddXtXp

dXvtXpc

dXptXvc

ddXptXqddXqtXptp

0

*

0*020

*02

0

*

0

*

,,,,

,,1,,1

,,,,,,,,,

ττγτξ

ξξ

τττξτττξξ

(1.13)


A expressão (1.13) é válida para um ponto ξ não pertencente ao contorno. No

caso de este ponto pertencer ao contorno, o primeiro termo do lado esquerdo de (1.13)

é multiplicado pelo coeficiente c(ξ), dado por (1.14).

( )

( )

( )∫∫

Γ

Γ

→Γ

Γ

=

ε

ε

ε

ε

εξ

ξ

ξdXq

dXq

imlc,

,

*

*

0 (1.14)

Em (1.14), εΓ é o contorno de um domínio εΩ , que corresponde a uma esfera

de raio infinitesimal ε em torno do ponto ξ, sendo εΓ a parcela do contorno de εΩ

localizada no interior do domínio original do problema Ω.

Para contornos tridimensionais, quando o contorno não é suave em ξ, o cálculo

pode se tornar bastante trabalhoso, portanto, evita-se colocar pontos fontes em pontos

angulosos de contornos tridimensionais. Sendo o contorno suave em ξ, não

apresentando “bico”, c(ξ) = 0,5.

Se o domínio modelado é caracterizado por duas dimensões (problemas

bidimensionais ou axissimétricos), εΩ não mais é uma esfera (torna-se um círculo,

representando um cilindro infinito, no problema bidimensional, e um toróide, no

problema axissimétrico), e a expressão (1.14) se simplifica em (1.15), sendo α o

ângulo interno formado pelas tangentes ao contorno do modelo no ponto ξ.

( )π

αξ2

=c (1.15)


1.2. Solução fundamental axissimétrica

Neste item é sucintamente exposto o procedimento apresentado por CZYGAN

e VON ESTORFF [24] para se chegar à solução fundamental para problemas

acústicos axissimétricos resolvidos pelo MEC.

Este procedimento consiste basicamente de uma integração em torno do eixo

de axissimetria para θ variando de 0 a 2π, de acordo com a figura 1.1. Em (1.16), a

divisão por 2π é meramente um artifício que terá sua utilidade mostrada ao final do

capítulo.

∫=

π

θπ

2

0

*3

*

21 dpp daxi (1.16)

Para substituir a solução fundamental tridimensional, dada por (1.11), em

(1.16), é necessário expressar a distância entre um ponto qualquer do anel de fontes

descrito por ξ (figura 1.1) e o ponto campo X em função apenas de variáveis presentes

no plano 1-2, que contém o contorno do modelo, e do ângulo θ, que será integrado,

desaparecendo da expressão. Deste modo, chega-se por simples geometria a (1.17).

( )222

2111

21 cos2 ξξθξ −++−= xxxR (1.17)

Figura 1.1 – variáveis necessárias para solução axissimétrica


Escrevendo-se r conforme (1.18), R passa a ser dado pela expressão (1.19).

( )222

2111

21 2 ξξξ −++−= xxxr (1.18)

( ) ( )θξθ cos12, 112 −+= xrrR (1.19)

A distância máxima d entre um ponto do anel de fontes e o ponto campo em

questão, dada por (1.20) é obtida para cos θ = – 1.

112 4 ξxrd += (1.20)

Substituindo-se a expressão (1.19) na solução fundamental dada por (1.11) e

esta em (1.16), pode-se proceder à integração indicada. A matemática desta integração

é um tanto árdua, não sendo aqui apresentada. Ela pode ser encontrada na referência

[24].

A expressão final a que se chega é apresentada em (1.21).

( )[ ] ( )( )( )[ ]

( )[ ]τξ

ττξτπ

−−+

−−

⎥⎦⎤

⎢⎣⎡ −−−−−

=

tcrxH

rtcHrtcxrtc

cpaxi

211

22211

222

*

4

414

x

x

(1.21)

A expressão (1.21) é não-nula somente para r < c(t – τ) < 2114 rx +ξ = d, ou

seja, apenas para pontos já atingidos pela onda emitida pelo ponto do anel de fontes

mais próximo a eles e ainda não atingidos pelo ponto do anel de fontes mais distante

deles. A interpretação deste fato conduz ao raciocínio de que o caso axissimétrico é

um caso intermediário entre o caso bidimensional e o tridimensional. No primeiro, o

ponto campo é afetado pelo ponto fonte desde a chegada da primeira onda até o

infinito, já que este ponto fonte na realidade representa uma linha infinita de fontes.

No caso tridimensional, o ponto campo é afetado apenas em um instante de tempo

determinado, já o ponto fonte representa apenas ele mesmo, e a solução fundamental é

um Delta de Dirac. Já para o caso axissimétrico, a solução fundamental é não-nula, ou

Formulação do MEC para problemas acústicos axissimétricos 15 seja, o ponto campo é afetado, durante um intervalo de tempo definido, compatível

com o fato de o ponto fonte na realidade representar uma circunferência formada por

fontes pontuais infinitesimais.

A derivada da solução fundamental axissimétrica na direção normal ao

contorno é dada por (1.22), onde o vetor n representa a normal ao contorno, de

componentes n1 e n2, sendo x1 e x2 as coordenadas do ponto campo, sempre em

relação aos eixos 1 e 2 da figura 1.1.

( )[ ] ( )( )

( )[ ] ( )( )( )[ ] ( )[ ]τξτ

τξτ

ξτξξ

π

−−+−−

⎭⎬⎫

⎩⎨⎧

⎥⎦⎤

⎢⎣⎡ −−−−−

−⎥⎦⎤

⎢⎣⎡ −−−−+

=

=∂∂

+∂

∂=

∂∂

=

tcrxHrtcH

rtcxrtc

xnrtcxxnxnc

nxp

nx

pn

pq axiaxiaxi

axi

211

2/3222

11222

2111

2221122211

2

22

*

11

***

4

21.

21.

4

x

x (1.22)

Nota-se que as funções Heaviside originárias de (1.21) não foram derivadas.

Fazendo-se analogia com o caso bidimensional, a derivação da função Heaviside em

relação ao tempo acaba produzindo um termo que inclui as condições iniciais e outro

que se soma à expressão (1.22) na integral ao longo do tempo, regularizando o

integrando. Deste modo, os problemas analisados utilizando-se o fluxo da solução

fundamental dado por (1.22) apresentam erro quando da consideração de condições

iniciais no problema em seu contorno.

O resultado da integral, com o integrando regularizado, equivale ao obtido

integrando-se diretamente a expressão (1.22), utilizando o conceito de partes finitas

quando há singularidades não-integráveis. No entanto, na formulação por partes

finitas, falta o termo que inclui as condições iniciais. No Capítulo 2, o procedimento

de integração utilizado neste trabalho é apresentado, e este aspecto volta a ser

discutido.

No caso em que o ponto fonte ξ ou o ponto campo X se localiza sobre o eixo

de axissimetria do modelo, a solução fundamental axissimétrica confunde-se com a

solução fundamental tridimensional. Esse aspecto só é possível se houver a divisão

por 2π presente em (1.16).

Implementação numérica do MEC 16

Capítulo 2

Implementação numérica do MEC

Neste capítulo, são abordados todos os aspectos relativos ao desenvolvimento

do esquema elaborado neste trabalho com o objetivo de se resolver problemas

acústicos axissimétricos, representados pela equação (1.1), através da implementação

numérica pelo Método dos Elementos de Contorno (MEC).

Primeiramente, mostra-se de que forma, partindo da equação integral (1.13),

pode-se resolver o problema numericamente. Este processo envolve a discretização do

contorno em elementos, bem como a discretização do tempo total de análise em

instantes de tempo definidos a partir de um intervalo de tempo de análise,

procedimento padrão para qualquer problema a ser resolvido no domínio do tempo

pelo MEC.

Em seguida, são apresentadas as aproximações utilizadas a fim de se

representar o comportamento da variável básica do problema e de sua derivada na

direção normal ao contorno (no caso da acústica, a pressão e o fluxo) ao longo do

espaço e do tempo em função de valores discretos nos nós do contorno (nós estes que

são determinados pela discretização espacial adotada) em certos instantes de tempo

(determinados em função da discretização temporal adotada), valores estes que

acabam se constituindo nas incógnitas do problema.

Tendo sido definidas as aproximações, fica definido o integrando das integrais

de convolução presentes na expressão (1.13). Parte-se então para a manipulação deste

integrando, a fim de se desenvolver a expressão analítica do resultado desta integral,

que posteriormente é integrado ao longo dos elementos, conforme será visto.

Estas expressões analíticas apresentam singularidades para alguns pontos do

contorno, sobre as quais é comentado na seqüência, apresentando-se o procedimento

de integração semi-analítica utilizado, e mostrando como é realizada a integração

numérica ao longo dos elementos do contorno, o que envolve uma transformação de

coordenadas a fim de se aumentar a precisão do método de integração. Desta forma,

são montadas as matrizes de convolução, bem como as matrizes correspondentes ao

tempo para o qual as incógnitas estão sendo calculadas.

Implementação numérica do MEC 17 Posteriormente são apresentadas as expressões analíticas integradas ao longo

do tempo relativas aos termos da equação integral que aparecem devido a fontes

pontuais geradoras de pressão presentes no meio acústico.

Neste trabalho, a integração analítica ao longo do tempo apresenta-se como

estratégia alternativa e mais precisa à integração numérica, através dos métodos de

Kutt [27] e Chebyshev [28], de acordo com o tipo de integral a ser calculada.

Por fim, aborda-se a implementação do “Esquema Teta”, que foi empregado

neste trabalho tendo como objetivo a estabilização da resposta ao longo do avanço no

tempo para alguns problemas com os quais se deparou.

2.1. Discretização e aproximações

No Capítulo 1, obteve-se a equação (1.13), que serve de base para a resolução

do problema acústico a partir apenas de incógnitas presentes no contorno do meio em

estudo. Por aquela equação, pode-se obter a solução em qualquer ponto, em qualquer

instante de tempo, se for conhecida a solução ao longo de todo o contorno do meio em

questão, em todo instante de tempo anterior a este.

O MEC tem como princípio a obtenção da solução em pontos definidos, no

contorno. Para isso, a equação (1.13) deve ser aplicada a estes pontos, de modo que

eles funcionem como pontos fontes. Estes pontos, chamados de nós funcionais, são

neste trabalho sempre posicionados nas extremidades dos elementos através dos quais

a geometria do contorno é discretizada. Os elementos com os quais se trabalha aqui

têm geometria sempre reta. Uma discretização espacial típica através de elementos

lineares é apresentada na figura 2.1.

(a) (b)

Figura 2.1 – (a) contorno original; (b) contorno discretizado por elementos lineares

Implementação numérica do MEC 18 De modo análogo, pode-se proceder com relação à discretização temporal.

Neste trabalho, adota-se sempre o valor do intervalo de tempo constante, ou seja,

define-se um intervalo de tempo, sendo o problema então resolvido para instantes de

tempo múltiplos do valor deste intervalo.

A equação integral (1.13) é aplicada considerando-se, de cada vez, um nó

funcional diferente do contorno como ponto fonte. Introduzindo-se as discretizações

empregadas, chega-se a um sistema de equações algébricas lineares, já que a cada

aplicação de (1.13) obtém-se uma equação que contém os valores da variável básica e

de sua derivada em todos os nós do contorno. A resolução do problema passa então

pela montagem e solução deste sistema de equações.

2.1.1. Funções de interpolação

Para as integrais presentes em (1.13) serem efetuadas, é necessário se definir

os valores da pressão p e do fluxo q em qualquer ponto, em qualquer instante de

tempo, em função dos valores nos nós da discretização espacial e temporal, o que é

feito através de funções de interpolação, conforme mostrado abaixo.

( ) ∑∑=

+

=

=NN

j

n

m

mj

jm pXttXp1

1

1

)()(, ηφ (2.1)

( ) ∑∑=

+

=

=NN

j

n

m

mj

jm qXttXq1

1

1

)()(, νθ (2.2)

Nas equações (2.1) e (2.2), X corresponde a um ponto qualquer do contorno e t

a um instante qualquer de tempo; NN é o número total de nós e n o número total de

intervalos de tempo decorridos até aquele instante; o índice m percorre os instantes de

tempo discretizados até o instante para o qual se está resolvendo o problema,

enquanto o índice j percorre todos os nós do contorno. pjm e qj

m representam,

respectivamente, os valores da pressão e do fluxo no nó j no tempo m.

As funções de interpolação φ m(t), θ m(t), η j(X) e ν j(X), como qualquer função

de interpolação, devem ter como requisito básico assumir valor unitário no nó ao qual

se referem e nulo em todos os demais. Neste trabalho, são empregadas funções de

interpolação sempre lineares no espaço. Para a interpolação ao longo do tempo, são

utilizadas funções lineares, para p, e constantes, para q.


Deste modo, as funções de interpolação podem ser representadas a partir das

expressões (2.3) a (2.5).

Funções de interpolação no tempo:

⎪⎪⎪

⎩

⎪⎪⎪

⎨

⎧

><

≤≤Δ

−

≤≤Δ

−

=

+−

++

−−

11

11

11

,0

,

,

)(

mm

mmm

mmm

m

tout

ttt

t

tttt

ττ

ττ

ττ

τφ (2.3)

⎩⎨⎧

><≤≤

=−

−

mm

mmm

touttt

τττ

τθ1

1

,0,1

)( (2.4)

Funções de interpolação no espaço:

( )

( )⎪⎪⎩

⎪⎪⎨

⎧

Γ∈−

Γ∈+==

ll

kkjj

X

XXX

,121

,121

)()(ξ

ξνη (2.5)

Nas expressões (2.5), Γk e Γl representam os elementos adjacentes ao nó j, e ξk

e ξl suas respectivas coordenadas naturais, conforme figura 2.3.

Conforme já mencionado, as funções φ, θ e η assumem comportamento linear

entre pontos adjacentes, enquanto a função ν assume um comportamento constante, o

que ilustram as figuras 2.2 e 2.3.

Figura 2.2 – funções de interpolação no tempo φ

m(t) e θ m(t)


Figura 2.3 – funções de interpolação no espaço ηj(X) e νj(X) ao longo dos elementos adjacentes ao nó j

Como as funções são sempre lineares ou constantes, elas assumem valores não

nulos apenas em elementos ou intervalos de tempo adjacentes ao ponto em relação ao

qual se referem. Portanto, os valores de p e q em um ponto qualquer em certo instante

de tempo são determinados em função apenas dos valores nos nós extremos do

elemento do qual este ponto faz parte, nos instantes de tempo discretizados

imediatamente posterior e anterior ao instante de tempo em questão.

2.1.2. Montagem do sistema de equações

Tendo as equações (2.1) e (2.2) totalmente definidas, substituem-se na

equação (1.13) os valores de p e q em função das aproximações, chegando-se à

expressão (2.6), em que se considera o ponto-fonte localizado no nó i. Não são

consideradas condições iniciais de pressão nem de sua derivada em relação ao tempo

em nenhum ponto do domínio.

( ) )1(1

1 1

)1(1

1 1

)1(1 ++

= =

++

= =

++ +=+ ∑∑∑∑ ni

n

m

NN

j

mj

mnij

n

m

NN

j

mj

mnij

nii SqGpHpc ξ (2.6)

Na equação (2.6):

∫ ∫Γ

++ Γ= )()(),;,()()(

0

1*

1)1( XddtXqXxXH axi

t

min

jmnij

n

ττφτξη (2.7)

∫ ∫Γ

++ Γ= )()(),;,()()(

0

1*

1)1( XddtXpXxXG axi

t

min

jmnij

n

ττθτξν (2.8)

τττξ dftXpSnt

sinS

NF

s

ni ∫∑ +

=

+ =0

1*

1

)1( )(),;,( (2.9)


Nas equações (2.7) a (2.9), p* e q* representam, respectivamente, a solução

fundamental axissimétrica para pressão e sua derivada na direção normal ao contorno

para um ponto X, no instante tn+1, devidas a uma fonte localizada em ξi no instante τ.

A função fs(τ) representa a variação ao longo do tempo da vibração emitida pela s-

ésima fonte pontual presente no domínio, sendo NF o número total de fontes. O

coeficiente c(ξi) é aquele descrito no Capítulo 1, dividido por 2π, já que a solução

fundamental axissimétrica, conforme já citado, também sofre esta divisão.

Aplicando-se a equação (2.6) aos NN nós do contorno, obtêm-se n+1 matrizes

G e H, de dimensão NNxNN, bem como um vetor S, de dimensão NN. Devido às

propriedades das funções de interpolação citadas anteriormente, e ilustradas nas

figuras 2.2 e 2.3, para se obter as matrizes H(n+1)m e G(n+1)m basta se proceder à

integração entre os instantes de tempo tm-1 e tm+1. Para se obter os elementos Hij e Gij

de cada uma dessas matrizes, basta realizar a integração ao longo dos elementos

adjacentes ao nó j.

A partir da equação (2.6), pode-se montar um sistema de equações através do

qual se obtém a solução para o instante de tempo n+1, chamado de instante “atual” de

tempo. Neste instante, já são conhecidos os valores dos vetores pm e qm para m <

n+1, que multiplicados, respectivamente, pelas matrizes H(n+1)m e G(n+1)m, realizam a

convolução do problema. Portanto, as incógnitas são pn+1 e qn+1, ou seja, 2NN valores

desconhecidos. Como o sistema apresenta NN equações, seria indeterminado, não

fossem as condições de contorno, que impõem que em todos os instantes, para cada

nó, ou a pressão ou o fluxo deve ser prescrito, gerando então NN valores

desconhecidos de pressão e fluxo.

Para a resolução do sistema, é formada uma matriz que multiplica somente as

incógnitas. Primeiramente, soma-se à diagonal da matriz H(n+1)(n+1), em cada linha i, a

parcela c(ξi), formando-se a matriz H(n+1)(n+1). Escolhendo-se, por exemplo, a matriz

H(n+1)(n+1) como base para a formação da matriz de incógnitas, são trocadas, mudando

o sinal, suas colunas correspondentes aos nós com p prescrito pelas correspondentes

colunas da matriz G(n+1)(n+1), já que nestes nós a incógnita é q. Estas novas matrizes

derivadas de H(n+1)(n+1) e G(n+1)(n+1) são então denominadas, respectivamente, A e B,

formando o sistema de equações indicado em (2.10).

n

n

m

NN

j

mj

mnij

n

m

NN

j

mj

mnij pHqG SByAx +−+= ∑∑∑∑

= =

+

= =

+

1 1

)1(

1 1

)1( (2.10)


Na equação (2.10), x contém todas as incógnitas no instante n+1 e y, os

valores prescritos neste mesmo instante.

2.1.3. Propriedade de translação

As matrizes H e G apresentam uma propriedade bastante útil, chamada

propriedade de translação. Esta propriedade pode ser representada através das

seguintes expressões:

nmlmln

nmlmln

GGHH

=

=++

++

))((

))((

(2.11)

As expressões (2.11) indicam que as matrizes H e G dependem apenas da

diferença entre o “instante atual” n (para o qual se está efetuando o cálculo) e o

instante m no qual a solução fundamental foi emitida, sendo l um inteiro genérico. A

partir de agora, portanto, elas podem passar a ser identificadas apenas por um índice,

que representa a diferença (n – m). A figura 2.4 elucida bem esta propriedade. Nela,

no eixo vertical estão representadas a função de interpolação e a solução fundamental.

Pode-se perceber que, em ambos os casos, a integral do produto destas duas funções,

presente na expressão (2.7), é a mesma. Portanto: Hnm = H(n-l)(m-l) = H(n-m).

Figura 2.4 – exemplificação da propriedade de translação


Pode-se tirar bastante proveito desta propriedade, visto que, a cada passo de

tempo, podem-se utilizar resultados obtidos em instantes de tempo anteriores,

bastando calcular parcelas das matrizes G e H relativas a funções de interpolações

não-nulas entre τ = t1 e τ = t2, ou seja, relativas a G(n-1), H(n) e H(n-1), sendo n+1 o

instante de tempo atual, correspondente ao passo de tempo n. A figura 2.5 ilustra o

que foi dito. No passo de tempo seguinte, soma-se mais uma parcela aos termos de

H(n). A parcela equivalente àquela somada no instante atual a H(n) foi somada no

instante anterior a H(n-1).

Aqui se verifica a restrição imposta pela não derivação das funções Heaviside

da solução fundamental axissimétrica, o que regularizaria a expressão a ser integrada.

Apenas ao se somarem as duas parcelas de H(n) calculadas em instantes de tempo

consecutivos chega-se aos valores corretos de seus termos. Deste modo, sendo o

instante de tempo atual n+1, a matriz H(n), que multiplica p1 (pressões iniciais), ainda

não está completa, portanto apresenta valores incorretos, inviabilizando a análise com

condições iniciais não-nulas no contorno. Comportamento análogo a este é

apresentado e detalhado, para o caso bidimensional, em MANSUR e CARRER [29].

Figura 2.5 – região onde é necessário efetuar as integrais de convolução em cada passo de tempo, com

suas respectivas funções de interpolação


2.2. Integração analítica no tempo

Neste item, são apresentados os procedimentos empregados para o cálculo das

integrais ao longo do tempo (integrais de convolução) presentes nas equações (2.7) e

(2.8). No presente trabalho, as integrais são desenvolvidas analiticamente, gerando

expressões que posteriormente são integradas ao longo dos elementos de contorno

para se chegar aos elementos das matrizes G e H.

Conforme mencionado no item 2.1.3, as integrais, a cada passo de tempo, são

executadas sempre entre τ = t1 e τ = t2. Para a matriz G, como a função de

interpolação no tempo de q é constante, calcula-se apenas uma expressão, chamada de

g, que após ser integrada no espaço, produz os elementos da matriz G(n-1), sendo n+1

o instante atual de tempo. Já para a matriz H, como a função de interpolação de p é

linear, são calculadas duas expressões, relativas às duas funções de interpolação não-

nulas neste intervalo. A expressão relativa a φ1(τ) é chamada de hI, enquanto a

expressão relativa a φ2(τ) é chamada de hF. Para tornar as expressões genéricas para

qualquer intervalo de tempo, os limites do intervalo no qual se processa a integração

são tratados, a partir de agora, como ti e tf, ao invés de t1 e t2.

2.2.1. Ponto fonte fora do eixo de axissimetria

Quando o ponto fonte se encontra fora do eixo de axissimetria, a solução

fundamental e sua derivada em relação à normal ao contorno são regidas pelas

expressões (1.21) e (1.22).

Primeiramente, apresenta-se aqui a mudança de variável de integração

realizada a fim de se facilitar a manipulação das expressões a serem integradas.

Originalmente, têm-se as expressões (2.12) a (2.14).

( )

( )[ ] ( )( )τ

ττθτξξ

drτtcξxrτtcπ

c

dtXptXg

f

i

f

i

t

t nn

t

t

finn

∫

∫

⎥⎦⎤

⎢⎣⎡ −−−−−

==

++

++

221

211

221

22

1*

1

414

)(),;,(,,

(2.12)


( )

( )[ ] ( )( )

( )[ ] ( )( )τ

τξξ

ττφτξξ

dt

t

rτtcξxrτtc

xnrτtcξxxnxn

πc

dtXqtXh

f

i

f

i

t

t

f

nn

n

t

t

iinnI

∫

∫

⎟⎟⎠

⎞⎜⎜⎝

⎛Δ

−

⎭⎬⎫

⎩⎨⎧

⎥⎦⎤

⎢⎣⎡ −−−−−

−⎥⎦⎤

⎢⎣⎡ −−−−+

=

==

++

+

++

2/322

12

1122

12

2111

221

21122211

2

1*

1

41

21

4

)(),;,(,,

(2.13)

( )

( )[ ] ( )( )

( )[ ] ( )( )τ

τξξ

ττφτξξ

dtt

rτtcξxrτtc

xnrτtcξxxnxn

πc

dtXqtXh

f

i

f

i

t

t

i

nn

n

t

t

finnF

∫

∫

⎟⎠

⎞⎜⎝

⎛Δ−

⎭⎬⎫

⎩⎨⎧

⎥⎦⎤

⎢⎣⎡ −−−−−

−⎥⎦⎤

⎢⎣⎡ −−−−+

=

==

++

+

++

2/322

12

1122

12

2111

221

21122211

2

1*

1

41

21

4

)(),;,(,,

(2.14)

Como se pode notar, não são colocadas no integrando as funções Heaviside

que multiplicam a solução fundamental e sua derivada nas equações (1.21) e (1.22).

Conforme será exposto mais adiante neste item, isto é compensado pela alteração nos

limites das integrais presentes nas expressões (2.12) a (2.14), de acordo com o caso no

qual se enquadra o ponto X do contorno a que se refere a integral.

Define-se, em seguida:

112 4 ξxrd += (2.15)

Em (2.15), d é a máxima distância entre o anel de fontes representado pelo

ponto fonte do contorno ξ e o ponto campo X, conforme mostrado no Capítulo 1. A

mudança de variável é então realizada:

( )τ−= +1ntcu (2.16)

Com as definições (2.15) e (2.16), as expressões (2.12) a (2.14) podem ser

reescritas da seguinte forma, mais útil com vistas à integração:

=

=

=

=


( )( )( )

duudruπ

tXgi

f

u

u

n ∫ −−=

2222221,,ξ (2.17)

( ) ( )[ ][ ]( )( )[ ] du

tcπudru

uuruCCtXh

i

f

u

u

fnI ∫ Δ−−

+−−−= 22/32222

2221 12

,,ξ (2.18)

( ) ( )[ ][ ]( )( )[ ] du

tcπudru

uuruCCtXh

i

f

u

u

inF ∫ Δ−−

−−−= 22/32222

2221 12

,,ξ (2.19)

Nas expressões (2.17) a (2.19):

( )( )fnf

ini

ttcuttcu

−=−=

(2.20)

( ) ( )[ ] 11112221111 ξξξξ rxnrxxnxnC

∂∂

=−+−= (2.21)

( )[ ]222112 ξ−+= xnxnC (2.22)

Conforme já foi dito, os limites de integração nem sempre são aqueles

apresentados nas expressões (2.17) a (2.19). Como foi mostrado no Capítulo 1, um

ponto campo é afetado por determinado ponto fonte apenas no intervalo de tempo

definido desde a chegada da excitação provocada por este ponto até a passagem da

excitação provocada pelo ponto pertencente ao anel de fontes mais distante dele, o

que é matematicamente descrito pelas funções Heaviside em (1.21) e (1.22). Portanto,

a solução fundamental p*, bem como sua derivada na direção normal ao contorno q*,

ficam definidas apenas neste intervalo, representado matematicamente na expressão

(2.23), na qual ele é transformado, sendo expresso em função da nova variável de

integração.

durcdt

crt nn ≤≤⇒−≤≤− ++ 11 τ (2.23)

=

=

Implementação numérica do MEC 27 Há a possibilidade de existência de seis casos diferentes de limites de

integração, de acordo com as distâncias entre o ponto campo e os pontos mais

próximo e mais distante a ele do anel de fontes representado pelo ponto fonte original

ξ, para cada qual sendo definido um intervalo de integração para as integrais presentes

nas expressões (2.17) a (2.19), conforme apresentado em [24] e ilustrado na figura

2.6:

− Caso 1: r > c(tn+1 – ti) = ui

− Caso 2: r > c(tn+1 – tf) = uf e d > c(tn+1 – ti) = ui ≥ r

− Caso 3: r > c(tn+1 – tf) = uf e c(tn+1 – ti) = ui ≥ d

− Caso 4: c(tn+1 – tf) = uf > r e d > c(tn+1 – ti) = ui

− Caso 5: d > c(tn+1 – tf) = uf ≥ r e c(tn+1 – ti) = ui ≥ d

− Caso 6: c(tn+1 – tf) = uf > d

ui

ui

ui

ui

ui

ui

uf

uf

uf

uf

uf

uf

r

r

r

r

d

dr

d

d

d

CASO 1

CASO 2

CASO 4

CASO 3

CASO 5

CASO 6

Figura 2.6 – representação dos seis casos de limites de integração

Implementação numérica do MEC 28 Os limites de integração correspondem à interseção entre o intervalo definido

por r e d e aquele definido por uf e ui (limites originais) na figura 2.6.

Para os casos 1 e 6, pode-se afirmar de imediato que as integrais têm resultado

nulo, visto que a interseção é nula. No caso 1, ainda não chegou a excitação

provocada pelo ponto mais próximo do anel de fontes (ponto ξ, à distância r do ponto

campo), enquanto no caso 6 já passou a influência provocada pelo ponto mais distante

do anel de fontes (à distância d do ponto-campo). Para os demais casos, é necessária a

integração, ao longo dos seguintes intervalos:

− Caso 2: r até ui

− Caso 3: r até d

− Caso 4: uf até ui

− Caso 5: uf até d

As integrais para montagem da matriz G, obtidas pela equação (2.17), são de

obtenção imediata, empregando-se o conceito de integral elíptica [28], sendo seus

resultados apresentados no Anexo A. No Anexo B, encontram-se maiores explicações

acerca do conceito de integral elíptica.

Já as integrais presentes nas equações (2.18) e (2.19), responsáveis pela

montagem da matriz H, têm desenvolvimento um pouco mais complexo, mesmo

porque, dependendo dos limites de integração, elas podem ser compreendidas

somente através do conceito de parte finita de integral [26], caso haja singularidades

não-integráveis. Maiores detalhes acerca deste conceito são fornecidos no Anexo C.

Para desenvolvê-las, o primeiro passo foi a separação do integrando em quatro

parcelas distintas (toma-se para exemplificar o procedimento, a partir de agora,

apenas a integral correspondente a hI, sendo aquela correspondente a hF análoga). Em

seguida, são apresentados as quatro parcelas e os procedimentos adotados para a

integração de cada uma delas.

No Anexo A, encontram-se os resultados da integral para cada uma das

parcelas separadamente, para os quatro casos (2 a 5).


Parcela 1:

( )( )[ ] tcπudru

uC f

Δ−−22/32222

1 12 (2.24)

Pode-se notar que a expressão (2.24) compõe-se de uma constante dividida

pelo fator que está elevado a 3/2 no denominador. Na realidade, então, o que se deseja

integrar é:

( )( )[ ] 2/32222

1udru −−

(2.25)

− Caso 2:

O integrando (2.25) apresenta uma singularidade não-integrável no limite

inferior de integração (r). Adota-se então o seguinte procedimento:

( )( )[ ] ( ) ( )du

rurAdu

rurAuA

udrudu

iii u

r

u

r

u

r∫∫∫ −

+−

−=

−−2/32/32/32222

)()()( (2.26)

Em (2.26):

( )( )[ ] 2/322

1)(ruud

uA+−

= (2.27)

A primeira parcela do segundo membro de (2.26) pode ser integrada

analiticamente de modo direto, enquanto para a segunda parcela utiliza-se o resultado

apresentado em (2.28) [26] e demonstrado no Anexo C, que corresponde à parte finita

da integral.

( ) 2/12/3 )()(2)(

rurAdu

rurA

i

u

r

i

−−

=−∫ (2.28)

= =

=


− Caso 3:

Neste caso, o integrando apresenta singularidades não-integráveis em ambos

os limites de integração. Porém, podem-se aqui aproveitar resultados de outros casos

pela utilização da seguinte propriedade:

∫∫∫ +=

d

u

u

r

d

r f

f

duufduufduuf )()()( (2.29)

Deste modo, basta somar as duas parcelas seguintes: aquela equivalente à do

caso 2, trocando-se ui por uf ; e aquela correspondente ao caso 5, apresentado mais

adiante.

− Caso 4:

Na realidade, esta integral pode ser desenvolvida normalmente, visto que o

integrando não apresenta singularidades no intervalo de integração. Porém, para se

aproveitar do resultado obtido para o caso 2, leva-se em conta a seguinte propriedade:

∫∫∫ −=

fii

f

u

r

u

r

u

u

duufduufduuf )()()( (2.30)

Portanto, toma-se o resultado obtido para o caso 2 e subtrai-se dele expressão

equivalente na qual se troca ui por uf.

− Caso 5:

Novamente, o integrando (2.25) apresenta uma singularidade não-integrável,

agora no limite superior de integração (d). Adota-se então o seguinte procedimento:

( )( )[ ] ( ) ( )du

uddBdu

udrBuB

udrudu

d

u

d

u

d

u fff

∫∫∫ −+

−−

=−−

2/32/32/32222

)()()( (2.31)

= =


Em (2.31):

( )( )[ ] 2/322

1)(udru

uB+−

= (2.32)

Assim como para o Caso 2, a primeira parcela do segundo membro de (2.31)

pode ser integrada analiticamente de modo direto, enquanto para a segunda parcela

utiliza-se o resultado apresentado em (2.33) [26], que corresponde à parte finita da

integral.

( ) 2/12/3 )()(2)(

f

d

uud

dBduuddB

f

−−

=−∫ (2.33)

Parcela 2:

( )( )[ ] tcπudru

uCΔ−−

22/32222

1 12 (2.34)

Pode-se notar que a expressão (2.34) compõe-se de uma constante

multiplicada por u e dividida pelo fator que está elevado a 3/2 no denominador. Na

realidade, então, o que se deseja integrar é:

( )( )[ ] 2/32222 udruu

−− (2.35)

Os procedimentos para todos os casos são análogos aos da Parcela 1, sendo

que as expressões A(u) e B(u) são agora multiplicadas por u.

Parcela 3:

( )( )( )[ ] tcudru

ruuC f

Δ−−

−−22/32222

222 1

π (2.36)

=


Em (2.36), há uma constante multiplicada por:

( ) ( ) 2/3222/122

1udru −−

(2.37)

− Caso 2:

Para esta parcela a singularidade apresentada pelo integrando (2.37) no limite

inferior r é integrável, visto que o fator que contém esta singularidade é elevado a ½.

Portanto, a integração é implementada analiticamente de modo usual.

− Casos 3 e 4:

Repete-se o procedimento adotado para as parcelas anteriores.

− Caso 5:

Da mesma forma, repete-se o procedimento para as parcelas anteriores,

fazendo-se uma ressalva sobre a mudança na expressão de B(u), apresentada em

(2.38).

( ) ( ) 2/32/122

1)(udru

uB+−

= (2.38)

Parcela 4:

( )( )( )[ ] tcudru

ruuCΔ−−

−−22/32222

222 1

π (2.39)

A expressão (2.39) corresponde a uma constante multiplicada por:

( ) ( ) 2/3222/122 udruu

−− (2.40)

Os procedimentos são os mesmos adotados para a Parcela 3, observando-se

que a expressão B(u) é multiplicada por u.


É importante ressaltar que para a integração numérica das expressões obtidas

desta forma e apresentadas no Anexo A ao longo dos elementos de contorno produzir

resultados mais precisos, é necessário se dividir os elementos em segmentos, sendo

que cada segmento passa a conter apenas pontos pertencentes ao mesmo caso. Isto é

facilmente observado em vista da diferença entre as expressões obtidas para cada

caso.

Para esta divisão procuram-se as interseções de quatro círculos, de raios ui e uf,

com centros no ponto-fonte ξ e no ponto mais distante do anel relativo a ele (trata-se,

na realidade, do ponto simétrico a ξ, no plano bidimensional que contém Γaxi, em

relação ao eixo de axissimetria), com o elemento sobre o qual se está integrando.

Cada segmento, determinado por dois pontos de interseção destes, ou por

extremidades do elemento, corresponde a um caso diferente.

2.2.2. Ponto fonte pertencente ao eixo de axissimetria

Conforme exposto no Capítulo 1, na eventualidade do ponto fonte se localizar

sobre o eixo de axissimetria, a solução fundamental e sua derivada na direção normal

se transformam naquelas obtidas para o caso tridimensional, já que ao invés de um

anel de fontes, este ponto representará apenas ele mesmo (como se fosse um anel de

raio nulo). Portanto, a excitação por ele emitida em determinado instante de tempo se

comporta como um delta de Dirac, assim como no caso tridimensional, afetando cada

ponto campo apenas em um determinado instante de tempo (ou por um intervalo de

tempo infinitesimal).

Introduzindo-se as expressões (1.11) e (1.12) nas integrais no tempo presentes

em (2.7) e (2.8), chega-se às expressões (2.41) e (2.42), nas quais já são efetuadas as

mudanças de variáveis apresentadas em (2.20).

( ) ( )[ ] ( ) [ ]duruπr

drtcπrctXg

i

f

f

i

u

u

t

t

fnn ∫∫ −=−−= δττθτδξ

41

4,, (2.41)


( )

( ) ( )∫∫

∫

−∂∂

+−−

∂∂

=

=⎟⎟⎠

⎞⎜⎜⎝

⎛⎟⎠⎞

⎜⎝⎛ −−+⎟

⎠⎞

⎜⎝⎛ −−−

∂∂

=

i

f

i

f

f

i

u

u

FI

u

u

FI

t

t

FInnnFI

duuruπrcn

rduuruπrn

r

dcrt

ccrt

rπrnrtXh

)(4

1)(4

1

)(1141,,

//2

//

φδφδ

ττφτδτδξ

&

&

(2.42)

O resultado final obtido efetuando-se as integrais em (2.41) e (2.42), com a

substituição das funções de interpolação presentes em (2.42), está apresentado nas

expressões (2.43) a (2.45).

( )πr

tXg n 41,, =ξ (2.43)

( )tπrcn

rtcur

πrnrtXh f

nI Δ∂∂

+Δ

−

∂∂

−=4

14

1,, 2ξ (2.44)

( )tπrcn

rtcru

πrnrtXh i

nF Δ∂∂

−Δ−

∂∂

−=4

14

1,, 2ξ (2.45)

Para esta situação, não valem aqueles seis casos apresentados no item 2.2.1,

devido à natureza da solução fundamental (r = d) para qualquer ponto campo X. São

definidos apenas três casos:

− Caso 1b: r > ui

− Caso 2b: ui ≥ r ≥ uf

− Caso 3b: uf > r

No caso 1b, a excitação emitida pela fonte em ti ainda não atingiu o ponto

campo à distância r; no caso 3b, a excitação emitida pela fonte em tf já passou pelo

ponto campo à distância r; o caso 2b apresenta uma situação intermediária, na qual o

ponto campo está recebendo uma excitação proveniente do ponto campo, emitida

entre ti e tf, sendo o único caso em que é necessária a integração. Portanto, cada

elemento sobre o qual se está realizando a integração no espaço deve ser dividido em

Implementação numérica do MEC 35 segmentos determinados pela interseção entre ele mesmo e dois círculos de raios ui e

uf com centro no ponto fonte ξ, sendo a integração efetuada somente no segmento

cujos pontos pertencem ao caso 2b.

2.3. Integração no espaço e singularidades

A fim de se realizar a montagem das matrizes G e H, as expressões obtidas

para as integrais de convolução (Anexo A e (2.43) a (2.45)) devem ser agora

introduzidas em (2.7) e (2.8). Executam-se então as integrais ao longo do contorno,

elemento por elemento.

Essas integrais são, neste trabalho, efetuadas numericamente através do

processo usual de Gauss-Legendre. No entanto, observa-se que o integrando pode

apresentar singularidades em alguns pontos, pela presença de limites tendendo a

infinito. Neste caso, apesar de a singularidade ser integrável, a integração numérica se

mostra menos precisa, de modo que estas singularidades devem ser tratadas de modo

especial, pela preservação da qualidade do resultado da integral.

Nos itens 2.3.1, 2.3.2 e 2.3.3, são apresentados os três tipos de singularidades

passíveis de ocorrência nas expressões do Anexo A e nas expressões (2.43) a (2.45),

bem como o tratamento dado a cada um deles. Utiliza-se o aqui chamado de

“processo semi-analítico”, no qual o integrando, originalmente de difícil integração

analítica, devido à complexidade da sua expressão, é desmembrado em uma parcela

mais simples, a ser integrada analiticamente e outra que não apresenta singularidade

do tipo limite infinito, a ser integrada numericamente. Isto é feito subtraindo-se e

somando-se o valor do fator não-singular do integrando no ponto onde há a

singularidade, multiplicado pelo fator singular.

A parcela não-singular apresenta em geral, como será visto, um forte gradiente

nas proximidades da antiga singularidade, sendo mais adequado o emprego de uma

transformada apresentada por TELLES [30], resumidamente tratada no item 2.3.4,

cujo objetivo é a maior concentração de pontos de integração próximos à região de

alto gradiente.


2.3.1. Singularidades em r = 0

As singularidades discutidas neste item aparecem ao se integrar ao longo do

elemento ao qual pertence o ponto fonte, tendo em vista que neste caso o integrando

inclui o ponto em que r = 0. O primeiro caso reportado é o das integrais elípticas do

primeiro tipo presentes nas expressões do Anexo A. Quando r = 0, pode-se inferir

imediatamente que ambos os argumentos da integral elíptica se igualam a 1. Neste

caso, a integral tende a infinito. É apresentado em seguida o procedimento empregado

para as expressões relativas a g (item A.1):

1) Primeiramente, subtrai-se e soma-se, da integral elíptica incompleta,

quando for o caso, a integral elíptica completa. Sabe-se que a diferença entre estas

duas integrais possui limite finito em r = 0, como se pode inferir da Figura 2.7. Este

procedimento é exposto em (2.46), onde são consideradas integrais elípticas presentes

nas expressões relativas a g (item A.1). Para simplificar a notação e facilitar a

visualização dos desenvolvimentos, os argumentos da integral elíptica são

simbolizados por a e b.

( ) ( ) ( ) ( )∫∫∫ΓΓΓ

−+

−⎥⎦⎤

⎢⎣⎡ −=

−

iii

drxl

rlbKd

drxl

rlbKd

baFd

drxl

rlbaFd 111

11,1,1

(2.46)

Figura 2.7 – gráfico da diferença entre as integrais elípticas do primeiro tipo incompleta e completa,

para ui = 0,6; ξ1 = 1,0; elemento a 45o em relação à horizontal (0 ≤ r ≤ ui – trecho cujos pontos se

enquadram no caso 2 – expressão (A.1)).


Em (2.46), refere-se à integração sobre o elemento Γi, de comprimento l, em

uma das extremidades do qual está localizado o ponto fonte. Adota-se uma

coordenada no espaço, r, variando de 0, na extremidade onde está o ponto fonte, até l,

na outra extremidade. A função de interpolação – (l – r)/l – é relativa ao nó onde se

localiza o ponto fonte. Para a função relativa ao outro nó – r/l –, pode-se integrar

numericamente, visto que, como essa função se anula em r = 0, o produto dela pela

integral elíptica possui limite finito neste ponto, conforme mostra a figura 2.8.

Figura 2.8 – produto da função de interpolação r/l pela integral elíptica incompleta de primeiro tipo,

com os mesmos dados numéricos da figura 2.7

2) A primeira parcela do segundo membro de (2.46) pode então ser calculada

numericamente.

O problema passa a ser o cálculo da segunda parcela, também singular em

r = 0. Pode-se mostrar que é possível se aproximar a integral elíptica completa de

primeiro tipo, para r próximo de zero, por um logaritmo, conforme equação (2.47),

apresentada em BREBBIA et al. [31] e figura 2.9. Portanto, somando-se e subtraindo-

se da segunda parcela do segundo membro de (2.46) esta aproximação, obtém-se a

expressão indicada em (2.48).

( )⎟⎟⎠

⎞⎜⎜⎝

⎛−−

−=⎟⎟⎠

⎞⎜⎜⎝

⎛ −22

2

22

22

16ln

2 rdr

rdd

drdK (2.47)


Figura 2.9 – gráfico da diferença entre os dois membros da equação (2.47), para os mesmos dados

numéricos do gráfico da figura 2.7.

( ) ( ) ( )

( )∫

∫∫

Γ

ΓΓ

−⎟⎟⎠

⎞⎜⎜⎝

⎛−−

−

−−

⎥⎦

⎤⎢⎣

⎡⎟⎟⎠

⎞⎜⎜⎝

⎛−−

+=−

i

ii

drxl

rlrd

r

rd

drxl

rlrd

rrd

dbKd

drxl

rlbKd

122

2

22

122

2

221

16ln

2

1

16ln

211

(2.48) A expressão (2.48) é introduzida em (2.46). A primeira parcela do segundo

membro de (2.48) pode ser integrada numericamente, sendo somada à primeira

parcela do segundo membro de (2.46). Utilizando-se a propriedade do logaritmo da

divisão, as parcelas devidas ao logaritmo do denominador do argumento do logaritmo

original são canceladas, já que não apresentam singularidade em r = 0. Até aqui,

então, tem-se a seguinte expressão:

( ) ( )

∫∫∫

Γ

ΓΓ

−

−−

−−

⎥⎦

⎤⎢⎣

⎡

−+=

−

i

ii

drxl

rlrd

r

drxl

rlrd

rbaFd

drxl

rlbaFd

122

1221

ln

ln,1,1

(2.49)

A segunda parcela do segundo membro de (2.49) pode ser integrada

analiticamente, pois possui uma singularidade integrável em r = 0.


3) Para facilitar esta integração, subtrai-se e soma-se do coeficiente do

logaritmo o seu valor em r = 0 (singularidade), obtendo-se então novamente uma

parcela a ser integrada numericamente e outra de fácil integração analítica, conforme

apresentado em (2.50).

( )( )

( )( )∫

∫∫

Γ

ΓΓ

==

+

+⎥⎦

⎤⎢⎣

⎡

==

−−

−=

−

−

i

ii

rdrrdrx

rdrrdrxx

lrl

rddrx

lrl

rdr

ln00

ln001ln

1

1122122

(2.50)

Aproveitando-se de que x1(r = 0) / d(r = 0) = 1/2, e substituindo a expressão

(2.50) em (2.49), escreve-se finalmente a expressão (2.51).

( ) ( ) ∫∫∫ΓΓΓ

−⎥⎦⎤

⎢⎣⎡ +

−=

−

iii

drrdrrxl

rlbaFd

drxl

rlbaFd 2

ln2

ln,1,111 (2.51)

Integra-se numericamente a primeira parcela do segundo membro em (2.51),

cujo gráfico é traçado na Figura 2.10, para os mesmos dados numéricos das figuras

2.7 a 2.9, e analiticamente a segunda parcela deste membro. Deste modo, tem-se o

resultado para a integral da integral elíptica incompleta de primeiro tipo ao longo do

elemento que contém r = 0. Para a integral completa o procedimento é o mesmo, a

partir do passo 2).

Figura 2.10 – integrando da primeira parcela do segundo membro de (2.51)


Para as integrais elípticas do primeiro tipo presentes no item A.2,

correspondentes à matriz H:

de (A.7) a (A.14) não é necessário este tratamento, visto que ao se integrar

no espaço são multiplicadas por C1, que é sempre nulo ao longo do elemento em que

está o ponto fonte, já que neste caso ∂r/∂n = 0;

para as demais expressões de A.2, o procedimento pode ser aquele

apresentado anteriormente, conforme mostra a expressão (2.52). Nesta expressão a

primeira parcela do segundo membro difere do primeiro membro de (2.46) apenas

pela constante que a multiplica; já a segunda parcela pode ser diretamente integrada

numericamente, visto que contém uma multiplicação por r, proporcionando um

integrando com variação semelhante à da função traçada no gráfico da Figura 2.8.

( )[ ]( ) ( ) ( )

( ) ( )∫

∫∫

Γ

ΓΓ

−−+

+−

=−

−

−+

i

ii

drl

rlbaFd

ruxn

drxl

rlbaFd

undrx

lrlbaF

rdduxnxn

f

ff

,sen4

,14

,

1

222

11

1122

22211

θξ

ξ

ξξ

(2.52) Na expressão (2.52), utilizaram-se as igualdades:

d 2 – r

2 = 4x1ξ1

x1 = r senθ

Outra singularidade que ocorre para r = 0 é encontrada nas expressões (A.7) a

(A.10), nos coeficientes que multiplicam as integrais elípticas de segundo tipo, em

que aparece um fator r no denominador. No entanto, essas expressões são sempre

multiplicadas por C1. Como aqui só se adotam elementos de geometria reta, sempre

que o ponto fonte faz parte do elemento, ∂r/∂n se anula e, portanto, C1 também se

anula. Deste modo, as expressões (A.7) e (A.10) não aparecem nas integrais ao longo

do elemento de contorno quando r = 0.

Deve ser feita uma ressalva que diz respeito às expressões (A.3), (A.9) e

(A.17). Nelas, aparece uma diferença entre integrais elípticas de primeiro tipo

incompletas. Esta diferença por si só já apresenta limite finito em r = 0, como se vê

na figura 2.11; portanto, para este caso, pode-se efetuar a integral numericamente de

modo direto.


Figura 2.11 – gráfico da diferença entre integrais elípticas existente em (A.3), (A.9) e (A.17), para os

mesmos dados numéricos das figuras 2.7 a 2.9, e uf = 0,2, traçado para 0 ≤ r ≤ uf e d ≥ ui – trecho cujos

pontos se enquadram no caso 4.

Por fim, trata-se das expressões referentes ao ponto fonte posicionado no eixo

de axissimetria. Para a expressão (2.43), a integral analítica em relação ao espaço é de

resolução imediata, conforme (2.53).

⎥⎦

⎤⎢⎣

⎡−=

−∫ 24coscos

41 2

0

xlxl

drl

rlrr

x

πθθ

π (2.53)

Em (2.53), o limite superior de integração x pode ser igual ao próprio

comprimento de elemento (l) ou a ui, dependendo do que for menor. A coordenada

dos pontos campo x1 é substituída por r cosθ, onde θ é a inclinação do elemento em

relação ao eixo horizontal (eixo 1). Esta expressão é relativa ao nó onde se localiza o

ponto fonte. Para o outro nó do elemento, basta trocar a função de interpolação (ao

invés de (l – r)/l, utiliza-se r / l).

Já para as expressões (2.44) e (2.45), não é necessário nenhum tratamento,

visto que há uma multiplicação por ∂r/∂n, que, como já foi citado, sempre se anula

para pontos pertencentes ao elemento que contém o ponto fonte, que é onde ocorre o

caso r = 0.


2.3.2. Singularidades nas frentes de onda da solução fundamental

Estas singularidades aparecem nas expressões do Anexo A quando r ou d

coincidem com alguma frente de onda (referente a ui ou uf). São sempre causadas

devido à presença no denominador da raiz quadrada de (d2 - uf/i2) ou (uf/i

2 - r2).

Podem ocorrer nas seguintes situações:

− Caso 2:

Se r = ui : expressões (A.7) e (A.11).

Se d = ui : expressões (A.7), (A.11), (A.15) e (A.19).

− Caso 3:

Nunca.

− Caso 4:

Se r = uf : expressões (A.9) e (A.13).

Se d = ui : expressões (A.9), (A.13), (A.17) e (A.21).

− Caso 5:

Se r = uf : expressões (A.10) e (A.14).

Se d = uf : expressões (A.10), (A.14), (A.18) e (A.22).

Como ilustração, é apresentado aqui o procedimento completo empregado em

duas situações, sendo as demais análogas a estas.

1a situação: expressão (A.7) com r = ui:

Nesta expressão, apenas a parcela indicada pela expressão (2.54) da integral

ao longo do contorno, proveniente da primeira parcela de (A.7), é singular em r = ui.

( )∫ −−

−−

Sl

j

ii

i dxCxxruurd

ud

0

1122222

22

)(η (2.54)

Na expressão (2.54), η

j(x) representa a função de interpolação no espaço

correspondente a qualquer dos nós extremos do elemento sobre o qual se está

integrando, sendo x a variável que percorre o elemento, no sentido de integração.

Deve ser observado que os limites da integral dizem respeito ao segmento no qual o

Implementação numérica do MEC 43 elemento de contorno foi dividido, sendo ls o comprimento deste segmento.

Novamente, aqui se utiliza o artifício da integral semi-analítica, através do qual se

chega à expressão (2.55).

( ) ∫∫∫ −+

−

−=

−−

−−

SSS l

i

i

l

i

i

l

j

ii

i dxru

uAdx

ru

uAxAdxxCxx

rurdu

ud

022

022

0

1122222

22 )()()()()(η

(2.55)

Em (2.55):

( )( )

)()()(4)()()()( '11

2/31

11222

22

iij

ii

ij

i

i uCuuxu

uAxCxxurd

udxA η

ξη

−

−=⇒−

−−=

(2.56)

e

1

1'1

)()(x

xCxC = (2.57)

A primeira parcela do segundo membro de (2.55) pode ser integrada

numericamente, já que o integrando possui limite finito em r = ui. Resta então a tarefa

de se integrar analiticamente a segunda parcela deste membro. Realiza-se para isso

uma mudança da variável de integração, em que x é substituído por r (distância entre

ponto fonte e ponto campo). No entanto, como se pode notar pela expressão (2.58), o

Jacobiano da transformação (derivada dx/dr) pode mudar de sinal no meio do

segmento (o ponto de mudança é aquele no qual o raio é mínimo). Deve-se identificar,

para cada segmento, se o raio só cresce (derivada positiva), só diminui (derivada

negativa), ou se a derivada muda de sinal. Observando-se o sinal do Jacobiano, o

sentido de integração é mantido, trocando apenas a variável de integração.

22rhr

rdxdr

−±= (2.58)

Em (2.58), hr é a menor distância entre o ponto fonte e a reta que contém o

segmento sobre o qual se está integrando. Procedendo-se a mudança de variável, a

integral analítica se torna elementar, como indica a expressão (2.59).


f

i

f

i

S r

rii

r

r ri

i

l

i

i

huhruA

hr

drr

ru

uAdx

ru

uA

⎥⎥⎦

⎤

⎢⎢⎣

⎡

−−

±=−

±

−=

−−∫∫ 22

221

2222

0

22sen)(

)()( (2.59)

Caso ocorra mudança de sinal no interior do segmento:

f

i

f

r

r

i

f

i

r

rri

ri

r

h ri

i

h

r ri

i

r

r ri

i

huhrsenuA

hr

rdrru

uA

hr

rdrru

uA

hr

rdrru

uA

⎥⎥⎦

⎤

⎢⎢⎣

⎡

−−

=

=−−

+−−

−=−

±

−

−

∫∫∫

22

221

222222222222

)(

)()()(

(2.60)

Expressões de A(x) e A(ui) ou A(uf) correspondentes às demais singularidades

na frente de onda em r são apresentadas no Anexo D (item D.1).

2a situação: expressão (A.7) com d = ui:

A singularidade está presente apenas na parcela (2.61), proveniente da

primeira parcela de (A.7).

( )∫ −−

−Sl

j

ii

i dxCxxudrdu

ru

0

1122222

22

)(η (2.61)

Utilizando-se o mesmo processo aplicado a (2.54), chega-se a (2.62).

( ) ∫∫∫ −+

−

−=

−−

− SSS l

i

i

l

i

i

l

j

ii

i dxud

uBdx

ud

uBxBdxCxx

udrdu

ru

022

022

0

1122222

22 )()()()(η

(2.62)

Em (2.62):

( )( )

)()()(4)()()()( '11

2/31

11222

22

iij

ii

ij

i

i uCuuxu

uBxCxxurd

ruxB η

ξη

−

=⇒−

−=

(2.63)

Implementação numérica do MEC 45 Para se efetuar a integral analítica da segunda parcela do segundo membro de

(2.62), troca-se a variável de integração de x para d. Como no procedimento anterior,

a derivada dd/dx (2.64) pode assumir valores positivos ou negativos, porém, neste

caso, não pode haver mudança de sinal no meio do segmento, visto que, em sua

extensão, d não pode ser igual a hd (menor distância entre o ponto mais afastado do

anel de fontes e a reta que contém o segmento), já que, para os casos 2 e 4, d deve ser

maior que ui e ui deve ser maior que hd, o mesmo acontecendo para o caso 5 em

relação a uf. A integral é calculada então conforme (2.65).

22dhd

ddxdd

−±= (2.64)

f

r

f

i

S r

ridi

r

r di

i

l

i

i udhduBddhd

dud

uBdx

ud

uB⎥⎦⎤

⎢⎣⎡ ⎟

⎠⎞⎜

⎝⎛ −+−±=

−

±

−=

− ∫∫ 222

2220

2ln)(

)()(

(2.65)

Expressões de B(x) e B(ui) ou B(uf) correspondentes às demais singularidades

na frente de onda em d são apresentadas no Anexo D (item D.2).

2.3.3. Singularidades para ponto campo no eixo de axissimetria

Em todas as expressões presentes no item A.2, há termos com (d2 - r2) ou

mesmo (d2 - r2)2 no denominador. Considerando que

11

22 4 ξxrd =− (2.66),

este termo se anula quando ξ1 ou x1 se anulam. No primeiro caso, o ponto fonte se

localiza no eixo de axissimetria, então se utilizam as expressões tridimensionais. Já o

segundo caso representaria um problema para a integração ao longo do contorno.

Porém, ao se realizar essas integrais, as expressões são sempre multiplicadas por x1.

Além disso, C1 contém um fator x1. Deste modo, o termo x1 no denominador,

proveniente de (d2 - r2) ou (d2 - r2)2, é sempre cancelado, visto que sempre que a

expressão contém (d2 - r2)2, ela é multiplicada também por C1.


2.3.4. Transformada polinomial de 3o grau para a integração

Os integrandos provenientes do processo semi-analítico, a serem integrados

numericamente, em muitos casos, apesar de não-singulares, apresentam um gradiente

elevado nas proximidades da singularidade do integrando original.

Para se obter uma maior precisão através da integração numérica, então, faz-se

uso de uma transformação de coordenadas proposta por TELLES [30], que visa à

concentração maior de pontos de integração na região de elevado gradiente.

As novas coordenadas naturais η dos pontos de integração são determinadas

através da expressão (2.67), em função das coordenadas naturais (γ ∈ [-1,1]) dos

pontos de integração do processo padrão de Gauss.

( ) dcba +++= γγγγη 23 (2.67)

O jacobiano da transformação é dado por (2.68).

( ) cbaG ++= γγγ 23 2 (2.68)

Os parâmetros a, b, c e d são determinados a partir de (2.69).

( )

( )

( )

2

2

31

3

13

1

γ

γ

γ

+=

−=

+=

−−=

−=

Qbd

Qrc

Qrb

Qra

(2.69)

Em (2.69), γ é o valor da coordenada para a qual ( ) ηγη = , sendo η a

coordenada natural do ponto do elemento mais próximo da singularidade.

( ) ( )r

pqqpqq21

3 323 32

+++−−+++−=

ηγ (2.70)


( ) ( )

( )( ) ( )[ ]2

2

3

1314213

1

211

21223

2121

η

ηηη

−+−+

=

⎥⎦

⎤⎢⎣

⎡−

+⎟⎟⎠

⎞⎜⎜⎝

⎛+

−−+

=

rrr

p

rrr

rq

(2.71)

Em (2.71) r é um parâmetro que varia de acordo com a distância mínima

entre a singularidade e o elemento sobre o qual se está integrando. No caso, esta

distância é sempre nula.

1)ln(0832,0893,0

)ln(24,085,00

=+=

+==

rDr

Drr

(2.72)

Em (2.72):

slRD min2= (2.73)

Rmin é a mínima distância entre a singularidade e o elemento e ls o

comprimento do segmento sobre o qual se está integrando.

Deste modo, ficam completamente definidas, pela equação (2.67), as

coordenadas naturais dos pontos que serão utilizados para a integração ao longo do

elemento de contorno. Os pesos são os mesmos dos respectivos pontos de Gauss.

se 0,0 ≤ D ≤ 0,05 se 0,05 ≤ D ≤ 1,3 se 1,3 ≤ D ≤ 3,618 se 3,618 ≤ D


2.4. Termos relativos a fontes pontuais

O vetor de termos conhecidos do sistema de equações representado pela

expressão (2.10) apresenta uma parcela (S(n+1)) que contém as contribuições, para

cada ponto fonte no contorno, das fontes pontuais presentes no modelo. Esta parcela,

correspondente a um nó i do contorno, devida a uma fonte localizada no ponto XS, é

calculada de acordo com a expressão (2.9).

Nesta expressão, aparece a solução fundamental do problema e a função que

representa o comportamento da excitação provocada pela fonte ao longo do tempo.

Considera-se aqui que esta função é definida através de valores discretos (fm), dados a

cada instante de tempo m de análise, interpolados linearmente para a obtenção de

valores em qualquer instante. Então, a expressão da função para um instante de tempo

τ qualquer é dada por (2.74), em que as funções de interpolação φ m são as mesmas

expressas em (2.3) e NT é o número total de instantes de tempo de análise.

∑=

=NT

m

mm ff1

)()( τφτ (2.74)

Do mesmo modo que ocorre com a função p(τ), o valor da função f(τ) em

qualquer instante é definido apenas pelos valores extremos deste intervalo, devido à

linearidade das funções de interpolação. Empregando-se a definição (2.74), a

expressão (2.9) pode ser escrita da maneira indicada em (2.75), considerando a

presença de apenas uma fonte pontual.

∑+

=

++ =1

1

)1()1(n

m

mmni

ni fwS (2.75)

Em (2.75):

ττφτξ dtXpwnt

minS

mni ∫

+

++ =

1

0

1*)1( )(),;,( (2.76)

Implementação numérica do MEC 49 Como em (2.76) as funções φ m(τ) assumem valores não-nulos apenas para τ ∈

[tm-1,tm+1], basta se efetuar a integração neste intervalo a fim de se obter mniw )1( + . Vale

aqui também a propriedade de translação indicada no item 2.1.3. Deste modo a cada

passo de tempo são calculados, do mesmo modo que para os elementos da matriz H,

dois termos, para o intervalo [t1,t2], sendo um deles (wI) somado a ni

ni ww =+ 1)1( e outro

(wF) somado a 12)1( −+ = ni

ni ww . Primeiramente, considerando a fonte localizada fora do

eixo de axissimetria, obtêm-se as expressões (2.77) e (2.78).

( )

( )( ) ( )( )τ

τ

ττφτξξ

dt

t

rτtcξxrτtcπ

c

dtXptXw

f

i

f

i

t

t

f

nn

t

t

iinnI

∫

∫

⎟⎟⎠

⎞⎜⎜⎝

⎛Δ

−

⎟⎠⎞

⎜⎝⎛ −−−−−

=

==

++

++

221

211

221

22

1*

1

414

)(),;,(,,

(2.77)

( )

( )( ) ( )( )τ

τ

ττφτξξ

dtt

rτtcξxrτtcπ

c

dtXptXw

f

i

f

i

t

t

i

nn

t

t

finnF

∫

∫⎟⎠

⎞⎜⎝

⎛Δ−

⎟⎠⎞

⎜⎝⎛ −−−−−

=

==

++

++

221

211

221

22

1*

1

414

)(),;,(,,

(2.78)

Procede-se então à mesma mudança de variáveis realizada para as expressões

de g, hi e hf, chegando-se a (2.79) e (2.80). Estas integrais podem ser desenvolvidas

analiticamente, para cada um dos mesmos seis casos indicados em 2.2.1, de acordo

com as distâncias entre o ponto fonte e a fonte pontual, obtendo-se as expressões

apresentadas no Anexo E.

( )( )( )

duudrutcπ

uutXw

i

f

u

u

fnI ∫ −−Δ

+−=

222222,,ξ (2.79)


( )( )( )

duudrutcπ

uutXw

i

f

u

u

inF ∫ −−Δ

−=

222222,,ξ (2.80)

Quando o ponto se localiza no eixo de axissimetria e a solução fundamental

passa então a se confundir com a solução fundamental tridimensional, têm-se as

seguintes expressões, válidas apenas para o caso 2b, de acordo com o exposto no item

2.2.2. Para os outros dois casos, têm-se valores nulos.

( )tcur

rtXw f

nI Δ

−=

πξ

41,, (2.81)

( )tcru

rtXw i

nF Δ−

=π

ξ4

1,, (2.82)


2.5. Esquema Teta

Como é mostrado no Capítulo 4, há problemas em que, utilizando-se o método

padrão de avanço no tempo descrito no item 2.1.2, pode ocorrer uma severa

instabilidade na resposta obtida, a partir de determinado tempo da análise. Esta

instabilidade, advinda do próprio método numérico empregado, além de limitar o

tempo total de análise para o qual se obtém resposta aceitável, pode estar

prejudicando até mesmo a resposta anterior a este instante, apesar de não se notar erro

relevante.

Esta dificuldade pode também se manifestar de acordo com o problema. Para

alguns casos, pode não aparecer instabilidade alguma, por um longo tempo, enquanto

para outros ela pode se pronunciar logo nos primeiros instantes da análise. Nos

problemas de domínio infinito, essencialmente, este problema raramente aparece, pela

dissipação da energia. Quando o meio é limitado, e ocorrem múltiplas reflexões em

seu interior, o erro acaba se acumulando, até eclodir na forma de instabilidade

pronunciada.

Para tentar solucionar este problema, é adotado aqui o chamado “Esquema

Teta”, apresentado por MANSUR et al. [32]. Neste esquema, a resposta, a cada passo

de tempo, inicialmente é calculada para um instante de tempo tn+θ = tn + θ.Δt, no

lugar de tn+1. Posteriormente, interpolam-se as respostas obtidas para p e q nesse

instante de tempo, de acordo com as funções de interpolações no tempo empregadas

(neste trabalho, linear para p e constante para q), chegando-se à solução para o

instante tn+1, de acordo com as expressões (2.83) e (2.84). Ressalta-se aqui que é

condição obrigatória para o esquema que θ seja sempre maior que a unidade, ou seja,

tn+θ > tn+1. O que se faz na realidade é uma suavização da resposta, levando-se em

conta para o cálculo no instante tn+1 a resposta no instante anterior e um instante

posterior, artificialmente incluído na análise. Esta suavização acaba evitando que a

instabilidade se pronuncie.

ni

ni

ni ppp

θθ

θθ 111 −

+= ++ (2.83)

θ++ = n

ini qq 1 (2.84)


As matrizes G e H passam então a ser calculadas para um tempo final tn+θ.

Portanto, nas expressões para o cálculo de seus termos, onde aparece tn+1, com a

adoção do Esquema Teta, aparecerá tn+θ, conforme expresso em (2.85).

)(

1 1

)(

1 1

)()()( θθθθθ +

= =

+

= =

+++ +−+= ∑∑∑∑ n

n

m

NN

j

mj

mnij

n

m

NN

j

mj

mnij

nn pHqG SyBxA (2.85)

Para as integrações realizadas no primeiro intervalo de tempo, a fim de se

montar as matrizes A(n+θ), B(n+θ) e a primeira parcela de G(n+θ)1 e H(n+θ)1, utiliza-se ti =

t1 e tf = t(1+θ). Para os demais instantes de tempo, utiliza-se, conforme procedimento

usual, ti = t1 e tf = t2. Em todos os casos, ao invés de tn+1, emprega-se tn+θ. O mesmo é

válido para os termos que formam o vetor S(n+θ).

Como é inferido através dos resultados mostrados no Capítulo 4, o Esquema

Teta mostrou-se eficaz para as aplicações com as quais este trabalho lidou. Porém,

pode apresentar um amortecimento numérico maior que o fornecido pelo avanço no

tempo padrão. Em alguns casos, no entanto, isto é necessário, sendo preferível um

amortecimento artificial um pouco maior à instabilidade da solução.

Acoplamento MEC-MEF 53

Capítulo 3

Acoplamento MEC-MEF

No presente capítulo são abordados alguns aspectos relativos ao acoplamento

de sistemas fisicamente distintos, ou fisicamente similares porém modelados por

métodos numéricos diferentes, sempre considerando problemas axissimétricos. O

acoplamento, nos casos aqui analisados, é feito através da interface entre os dois

sistemas.

Os métodos numéricos com os quais se trabalha aqui são o Método dos

Elementos de Contorno (MEC) e o Método dos Elementos Finitos (MEF). Os meios

físicos podem ser, por sua vez, acústicos ou elásticos. Um meio acústico, neste

trabalho, pode ser modelado tanto pelo MEF quanto pelo MEC; já um meio elástico é

modelado apenas pelo MEF.

Portanto, são possíveis neste caso seis tipos de acoplamento:

1. acústico-acústico MEC-MEF;

2. acústico-acústico MEC-MEC;

3. acústico-acústico MEF-MEF;

4. acústico-elástico MEC-MEF;

5. acústico-elástico MEF-MEF.

6. elástico-elástico MEF-MEF.

Serão enfocados aqui o primeiro e o quarto tipo, por tratarem de acoplamento

entre métodos numéricos distintos. Os demais apresentam procedimentos análogos

aos aqui detalhados.

Para implementação numérica do acoplamento é empregado sempre o

algoritmo iterativo, através do qual cada sistema é resolvido isoladamente, mantendo

suas características de resolução originais. A cada iteração, são verificadas as

variáveis que as condições de interface exijam que sejam iguais, até que se atinja a

convergência.

O acoplamento, além das vantagens citadas na introdução, permite a adoção

de condições inicias não-nulas em determinada região do domínio. Como foi

mencionado, o esquema desenvolvido para o MEC não comporta este tipo de análise.

Pode-se então modelar a região em questão pelo MEF, acoplando com a malha do

MEC.


3.1. Expressões obtidas para o MEF

Neste item, é brevemente apresentada a formulação do MEF, para problemas

acústicos e elásticos. Em ambos os casos, utiliza-se o método dos resíduos ponderados

para obtenção das equações finais. Ressalta-se que só são considerados meios com

comportamento linear.

3.1.1. Meio acústico

A equação básica da propagação de ondas em um meio acústico homogêneo é

aqui escrita na forma apresentada pela equação (3.1), mais adequada para a análise

por elementos finitos e também mais geral, já que a velocidade de propagação da

onda primária c, antes a única propriedade descritiva do meio, é desmembrada em

duas outras propriedades, a densidade ρ e o módulo de compressibilidade K.

02 =∇−+ pKpp &&& ςρ (3.1)

Pode-se notar que na equação (3.1) é considerado o amortecimento do meio,

através do coeficiente de amortecimento viscoso ζ. Esta equação pode ser escrita para

cada elemento, considerando-se o meio homogêneo em seu interior. Porém, ao se

escrever a equação (3.1), ao invés da equação (1.1), faz-se a correta consideração da

variação das propriedades físicas entre elementos distintos.

Deve-se ressaltar que é desprezada na equação (3.1) a presença de fontes

geradoras de pressão no interior do elemento.

Aplicando-se o método dos resíduos ponderados à equação (3.1), com

formulação fraca e funções de ponderação w para o domínio Ω, w para Γ1 (parcela do

contorno Γ em que é prescrita a pressão p ) e w para Γ2 (parcela em que é prescrito o

fluxo q , ou seja, a derivada da pressão na direção normal ao contorno), chega-se a

(3.2).

( ) ( ) ( ) 0ˆ

21

212 =Γ−+Γ−+Ω∇−+ ∫∫∫

ΓΓΩ

dwffdwppwdpKpp &&& ςρ (3.2)

Em (3.2), f = Kq. Na seqüência, a última parcela do termo entre parênteses no

primeiro membro de (3.2) é expandida e integrada por partes, conforme (3.3).

Acoplamento MEC-MEF 55 Assume-se que a função de ponderação é nula em Γ1, eliminando-se a integral nesta

parcela do contorno. Maiores detalhes acerca do procedimento empregado podem ser

encontrados em ZIENKIEWICZ e TAYLOR [6].

( ) ∫∫∫ΓΩΩ

Γ⎟⎠⎞

⎜⎝⎛

∂∂

+Ω∇⋅∇−=Ω∇ dnpKwdpwKdwpK .T2 (3.3)

Substituindo-se (3.3) em (3.2) e igualando w a w, pode-se escrever a

expressão (3.4).

( ) ( ) ( ) ( )∫∫∫∫

ΓΩΩΩ

Γ=Ω+Ω+Ω∇⋅∇2

2T dqwKdpwdpwdpwK &&& ρς (3.4)

Utilizando-se a formulação de Galerkin, as funções de ponderação w são

interpoladas pelas mesmas funções empregadas para se aproximar a variável básica p

no espaço, conforme descrito por (3.5) e (3.6).

∑=

=N

jjj tpXNtXp

1

)()(),( (3.5)

∑=

=N

jjj twXNtXw

1

)()(),( (3.6)

Em (3.5) e (3.6), N é o número de nós do elemento, X é um ponto genérico no

interior do elemento e t um instante de tempo qualquer da análise. Nj(X) representa a

função de interpolação relativa ao nó j, que deve apresentar valor unitário neste nó e

nulo nos demais. São empregadas sempre, neste trabalho, funções de interpolação

bilineares. pj(t) e wj(t) contêm o valor, respectivamente, da pressão e da função de

ponderação no nó j no instante de tempo t. Expandindo-se a expressão (3.4) e

empregando-se as aproximações contidas em (3.5) e (3.6), chega-se a (3.7).

∫ ∑∫ ∑ ∑

∫ ∑ ∑∫ ∑ ∑

Γ =Ω = =

Ω = =Ω = =

Γ=Ω∇∇+

+Ω+Ω

21

21 1

T

1 11 1

N

iii

N

i

N

jjjii

j

N

i

N

jjiij

N

i

N

jjii

dqwNKdNpNwK

dpNwNdpNwN &&& ςρ

(3.7)

Acoplamento MEC-MEF 56 Tomando-se a equação (3.7) N vezes (k = 1..N) e, para cada k, adotando-se

arbitrariamente wi = 1 para i = k, e wi = 0 para i ≠ k, chega-se à expressão (3.8).

∫∑∫

∑∫∑∫

Γ= Ω

= Ω= Ω

Γ=Ω∇∇+

+Ω+Ω

21

T

11

dfNpdNKN

pdNNpdNN

i

N

j

jji

N

j

jji

N

j

jji &&& ςρ

i = 1..N (3.8)

Em (3.8), qKf = . Tem-se, a partir de (3.8), um sistema de equações, para

cada elemento, que pode ser escrito através da expressão (3.9).

nnnn FKPPCPM =++ &&& (3.9)

Em (3.9) Pn é o vetor de pressões nodais no instante de tempo nt ; Fn é o vetor

de fluxos equivalentes nodais atuantes neste instante; e M, C e K são as matrizes de

massa, amortecimento e rigidez do elemento, respectivamente. Obtendo-se as

matrizes para cada elemento, pode-se chegar à matriz global do sistema por simples

agrupamento dessas matrizes.

Cada elemento das matrizes de massa M, C e K é fornecido pelas expressões

(3.10) a (3.12).

∫Ω

Ω= dNNm jiij ρ (3.10)

∫Ω

Ω= dNNc jiij ς (3.11)

∫Ω

Ω= dKk jTiij BB (3.12)

Nas expressões (3.10) a (3.12), Ni representa a própria função de interpolação

relativa ao nó i, enquanto Bi é um vetor dado por (3.13).

i

i

i

i

i NzNyNxN

∇≡⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

∂∂∂∂∂∂

=///

B (3.13)


Cada elemento do vetor de fluxos equivalentes nodais atuantes é fornecido por

(3.14).

∫Γ

Γ=2

2dfNf ii (3.14)

Até aqui, todo o desenvolvimento foi efetuado para o caso genérico

tridimensional. Considerando-se agora a axissimetria do problema, as expressões

(3.10) a (3.12), além de (3.14), são reescritas nas expressões (3.17) a (3.20), onde se

adotam as coordenadas r (raio em relação ao eixo de axissimetria) e z (altura). Esta

passagem é realizada através da substituição do domínio e contorno infinitesimais

originais dΩ e dΓ pelas expressões (3.15) e (3.16), seguida de integração ao longo de

θ.

axiddrdzdrdrd Ω==Ω θθ (3.15)

axiddrd Γ=Γ θ (3.16)

Nessas expressões, o domínio infinitesimal dΩaxi representa uma área (dr.dz),

que, ao ser girada em torno do eixo de axissimetria, gera um toro infinitesimal, que

pode representar o domínio infinitesimal tridimensional do qual se estava tratando até

aqui.

∫Ω

Ω=axi

axijiij drNNm ρπ2 (3.17)

∫Ω

Ω=axi

axijiij drNNc ςπ2 (3.18)

∫Ω

Ω=axi

axijT

iij drKk BBπ2 (3.19)

∫Γ

Γ=axi

drfNf ii 22π (3.20)

Em (3.19):

ii

ii zN

rNNB ∇≡⎥

⎦

⎤⎢⎣

⎡∂∂∂∂

=//

(3.21)


Montadas as matrizes de massa, amortecimento e rigidez, trata-se do avanço

no tempo. É escolhido para efetuar este avanço o Método de Newmark, através do

qual os vetores de aceleração e velocidade são expressos, em função do vetor de

deslocamentos, pelas expressões (3.22) e (3.23), respectivamente. Para os exemplos

estudados neste trabalho, os parâmetros ajustáveis γ e β foram tomados como 1/2 e

1/4, respectivamente, o que equivale à regra trapezoidal.

( ) nnn1n1n ppppp &&&&& ⎟⎟⎠

⎞⎜⎜⎝

⎛−−−−= ++ 1

β21

βΔt1

βΔt1

2 (3.22)

( ) nnn1n1n β2γ1

βγ1

βΔtγ ppppp &&&& tΔ⎟⎟

⎠

⎞⎜⎜⎝

⎛−+⎟⎟

⎠

⎞⎜⎜⎝

⎛−−−= ++ (3.23)

Em (3.22) e (3.23), p representa o vetor de pressões em um determinado

instante de tempo caracterizado pelo índice subscrito; Δt é o intervalo de tempo da

análise por elementos finitos.

Introduzindo-se estas expressões em (3.9), escrita para o tempo tn+1, chega-se

ao sistema representado pela expressão (3.24).

BpA 1nˆˆ =+ (3.24)

Em (3.24) A e B representam, respectivamente, a matriz e o vetor efetivos

do sistema, dados por (3.25) e (3.26).

KCMA +⎟⎟⎠

⎞⎜⎜⎝

⎛+⎟⎟

⎠

⎞⎜⎜⎝

⎛=

βΔtγ

βΔt1ˆ

2 (3.25)

⎥⎦

⎤⎢⎣

⎡⎟⎟⎠

⎞⎜⎜⎝

⎛−−⎟⎟

⎠

⎞⎜⎜⎝

⎛−−+

+⎥⎦

⎤⎢⎣

⎡⎟⎟⎠

⎞⎜⎜⎝

⎛−+++= +

nnn

nnn21

β2γ1

βγ1

βΔtγC

1β2

1βΔt1

βΔt1ˆ

ppp

pppMFB

&&&

&&&n

(3.26)

Através deste esquema, o problema é solucionado a cada passo de tempo,

tendo já como conhecidos os valores da pressão e suas derivadas primeira e segunda

no tempo em toda a malha para o instante de tempo imediatamente anterior.


3.1.2. Meio elástico

Aqui são deduzidas as expressões utilizadas no MEF quando da análise de um

meio com propriedades elásticas. Enquanto que em um meio acústico a variável

básica (pressão) é um escalar, para um meio elástico ela se constitui de um vetor,

composto dos deslocamentos em todas as direções, de acordo com a dimensão do

problema, tornando a matemática do problema um pouco mais complexa.

No desenvolvimento aqui exposto, tomam-se como ponto de partida as

equações de equilíbrio da elasticidade expressas por (3.27), às quais é aplicado o

método dos resíduos ponderados. Estas equações apresentam como variáveis também

as tensões. Através de algumas manipulações apresentadas neste item, no entanto,

chega-se às expressões finais tendo como incógnitas apenas os deslocamentos.

Ressalta-se que, por razões de praticidade, o desenvolvimento é feito para o

caso bidimensional (elasticidade plana) e que, por procedimentos análogos, obtêm-se

as expressões para o caso axissimétrico.

0

0

2

2

2

2

=∂

∂−

∂

∂−+

∂

∂+

∂

∂

=∂

∂−

∂∂

−+∂

∂+

∂∂

tu

tu

byx

tu

tu

byx

yyy

yxy

xxx

xyx

ςρστ

ςρτσ

(3.27)

Trata-se primeiro das três primeiras parcelas das expressões (3.27). Elas são

multiplicadas pelas funções de ponderação e integradas no domínio, conforme

exposto em (3.28), na qual já está feita a integração por partes.

( ) ( )[ ]

( )∫ ∫

∫

∫

Ω Ω

Γ

Ω

Ω++Ω⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎟⎠

⎞⎜⎜⎝

⎛∂

∂+

∂

∂+⎟⎟

⎠

⎞⎜⎜⎝

⎛∂

∂+

∂∂

−

−Γ+++=

=Ω⎥⎥⎦

⎤

⎢⎢⎣

⎡⎟⎟⎠

⎞⎜⎜⎝

⎛+

∂

∂+

∂

∂+⎟⎟

⎠

⎞⎜⎜⎝

⎛+

∂

∂+

∂∂

dwbwbdy

wx

wy

wx

w

dwnnwnn

dwbyx

wbyx

yyxxy

yy

xyx

xyx

x

yyyxxyxyxyxx

yyyxy

xxxyx

σττσ

σττσ

σττσ

(3.28) Nota-se que a função de ponderação aqui se constitui de um vetor w de

componentes wx e wy. bx e by são as componentes do vetor b representativo das forças

Acoplamento MEC-MEF 60 de volume presentes no meio, enquanto o vetor n, de componentes nx e ny, representa

o vetor normal ao contorno.

A expressão (3.28) pode ser escrita na forma vetorial, conforme (3.29).

( ) ∫∫∫

ΓΩΩ

Ω+Ω+Ω− ddd TTT wtwbσwL (3.29)

Em (3.29), σ é o vetor que contém as tensões normais σx e σy e a tensão

cisalhante τxy, enquanto o vetor t contém as componentes das forças de superfície. A

matriz de transformação L é dada por (3.30).

⎥⎥⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢⎢⎢⎢

⎣

⎡

∂∂

∂∂

∂∂

∂∂

=

xy

y

x0

0

L (3.30)

São utilizadas então as relações constitutivas σ = D.ε para se chegar à

expressão (3.31). D é a matriz constitutiva expressa por (3.32), para estado plano de

tensão, e (3.33), para estado plano de deformação, enquanto ε é o vetor que contém as

deformações εx, εy e εxy, dado pela expressão ε = Lu, sendo u o vetor que contém os

deslocamentos ux e uy.

( ) ( ) ∫∫∫

ΓΩΩ

Ω+Ω+Ω− ddd TTT wtwbuDw LL (3.31)

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−−

=2/)1(00

0101

1 2

νν

ν

νED (3.32)

⎥⎥⎥

⎦

⎤

⎢⎢⎢

⎣

⎡

−−−

−

−−−

=)1(2/)21(00

01)1/(0)1/(1

)21)(1()1(

νννν

νν

νννED (3.33)

Obtida então a expressão (3.31), o passo seguinte é o emprego da aproximação

numérica. Assim como para o meio acústico, utiliza-se aqui também a formulação de

Acoplamento MEC-MEF 61 Galerkin. Os deslocamentos, bem como as funções de ponderação, são expressos,em

qualquer ponto no interior do elemento, em função dos valores nodais, segundo as

expressões (3.34) e (3.35), nas quais N é o número de nós do elemento em questão.

( ) ∑=

=N

jtXtX

1)()(, jj uNu (3.34)

( ) ∑=

=N

jtXtX

1)()(, jj wNw (3.35)

A matriz de funções de interpolação Ni é dada por (3.36), onde Ni é a função

de interpolação relativa ao nó i.

⎥⎦

⎤⎢⎣

⎡=

i

ii N

N0

0N (3.36)

As mesmas aproximações são utilizadas para a velocidade e a aceleração

presentes nas expressões (3.27). Introduzindo-se as expressões (3.34) e (3.35) em

(3.31), e esta em (3.27), além de na última, serem aplicadas as aproximações para

velocidade e aceleração, chega-se à expressão (3.37). Aqui já é considerada a

ponderação no contorno, como no caso acústico, admitindo que as funções de

ponderação anulam-se na parcela do contorno em que, no caso, o deslocamento é

prescrito (Γ1).

∫ ∑∫ ∑∫ ∑∑

∫ ∑ ∑∫ ∑ ∑

Γ =Ω =Ω ==

Ω = =Ω = =

Γ+Ω=Ω⎟⎟

⎠

⎞

⎜⎜

⎝

⎛⎟⎟⎠

⎞⎜⎜⎝

⎛+

+Ω+Ω

21

2111

1 11 1

N

iii

TN

iii

TN

jjj

TN

iii

j

N

i

N

jiijj

N

i

N

jjii

ddd

dd

wNtwNbuNDwN

uNwNuNwN

LL

&&& ςρ

(3.37)

Toma-se então a expressão (3.37) N vezes (k = 1 .. N), de modo que para cada

valor de k toma-se primeiramente a equação (3.37) com ⎥⎦

⎤⎢⎣

⎡=

01

iw e depois ⎥⎦

⎤⎢⎣

⎡=

10

iw

para i = k, mantendo-se sempre ⎥⎦

⎤⎢⎣

⎡=

00

iw para i ≠ k. Deste modo, é obtido um sistema

de 2N equações, descrito pela expressão (3.38) e, de forma mais resumida, por (3.39).


( ) ( ) ∫∫∑∫

∑∫∑∫

ΓΩ= Ω

= Ω= Ω

Γ+Ω=Ω+

+Ω+Ω

2

21

11

ddd

dd

Ti

Ti

N

jjj

Ti

N

jjji

N

jjji

tNbNuNDN

uNNuNN

LL

&&& ςρ

(i = 1..N) (3.38)

nnnn RKUUCUM =++ &&& (3.39)

Em (3.39), o vetor Un contém os deslocamentos nodais do elemento no

instante de tempo n, enquanto Rn contém as cargas nodais equivalentes relativas a

este instante.

As sub-matrizes de massa, amortecimento e rigidez são calculadas pelas

expressões (3.40) a (3.42). Percorrendo-se todos os nós, são montadas as matrizes

globais. O sub-vetor de cargas nodais é dado por (3.43).

∫Ω

= dΩjT

iij NNM ρ (3.40)

∫Ω

= dΩjT

iij NNC ς (3.41)

∫Ω

Ω= djT

iij BDBK (3.42)

∫∫ΓΩ

Γ+Ω=2

2dd Ti

Tii tNbNR (3.43)

A matriz Bi é dada por (3.44).

i

ii

i

i

i

xN

yN

yN

xN

NB L≡

⎥⎥⎥⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢⎢⎢⎢

⎣

⎡

∂∂

∂∂

∂∂

∂∂

= 0

0

(3.44)

Para o caso axissimétrico, o procedimento é semelhante, porém parte-se de

equações distintas, já deduzidas para esta situação, conforme demonstrado em

ZIENKIEWICZ e TAYLOR [6]. O comportamento do sistema neste caso pode ser

descrito também pelo deslocamento em apenas duas direções (ur e uz, utilizando o

mesmo sistema de coordenadas do caso acústico); portanto, o problema continua

Acoplamento MEC-MEF 63 sendo modelado apenas em duas dimensões. Porém, devem-se levar em conta quatro

componentes de tensão (σT = [σr σz σθ τrz]), ou deformação (εT = [εr εz εθ γrz]), ao

contrário das três consideradas para o caso bidimensional (neste caso, ou a tensão, no

estado plano de tensão, ou a deformação, no estado plano de deformação, na direção z

sempre se anula e a equação correspondente a esta direção fica desacoplada do

sistema). Com isso, a matriz constitutiva D e a matriz B são modificadas, de acordo

com (3.47) e (3.48).

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

−−

−−

−+=

νννν

νννννν

νν2/1000010101

)21)(1(ED (3.45)

⎥⎥⎥⎥

⎦

⎤

⎢⎢⎢⎢

⎣

⎡

∂∂∂∂

∂∂∂∂

=

rNzNrN

zNrN

ii

i

i

i

i

//0//0

0/

B (3.46)

As matrizes M, C e K são calculadas para o caso axissimétrico a partir das

equações (3.47) a (3.49). Nestas expressões está incluída a integração ao longo de um

anel em volta do eixo de axissimetria. Do mesmo modo, o vetor de cargas nodais

passa a ser calculado a partir de (3.50).

∫Ω

Ω=axi

axijT

iij drNNM ρπ2 (3.47)

∫Ω

Ω=axi

axijT

iij drNNC ςπ2 (3.48)

∫Ω

Ω=axi

axijT

iij drBDBK π2 (3.49)

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡Γ+Ω= ∫∫

ΓΩ axiaxi

axiTiaxi

Tii drdr tNbNR π2 (3.50)

O avanço no tempo para o caso elástico é realizado da mesma forma que para

o caso acústico (método de Newmark com regra trapezoidal), conforme expressões

(3.22) a (3.26).


3.2. Equações governantes do acoplamento

Para que seja fisicamente válido o acoplamento entre sistemas resolvidos

independentemente, devem, na interface, ser respeitadas determinadas condições, que

relacionam variáveis destes sistemas.

3.2.1. Acoplamento acústico-acústico MEC-MEF

Neste item, é tratado o acoplamento entre sistemas fisicamente similares, do

tipo acústico-acústico, ou fluido-fluido, modelados por métodos numéricos distintos:

um deles é modelado através do MEC e o outro através do MEF.

Neste caso, devem ser obedecidas as seguintes condições ao longo da interface

de acoplamento ΓI:

( ) ( )iI

FiI

C PP = (3.51)

( ) ( )iI

FiI

C QQ ~~−= (3.52)

Nas equações (3.51) e (3.52) P e Q representam vetores contendo,

respectivamente, os valores de pressão e fluxo. O índice superior esquerdo I indica

que neste vetor estão contidos apenas os valores na interface ΓI; o índice inferior

esquerdo indica o método empregado na obtenção deste vetor (C para o MEC e F para

o MEF); por fim, o índice fora dos parênteses informa o nó i da interface, o que

significa que estas condições devem ser respeitadas em cada nó.

O fluxo ( )iIC Q~ em (3.52), que representa o fluxo nodal equivalente no nó i

obtido pelo MEC, é calculado através dos valores ( )iI

C Q originalmente fornecidos

pelo método. Estes representam o valor do fluxo distribuído ao longo do contorno em

determinado nó. A transformação é obtida pela expressão (3.53).

( ) ( ) ( ) ( ) ( )i

jI

CjCTF

i

jI

CjCTFi

IC

j

j

j

jdd

⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛Γ+

⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛Γ= ∫∫

+

+

Γ

+Γ+

Γ

Γ

1

1 11 rr~ QNNQNNQ (3.53)


Em (3.53) os índices Γj e Γj+1 relacionam-se aos elementos de contorno

adjacentes ao nó i. ( )j

IC ΓQ e ( )

1+Γ j

IC Q são agora vetores que contêm, cada um, os

valores obtidos diretamente pelo MEC para o fluxo distribuído nos nós extremos de

seu respectivo elemento. Por fim, as variáveis N são vetores contendo os valores das

funções de interpolação relativas aos nós extremos do elemento ao qual eles se

referem. O sub-índice destes vetores indica o método do qual se está tratando. A

expressão (3.54) exemplifica um vetor N relativo a um elemento j, de nós extremos i e

i+1, com as funções de interpolação utilizadas no MEC, sendo X um ponto qualquer

do contorno que está sendo integrado. Cada integral presente em (3.53) fornece um

vetor com os fluxos nodais equivalentes nos dois nós extremos do elemento. Para se

obter ( )iIC Q~ toma-se apenas o valor relativo ao nó extremo i.

( ) [ ])()( 1 XNXN iC

iCjC

+=N (3.54)

3.2.2. Acoplamento acústico-elástico MEC-MEF

Neste item, é tratado o acoplamento acústico-elástico entre sistemas

discretizados pelo MEF (elástico) e pelo MEC (acústico). O exemplo clássico deste

tipo de acoplamento é o acoplamento fluido-estrutura, já que o fluido pode ser

considerado um meio acústico (propagam-se apenas ondas de compressão, já que o

material do meio não resiste ao cisalhamento), enquanto a estrutura é composta de

material elástico (metal, concreto, polímeros, entre outros diversos).

Ao longo da interface ΓI devem ser obedecidas as condições (3.55) a (3.57):

( ) ( )iI

CiI

FN QUρ1

=&& (3.55)

( ) 0=iI

FT F (3.56)

( ) ( )iI

CiI

F FF −= (3.57)

As notações são equivalentes àquelas empregadas no item 3.2.1. O operador

N, aplicado a um vetor, extrai deste a sua componente na direção normal ao contorno,

Acoplamento MEC-MEF 66 por simples decomposição vetorial, enquanto através do operador T, do mesmo modo,

é obtida a componente na direção tangencial.

O vetor U&& contém as acelerações obtidas pelo MEF, enquanto F apresenta as

cargas nodais equivalentes fornecidas por este método ou, indiretamente, através do

MEC. O vetor FIC , que contém valores de força nodal equivalente, é calculado a

partir dos valores de pressão nodal obtidos na interface diretamente através do MEC.

A transformação é realizada através da expressão (3.58).

( ) ( ) ( ) ( ) ( )i

jI

CjCTF

i

jI

CjCTFi

IC

j

j

j

jdd

⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛Γ+

⎟⎟⎟

⎠

⎞

⎜⎜⎜

⎝

⎛Γ= ∫∫

+

+

Γ

+Γ+

Γ

Γ

1

1 11 rr PnNNPnNNF (3.58)

Pode-se observar que, em relação à expressão (3.53), aparece o vetor n, que

exprime a direção da normal ao elemento de contorno. Como são utilizados sempre

elementos lineares, este vetor é constante ao longo deles, sendo dado pela expressão

(3.59). Nesta expressão, θ é o ângulo que o elemento faz com a horizontal,

respeitando o sentido de integração.

⎥⎦

⎤⎢⎣

⎡−

=θ

θcos

senn (3.59)

Como a pressão é um escalar, ao ser integrada ao longo do elemento de

contorno, fornece uma força normal ao elemento, de acordo com a própria definição

de pressão. O vetor n tem como objetivo decompor esta força nas direções das

coordenadas axissimétricas r e z, gerando como resultado uma grandeza vetorial.

Do mesmo modo, como as forças nodais são grandezas vetoriais, NF não mais

se apresenta como um vetor, mas sim como uma matriz, dada por (3.60). Nesta

expressão, exemplifica-se uma matriz relativa a um elemento j, de nós extremos i e

i+1. As integrais em (3.58), de modo equivalente às de (3.53), fornecem um vetor

com as componentes nas duas direções da força equivalente nodal para os dois nós

extremos do elemento. Para se chegar a ( )iI

C F tomam-se apenas os valores relativos

ao nó extremo i.

( ) ⎥⎦

⎤⎢⎣

⎡=

+

+

)(0)(00)(0)(1

1

XNXNXNXN

iF

iF

iF

iF

jFN (3.60)


3.3. Procedimentos numéricos para o acoplamento iterativo

Neste item é descrito passo a passo o processo de resolução de um problema

de acoplamento através do método iterativo, introduzido por SOARES JR. e VON

ESTORFF [16]. São apresentados separadamente os procedimentos para o

acoplamento acústico-acústico e acústico-elástico entre sistemas resolvidos pelo MEF

e pelo MEC.

O princípio do método estabelece a resolução independente dos meios

acoplados, verificando-se a cada iteração se a convergência é atingida. Podem ser

adotadas para cada meio diferentes discretizações temporais e espaciais,

flexibilizando o modelo utilizado, de modo a se tomar partido das vantagens

oferecidas por cada método numérico empregado.

3.3.1. Acoplamento acústico-acústico MEC-MEF

No primeiro instante de tempo, calculam-se as atrizes M, C e K para o sistema

resolvido pelo MEF, bem como as matrizes A e B relativas ao sistema modelado

através do MEC e as primeiras parcelas das primeiras matrizes de convolução G e H e

do vetor S (equação (2.10)).

É resolvido então, para o primeiro passo de tempo, em primeiro lugar, o meio

modelado pelo MEF, arbitrando-se na interface de acoplamento fluxo nodal

equivalente nulo. Através da aplicação de (3.24) o sistema é resolvido, considerando-

se as condições iniciais. Com isso, são obtidas as pressões 1+kt

IF F

P na interface, onde tF

representa o instante de tempo da discretização de elementos finitos para o qual estão

sendo calculadas (no caso, t é igual ao segundo instante de tempo t2), e k o número da

iteração (no caso, k é igual a 1). Para acelerar a convergência, adota-se um parâmetro

de relaxação, produzindo-se a expressão (3.61), segundo a qual o valor no instante de

interesse é ponderado pelo valor obtido pela iteração anterior (no caso de ser a

primeira iteração, pelo valor final no instante de tempo anterior). Segundo os testes

realizados, o valor de α = 0,5 foi considerado adequado. Mais detalhes sobre este

parâmetro podem ser obtidos em ELLEITHY et al. [33], para o caso elastoestático.

( ) k

tI

Fkt

IF

kt

IF FFF

PPP αα −+= ++ 111 (3.61)


Aqui, entra um importante aspecto relativo às discretizações espacial e

temporal adotadas. Caso se utilizem discretizações espaciais distintas, os valores

obtidos pelo MEF devem ser interpolados para se obter os valores nos nós da interface

pertencentes ao modelo do MEC ( 1+kt

IC F

P ). Caso as discretizações temporais também

sejam diferentes, é necessária uma extrapolação no tempo para se obter os valores de 1+k

tI

C CP , onde tC corresponde ao instante de tempo da discretização de elementos de

contorno para o qual se está resolvendo o problema. No caso, considera-se o intervalo

de tempo do MEF menor que o do MEC, o que é em geral válido para problemas com

discretização espacial semelhante. Caso o passo de tempo do MEF seja maior que o

do MEC, o procedimento é análogo.

Esta extrapolação é representada pela equação (3.62), para a qual se considera

variação linear da pressão entre dois instantes de tempo consecutivos (figura 3.1).

( )

ηη−

−= Δ−

++

1

11 CCF

C

ttI

Ckt

ICk

tI

C

PPP (3.62)

Figura 3.1 – extrapolação da pressão obtida pelo MEF

Em (3.62), η = (tC – tF)/ΔtC, onde ΔtC é o intervalo de tempo da análise de

elementos de contorno.

As pressões obtidas na interface são impostas como condições de contorno

para o meio modelado pelo MEC nos nós correspondentes. Então este sistema é

resolvido, através de (2.10). Os valores calculados para o fluxo na interface de

acoplamento são em seguida trazidos de volta para o instante de tempo “atual” de

elementos finitos (tF). Como o fluxo é considerado constante no tempo para cada

1+kt

IC C

P

1+kt

IC F

P

CC ttI

C Δ−P

Acoplamento MEC-MEF 69 intervalo, têm-se imediatamente os valores para o fluxo distribuído, em cada nó, de

acordo com (3.63) e ilustrado na figura 3.2.

11 ++ = k

tI

Ckt

IC CF

QQ (3.63)

Figura 3.2 – interpolação do fluxo obtido pelo MEC

Estes valores são, a seguir, transformados de acordo com a expressão (3.53)

para valores de fluxo nodal equivalente, sendo estes então comparados com aqueles

arbitrados como condição de contorno no MEF, após a interpolação espacial, se

necessária. Aqui se faz necessário o teste de convergência.

O teste de convergência neste trabalho é realizado através do erro relativo

entre as normas dos vetores que contêm os valores da interface a serem testados. A

condição para que haja convergência é expressa por (3.64), sendo ε um parâmetro de

convergência escolhido.

εε ≤−

≤−

+

++

+

++

1

11

1

11

~

~~

~

~~

kt

IC

kt

IF

kt

IC

kt

IF

kt

IF

kt

IC

F

FF

F

FF eQ

QQ

Q

QQ (3.64)

Caso a convergência não seja atingida, deve ser efetuada uma nova iteração,

tendo como ponto de partida a adoção de condições de contorno na interface para o

sistema modelado pelo MEF. Aqui, segundo os testes realizados, opta-se novamente

pela utilização de um parâmetro de relaxação α , fazendo com que seja tomada como

condição de contorno na interface uma ponderação entre os fluxos nodais

equivalentes obtidos na iteração anterior por elementos de contorno e os fluxos

11 ++ = kt

IC

kt

IC CF

QQ

Acoplamento MEC-MEF 70 tomados naquela iteração como condição de contorno para elementos finitos, segundo

a equação (3.65).

( ) 112 ~1~~ +++ −+= kt

IF

kt

IC

kt

IF FFF

QQQ αα (3.65)

Segue-se então com a resolução do sistema pelo MEF, sendo o processo

repetido até que se atinja a convergência. Então, passa-se para o próximo instante de

tempo, sendo adotado inicialmente como condição de contorno na interface o fluxo

nodal obtido pela expressão (3.65) após a última iteração.

3.3.2. Acoplamento acústico-elástico MEC-MEF

Os procedimentos no primeiro passo de tempo são os mesmos executados para

o acoplamento acústico-acústico.

Também é inicialmente resolvido o sistema modelado pelo MEF. Neste caso,

são arbitradas de início condições de contorno de forças nodais nulas na interface de

acoplamento. Após a solução, são obtidos os valores de deslocamento, velocidade e

acelerações na interface. Novamente, é empregado um parâmetro de relaxação α, de

acordo com a expressão (3.66), para a aceleração, que é a variável de interesse.

( ) kt

IF

kt

IF

kt

IF FFF

UUU &&&&&& αα −+= ++ 111 (3.66)

No caso da primeira iteração em determinado instante de tempo, toma-se

como valor da aceleração na iteração anterior o valor final obtido para o instante de

tempo anterior. No primeiro passo de tempo, toma-se a aceleração inicial, calculada a

partir da velocidade e do deslocamento iniciais.

O vetor obtido em (3.66) é então aplicado à equação (3.55), através da qual se

chega a valores de fluxo através da interface. Novamente, aqui é necessário

considerar-se a diferença nas discretizações espacial e temporal. Caso a primeira seja

distinta, ou seja, os nós dos dois modelos acoplados não coincidam, deve-se realizar

uma interpolação para que sejam calculados os valores nos nós de interface do

sistema modelado pelo MEC. Além disso, no caso de a discretização temporal ser

distinta, os valores devem ser extrapolados (considerando intervalo de tempo do MEF

Acoplamento MEC-MEF 71 menor que o do MEC). Como o fluxo é considerado constante entre dois instantes de

tempo consecutivos, é utilizada a expressão (3.67).

11 ++ = k

tI

Ckt

IC FC

QQ (3.67) Estes valores são em seguida adotados como condições de contorno na

interface de acoplamento para a resolução do sistema modelado pelo MEC. Deste

modo, são obtidos valores de pressão na interface. Empregando-se a equação (3.58),

esses valores são transformados em forças nodais equivalentes, obtidas pelo MEC.

Novamente, então, esses valores devem ser trazidos para o instante de tempo da

análise por elementos finitos, de acordo com a expressão (3.68), que considera

variação linear de pressão e, em conseqüência, de forças nodais equivalentes, entre

instantes consecutivos de tempo da análise.

( ) 11 1 +

Δ−+ −+= k

tI

CttI

Ckt

IC CCCF

FFF ηη (3.68) Em (3.68), η tem o mesmo significado daquele presente em (3.62).

Finalmente, os valores resultantes de (3.68) podem ser comparados com

aqueles adotados ao início da iteração como condições de contorno para o sistema

resolvido pelo MEF, através de teste de convergência semelhante àquele indicado por

(3.64), trocando-se apenas a variável a ser testada (F ao invés de Q~ ).

Caso a convergência não seja atingida, uma nova iteração é realizada.

Novamente, de acordo com os testes efetuados, é utilizado um parâmetro de relaxação

α para a adoção das condições de contorno da interface, para o meio modelado pelo

MEF, na próxima iteração, de acordo com a expressão (3.69).

( ) 112 1 +++ −+= k

tI

Fkt

IC

kt

IF FFF

FFF αα (3.69) Segue-se então com a resolução do meio modelado pelo MEC, sendo o

processo repetido até que se atinja a convergência. Passa-se em seguida para o

instante de tempo subseqüente, adotando-se inicialmente como condição de contorno

na interface as forças nodais obtidas pela expressão (3.69), após a última iteração do

instante anterior. Deve-se verificar a cada mudança de instante de tempo de resolução

se este instante coincide com um instante de resolução de ambos os meios ou de

apenas um deles.

Aplicações 72

Capítulo 4

Aplicações

Neste capítulo são apresentadas algumas aplicações nas quais os métodos aqui

desenvolvidos podem ser empregados. Este emprego se dá através de rotinas em

FORTRAN implementadas ao longo deste trabalho.

Primeiramente, trata-se de aplicações em meios exclusivamente acústicos, de

características axissimétricas, discretizados somente por elementos de contorno, nos

quais se aplicam os procedimentos apresentados no Capítulo 2.

Em seguida, faz-se uso do acoplamento acústico-acústico discutido no

Capítulo 3, para se resolver problemas que dizem respeito a meios exclusivamente

acústicos, nos quais se deseja discretizar uma região por elementos finitos e outra por

elementos de contorno. Nestes casos, incluem-se problemas não-homogêneos, em que

se podem distinguir regiões no interior das quais há homogeneidade.

Por fim, são apresentados exemplos de aplicação do acoplamento iterativo

acústico-elástico para problemas axissimétricos, conforme explicitado no Capítulo 3,

em que há, no meio em análise, uma região composta de material com

comportamento acústico (fluido, por exemplo), e outra composta de material elástico

(algum material estrutural, como aço ou concreto, por exemplo).

4.1. Meio discretizado exclusivamente pelo MEC

Aqui são apresentadas três aplicações nas quais os conceitos expostos no

Capítulo 2 podem ser empregados.

A primeira aplicação é em uma barra prismática circular, composta de material

acústico, submetida a um fluxo em uma de suas extremidades e com pressão nula na

outra, que na realidade, como será mais explicado adiante, pode ser equiparada a uma

barra composta de material elástico engastada em uma extremidade e submetida a

carregamento normal na outra. Para esta aplicação, são realizadas diversas análises

relativas ao método, que servem então de orientação para as aplicações seguintes.

Para esta mesma barra, em seguida, é feita a análise trocando-se as condições de

contorno das extremidades: onde se prescrevia o fluxo diferente de zero, passa-se a

prescrever pressão, e onde se prescrevia pressão nula, passa-se a prescrever fluxo

Aplicações 73 nulo, o que pode ser equiparado a uma barra elástica com deslocamento prescrito não-

nulo em uma de suas extremidades e livre na outra. Para este caso, apresentado ao

final do item 4.1, é feita uma análise acerca do refinamento da malha empregada na

análise.

A segunda aplicação corresponde a uma cavidade esférica em meio infinito,

submetida a um fluxo ao longo de sua superfície.

A terceira e última aplicação representa uma placa circular de material

acústico que, através da redução de espessura, é aproximada a uma membrana, a qual

se encontra submetida a um fluxo ao longo de uma coroa circular em sua superfície.

Os resultados obtidos pelo programa implementado são comparados com

resultados analíticos conhecidos para estas três aplicações, verificando-se então a

validade do método aqui desenvolvido.

4.1.1. Barra prismática circular submetida a fluxo não-nulo em uma

extremidade

Neste item, é feita a análise através do MEC para meios acústicos

axissimétricos de uma barra prismática circular, constituída de material acústico,

submetida a um fluxo não-nulo em uma de suas extremidades e a pressão nula na

outra.

São analisadas três situações:

− Situação 1: corresponde a uma barra circular vazada, sendo sua seção

transversal uma coroa circular;

− Situação 2: diferencia-se da primeira apenas pelo fato de o cilindro vazado

no interior da barra possuir um raio muito grande, ou seja, a região maciça da barra se

encontra distante do eixo de axissimetria, o que faz com que este caso se aproxime do

caso bidimensional;

− Situação 3: representa uma barra de seção cheia, com pontos do contorno

sobre o eixo de axissimetria.

Através destas três situações, esquematizadas na figura 4.1, chega-se a uma

conclusão acerca da influência da proximidade do eixo de axissimetria em relação aos

Aplicações 74 elementos de contorno na resposta obtida, além de se validar as expressões analíticas

desenvolvidas para todos os casos apresentados no item 2.2.

É apresentada a resposta para pressão em um ponto do contorno na

extremidade da barra onde se aplica o fluxo e em um ponto interno a meio

comprimento da barra. Também se analisa o comportamento do fluxo na extremidade

da barra onde a pressão é nula.

Posteriormente, são feitas considerações acerca do intervalo de tempo de

análise empregado, verificando-se o mais adequado.

Em seguida, a análise é estendida a fim de se verificar a estabilidade do

método. Neste ponto, torna-se necessária a adoção do Esquema Teta apresentado no

item 2.5, estudando-se qual o valor de θ mais adequado para esta situação.

Por fim, faz-se uma comparação entre a análise realizada através da integração

analítica no tempo, desenvolvida neste trabalho, e através de integração numérica ao

longo do tempo, através dos métodos de Gauss-Chebyshev [28] e Kutt [27],

implementada por CZYGAN e VON ESTORFF [24] quando da resolução deste tipo

de problema.

b

ac

q(t)

A

B

C

ac

(a) (b) (c)

Figura 4.1 – (a) esquema do contorno axissimétrico; (b) seção transversal, para situações 1, 2 e 3; (c)

malha empregada na análise

As dimensões da barra são as seguintes: a = 6 m; b = 12 m; c = 1 m (situação

1), 100000 m (situação 2) e 0 m (situação 3). A função q(t) representativa do fluxo

Aplicações 75 prescrito ao longo do tempo é uma função Heaviside unitária, começando em t = 0,

como se fosse uma carga de impacto aplicada à barra neste instante (q(t) = H(t-0)),

medida em kN/m3.

Foram empregados para a discretização do contorno axissimétrico elementos

de 1,5 m, o que produz uma malha de 24 elementos para as situações 1 e 2, e de 16

elementos para a situação 3.

Nos cantos do contorno, são adotados nós duplos, a fim de se representar

corretamente as condições de contorno de cada aresta que conflui para cada canto.

Considera-se para o material representativo do meio acústico uma velocidade

de propagação de onda c = 1,0 m/s.

O intervalo de tempo de análise varia, de acordo com o parâmetro β adotado,

sendo β dado pela expressão (4.1).

l

tcΔ=β (4.1)

Em (4.1), c é a velocidade de propagação da onda primária no meio, Δt o

intervalo de tempo de análise e l o tamanho do elemento de contorno.

Nas figuras 4.2 são traçados os gráficos correspondentes à situação 1, para seis

valores diferentes de β, para os resultados nos pontos A, B e C representados na

figura 4.1. Eles são comparados com os gráficos correspondentes à solução analítica

do problema, apresentada no Anexo F. Nas figuras 4.3, os mesmos gráficos são

traçados para as situações 2 e 3.

O tempo de análise representado nestes gráficos é de 96 s, suficiente para a

vibração provocada pelo fluxo aplicado na extremidade percorrer a barra

longitudinalmente em ida e volta por quatro vezes, produzindo dois ciclos de resposta

em todos os pontos da barra.

Estes ciclos podem ser interpretados fisicamente considerando que, a cada vez

que a onda reflete na extremidade onde a pressão é nula, volta com mesma fase (onda

de compressão volta como onda de compressão e vice-versa), o contrário ocorrendo

na extremidade onde o fluxo é aplicado. Aplica-se um fluxo para dentro do meio

(onda de compressão, derivada em relação à normal negativa), sendo a compressão

representada no gráfico com sinal negativo.

Aplicações 76

0 20 40 60 80 100tempo (s)

-60

-40

-20

0

20

40

60pr

essã

o (k

Pa)

0 20 40 60 80 100tempo (s)

-175

-150

-125

-100

-75

-50

-25

0

25

50

75

100

125

150

175

pres

são

(kPa

)(a) β = 0,1 – pressões em A e B (b) β = 0,1 – fluxo em C

0 20 40 60 80 100tempo (s)

-25

-20

-15

-10

-5

0

5

10

15

pres

são

(kP

a)

0 20 40 60 80 100tempo (s)

-20

-10

0

10

20

pres

ão (k

Pa)

(c) β = 0,2 – pressões em A e B (d) β = 0,2 – fluxo em C

0 20 40 60 80 100tempo (s)

-25

-20

-15

-10

-5

0

5

pres

são

(kP

a)

0 20 40 60 80 100tempo (s)

-1

0

1

2

3

pres

são

(kP

a)

(e) β = 0,5 – pressões em A e B (f) β = 0,5 – fluxo em C

PONTO B

PONTO A

Aplicações 77

0 20 40 60 80 100tempo (s)

-25

-20

-15

-10

-5

0

5

pres

são

(kP

a)

0 20 40 60 80 100tempo (s)

-0.5

0

0.5

1

1.5

2

2.5

pres

são

(kP

a)

(g) β = 0,6 – pressões em A e B (h) β = 0,6 – fluxo em C

0 20 40 60 80 100tempo (s)

-25

-20

-15

-10

-5

0

5

pres

são

(kP

a)

0 20 40 60 80 100tempo (s)

-0.5

0

0.5

1

1.5

2

2.5

pres

são

(kP

a)

(i) β = 0,8 – pressões em A e B (j) β = 0,8 – fluxo em C

0 20 40 60 80 100tempo (s)

-25

-20

-15

-10

-5

0

5

pres

são

(kP

a)

0 20 40 60 80 100tempo (s)

-0.5

0

0.5

1

1.5

2

2.5

pres

são

(kP

a)

(k) β = 1,0 – pressões em A e B (l) β = 1,0 – fluxo em C

Figura 4.2 – comparação entre o resultado da situação 1 (em preto) e o analítico (em vermelho)

Aplicações 78

0 20 40 60 80 100tempo (s)

-25

-20

-15

-10

-5

0

5pr

essã

o (k

Pa)

0 20 40 60 80 100tempo (s)

-8

-4

0

4

8

pres

são

(kP

a)

(a) β = 0,1 – pressões em A e B (b) β = 0,1 – fluxo em C

0 20 40 60 80 100tempo (s)

-25

-20

-15

-10

-5

0

5

pres

são

(kP

a)

0 20 40 60 80 100tempo (s)

-4

-2

0

2

4pr

essã

o (k

Pa)

(c) β = 0,2 – pressões em A e B (d) β = 0,2 – fluxo em C

0 20 40 60 80 100tempo (s)

-25

-20

-15

-10

-5

0

5

pres

são

(kP

a)

0 20 40 60 80 100tempo (s)

-1

0

1

2

3

pres

são

(kP

a)

(e) β = 0,5 – pressões em A e B (f) β = 0,5 – fluxo em C

PONTO BPONTO A

Aplicações 79

0 20 40 60 80 100tempo (s)

-25

-20

-15

-10

-5

0

pres

são

(kP

a)

0 20 40 60 80 100tempo (s)

-0.5

0

0.5

1

1.5

2

2.5

pres

são

(kP

a)

(g) β = 0,6 – pressões em A e B (h) β = 0,6 – fluxo em C

0 20 40 60 80 100tempo (s)

-25

-20

-15

-10

-5

0

5

pres

são

(kP

a)

0 20 40 60 80 100tempo (s)

-0.5

0

0.5

1

1.5

2

2.5

pres

são

(kP

a)

(i) β = 0,8 – pressões em A e B (j) β = 0,8 – fluxo em C

0 20 40 60 80 100tempo (s)

-25

-20

-15

-10

-5

0

5

pres

são

(kP

a)

0 20 40 60 80 100tempo (s)

-0.5

0

0.5

1

1.5

2

2.5

pres

são

(kP

a)

(k) β = 1,0 – pressões em A e B (l) β = 1,0 – fluxo em C

Figura 4.3 – comparação entre o resultado das situações 2 (em preto, linha cheia) e 3 (em preto, linha

pontilhada) e o analítico (em vermelho)

Aplicações 80 Pode-se observar de imediato que para β = 0,1 e β = 0,2, o resultado se torna

instável já para um tempo de análise relativamente pequeno. Isto ocorre pois quanto

menor o valor de β, menor a distância que a onda percorre em um intervalo de tempo.

Sendo assim, o segmento do elemento de contorno a ser integrado concentra em boa

parte de sua extensão os gradientes altos das expressões integradas no tempo,

prejudicando a integração numérica ao longo de seu comprimento. Além disso,

devido à interpolação no espaço, a onda, em determinado instante de tempo, atinge

nós que não deveria, já que a função de interpolação propaga a onda de um nó ao

outro do elemento. Com valores pequenos para β, este aspecto é acentuado, já que são

necessários mais passos de tempo para que de fato a onda atinja o nó em questão,

agravando o problema causado por esta propagação artificial.

À medida que se aumenta o valor de β, este aspecto de instabilidade é

melhorado. Portanto, para um resultado mais preciso, não se pode refinar o intervalo

de tempo da análise arbitrariamente, devendo-se observar o β correspondente a ele.

Ou então, refinam-se concomitantemente o intervalo de tempo e o tamanho do

elemento. Para β = 0,5, percebe-se que, apesar de a solução para pressão ainda não

apresentar acentuada instabilidade, esta deve estar próxima de ocorrer, visto que o

fluxo em C já apresenta fortes oscilações.

Pode-se notar também que sempre o fluxo em C apresenta um aspecto mais

irregular que as pressões em A e B. Provavelmente isso ocorre pela própria natureza

da resposta analítica para o fluxo em C. Esta resposta apresenta saltos em

determinados pontos (derivadas indeterminadas), o que torna mais difícil sua

aproximação numérica.

É importante ressaltar também que os resultados para as situações 2 e 3 se

mostram mais bem comportados do que aqueles para a situação 1. Este aspecto é um

indício de que a integração pode estar perdendo precisão e, conseqüentemente,

prejudicando a qualidade da solução para contornos mais próximos do eixo de

axissimetria, de acordo com as expressões analíticas obtidas através da integração ao

longo do tempo.

Deve-se lembrar que a solução axissimétrica, para a situação 2, em que o

domínio axissimétrico encontra-se muito afastado do eixo de axissimetria, acaba

recaindo na solução bidimensional, pelo menos durante o tempo em que não há

influência, sobre este contorno, dos pontos mais distantes pertencentes ao anel de

Aplicações 81 fontes axissimétrico. Foram comparados os resultados desta situação com resultados

obtidos por programas bidimensionais já elaborados e as respostas foram idênticas,

até a precisão adotada.

O próximo passo foi estender a análise até que a solução se instabilize,

acompanhando o comportamento da solução em função do β adotado na análise. Nas

figuras 4.4, são apresentados os resultados para a pressão no ponto A nas situações 1 e

2, para os três maiores valores de β apresentados anteriormente, já que foi verificado

que para a situação 3 o comportamento é semelhante ao observado para a situação 2, e

para β menor a solução já se apresentava bastante deteriorada para instantes de tempo

menores. O tempo máximo de análise é de 576 s (12 ciclos de resposta).

Aplicações 82

0 200 400 600tempo (s)

-25

-20

-15

-10

-5

0

pres

são

(kP

a)

0 100 200 300 400tempo (s)

-40

-30

-20

-10

0

10

pres

são

(kP

a)

(a) β = 0,6 – situação 2 (b) β = 0,6 – situação 1

0 200 400 600tempo (s)

-25

-20

-15

-10

-5

0

pres

são

(kP

a)

0 200 400 600tempo (s)

-25

-20

-15

-10

-5

0

pres

são

(kP

a)

(c) β = 0,8 – situação 2 (d) β = 0,8 – situação 1

0 200 400 600tempo (s)

-25

-20

-15

-10

-5

0

pres

são

(kP

a)

0 50 100 150 200 250tempo (s)

-30

-20

-10

0

10

pres

são

(kP

a)

(e) β = 1,0 – situação 2 (f) β = 1,0 – situação 1

Figura 4.4 – pressão no ponto A para análise estendida

Para a situação 2 (contorno afastado do eixo) a solução se estabiliza para todos

os β tomados acima de 0,6. Já para a situação 1, apenas para β = 0,8, a solução se

mostra estável até o tempo final de análise. Além da proximidade do contorno,

Aplicações 83 provavelmente as expressões obtidas para o caso 3 e 5 a que se refere o item 2.2

possuem um comportamento mais irregular em relação às demais. Na situação 2,

durante todo o tempo de análise estendido (576 s), estes casos nunca aparecem, em

virtude do afastamento em relação ao eixo de axissimetria.

Nas figuras 4.4, pode-se notar o amortecimento considerável presente na

solução numérica à medida que a análise avança no tempo. Este amortecimento é

inerente ao MEC.

Na figura 4.5, é feita uma comparação entre os amortecimentos, para a

situação 2, de acordo com o β adotado. Por ela, deduz-se que este amortecimento é

maior para valores menores de β.

0 200 400 600tempo (s)

-25

-20

-15

-10

-5

0

pres

são

(kP

a)

beta = 0,6beta = 0,8beta = 1,0

Figura 4.5 – pressão no ponto A para análise estendida, situação 2 – comparação entre amortecimentos

numéricos

Por esta conclusão, teoricamente o β mais adequado seria aquele igual à

unidade. No entanto, como se pode inferir das figuras 4.4, este valor pode apresentar

instabilidade dependendo do problema. Deste modo, a escolha do β a ser empregado

na análise depende do problema em questão.

Para se tentar eliminar as instabilidades presentes, uma alternativa é a adoção

do Esquema Teta, apresentado por MANSUR et al. [32] e resumidamente explicado

no item 2.5. Foi então executada a situação 1, que se mostra a mais problemática em

termos de instabilidade, para os valores de β iguais a 0,6 e 0,8, para três valores de θ

(1,05; 1,10 e 1,20), a fim de se verificar o comportamento do amortecimento

Aplicações 84 numérico com a variação de θ, sendo os resultados para pressão no ponto A

apresentados na figura 4.6.

0 200 400 600tempo (s)

-25

-20

-15

-10

-5

0

pres

são

(kP

a)

0 200 400 600tempo (s)

-25

-20

-15

-10

-5

0

pres

são

(kP

a)

(a) β = 0,6 - θ = 1,05 (b) β = 0,8 - θ = 1,05

0 200 400 600tempo (s)

-25

-20

-15

-10

-5

0

pres

são

(kP

a)

0 200 400 600tempo (s)

-25

-20

-15

-10

-5

0

pres

são

(kP

a)

(c) β = 0,6 - θ = 1,10 (d) β = 0,8 - θ = 1,10

0 200 400 600tempo (s)

-25

-20

-15

-10

-5

0

pres

são

(kP

a)

0 200 400 600tempo (s)

-25

-20

-15

-10

-5

0

pres

são

(kP

a)

(e) β = 0,6 - θ = 1,20 (f) β = 0,8 - θ = 1,20

Figura 4.6 – variação da resposta para pressão no ponto A com a mudança de θ

Aplicações 85

À medida que se aumenta o valor de θ, aumenta também o amortecimento

numérico experimentado pela solução, ou seja, é como se a estabilização se desse às

custas deste amortecimento. Pode-se notar também que com o aumento de θ o

amortecimento para β = 0,8 é maior que para β = 0,6, ao contrário do que acontece

para θ = 1,0 (esquema padrão, figura 4.5).

A última análise relativa a esta aplicação corresponde à comparação entre o

procedimento de integração analítica no tempo, conforme descrito no item 2.2, e o

procedimento numérico adotado em [24], para a situação 1, com θ = 1,10 e β = 0,6.

Para implementação numérica destas integrais são empregados dois métodos

distintos, de acordo com CZYGAN e VON ESTORFF [24]:

Para as integrais necessárias à montagem das matrizes G (expressão 2.12)

é empregado o método de Gauss-Chebyshev, sobre o qual maiores detalhes podem ser

encontrados em ABRAMOWITZ e STEGUN [28]. São tomados sempre 20 pontos

para o cálculo da integral.

Para as integrais necessárias à montagem das matrizes H (expressão 2.12)

é empregado o método de Kutt, sobre o qual maiores detalhes podem ser encontrados

em KUTT [27]. Este método é especialmente utilizado quando do cálculo de integrais

que só existem no sentido de parte finita, conforme exposto no Anexo C. Para isso,

são utilizados abscissas e pesos específicos, de acordo com o método. Neste trabalho,

tomam-se sempre 12 pontos de Kutt para o cálculo das integrais.

Na figura 4.7, é feita uma comparação entre os resultados obtidos para a

pressão no ponto A, com θ = 1,10 e β = 0,6, através das expressões analíticas das

integrais no tempo e por integração numérica. Pela análise desta figura, pode-se notar

que à medida que a análise avança no tempo, a solução obtida calculando-se as

integrais no tempo numericamente apresenta uma defasagem de período e redução de

amplitude em relação à solução obtida empregando-se as expressões analíticas para

estas mesmas integrais.

Através de testes realizados, notou-se que a principal fonte de erro na

integração numérica é a integração através do método de Gauss-Chebyshev utilizada

para se efetuar as integrais responsáveis pela montagem da matriz G.

Por esta análise, pode-se concluir, como já era esperado, que o cálculo

analítico das integrais ao longo do tempo apresenta resultados melhores que a

Aplicações 86 integração numérica, melhoria esta que se torna considerável à medida que se estende

o tempo de análise.

0 200 400 600tempo (s)

-25

-20

-15

-10

-5

0

pres

são

(kP

a)

integração numéricaintegração analíticasolução analítica

Figura 4.7 – resposta para pressão no ponto A para integração no tempo numérica e analítica

É bom ressaltar que esta aplicação apresentada no item 4.1 é um tanto quanto

problemática, visto que nela ocorrem sucessivas reflexões nas extremidades da barra.

Essas múltiplas reflexões podem acentuar alguns aspectos relativos ao método, como

por exemplo a instabilidade verificada em alguns casos e, por outro lado, o

amortecimento artificial excessivo verificado em outros.

Aplicações 87

4.1.2. Cavidade em meio acústico infinito

Neste item, aplica-se o método desenvolvido a uma cavidade esférica inserida

em um meio acústico infinito, com velocidade de propagação c = 1000 m/s. O modelo

é apresentado na figura 4.8.

Esta cavidade é submetida a um fluxo q(t) em sua superfície, descrito ao longo

do tempo por uma função Heaviside unitária, em kN/m3 (q(t) = H(t-0)).

O raio da cavidade é tomado como a = 1 m.

Para a análise, são apresentados os resultados para pressão ao longo do tempo

tanto em pontos do contorno quanto em pontos internos ao domínio, representados na

figura 4.8.

O contorno axissimétrico foi discretizado inicialmente por oito elementos. Foi

escolhido um intervalo de tempo de análise de 0,2344 ms, o que fornece um β

aproximadamente igual a 0,6. O tempo total de análise empregado corresponde a 64

passos de tempo.

Primeiramente, o modelo é executado com θ = 1,0.

a

2 3 4 5 61

78

(a) (b)

Figura 4.8 – (a) esquema da cavidade esférica em meio infinito; (b) malha de oito elementos

empregada

Aplicações 88

0 0.004 0.008 0.012 0.016tempo (s)

0

-0.4

-0.8

-1.2

pres

são

(kP

a)

ponto 8ponto 7ponto 1analítico

Figura 4.9 – resultados para 8 elementos e θ = 1,0 – pontos do contorno

Na figura 4.9 são traçados os gráficos para a pressão em pontos pertencentes

ao contorno (ponto 1, no eixo horizontal, ponto 7, a 45 graus, e ponto 8, no eixo de

axissimetria). Os resultados para os pontos 1 e 7 aproximam-se bem mais do analítico

em relação ao ponto 8.

A análise não apresenta sinais de instabilidade, mesmo quando o tempo total é

aumentado. Este aspecto já era esperado, visto que se trata de um meio infinito, no

qual não ocorrem reflexões. Deste modo, a energia se espalha e não há acúmulo de

erros de precisão das integrações efetuadas.

Pela figura 4.10, pode-se notar que ao se empregar o Esquema Teta com θ =

1,1 o resultado praticamente não se altera, o que mostra que esta aplicação não

apresenta foco significativo de instabilidade.

Na figura 4.11, é mostrada a comparação entre resultados para a malha de 8

elementos descrita anteriormente e para uma malha com 16 elementos de contorno.

Nesta malha, o intervalo de tempo adotado foi de 0,1172 ms, mantendo o valor de β

bem próximo a 0,6.

É nítido que a malha de 16 elementos apresenta resultados com menos

oscilações. A precisão obtida para o ponto 8 melhorou com o refinamento. Os

resultados para este ponto, localizado sobre o eixo de axissimetria, parecem bastante

sensíveis aos "bicos" formados pela discretização com elementos retos. Uma

Aplicações 89 alternativa a ser testada é o emprego de elementos quadráticos para discretização da

geometria, pois neste caso, possivelmente, bons resultados para pontos do tipo do

ponto 8 podem ser obtidos sem que haja necessidade de refinamento excessivo.

0 0.004 0.008 0.012 0.016tempo (s)

0

-0.2

-0.4

-0.6

-0.8

-1pr

essã

o (k

Pa)

teta = 1,0teta = 1,1

Figura 4.10 – resultados para o ponto 1, variando-se o valor de θ

0 0.004 0.008 0.012 0.016tempo (s)

0

-0.4

-0.8

-1.2

pres

são

(kP

a)

ponto 8 - 8 elementosponto 1 - 8 elementosanalíticoponto 8 - 16 elementosponto 1 - 16 elementos

Figura 4.11 – comparação entre malhas

PONTO 1 PONTO 8

Aplicações 90 Por fim, na figura 4.12, são traçados gráficos para pressão ao longo do tempo

para os pontos internos com uma malha de 16 elementos, comparando os resultados

obtidos com os respectivos resultados analíticos.

Como se pode observar, os resultados apresentam precisão bastante aceitável.

Figura 4.12 – pontos internos (em preto), comparação com resultado analítico (em vermelho)

0 0.004 0.008 0.012 0.016tempo (s)

0

-0.2

-0.4

-0.6

-0.8

pres

são

(kP

a)

PONTO 2

PONTO 3

PONTO 4

PONTO 5

PONTO 6

Aplicações 91

4.1.3. Membrana submetida a carga impulsiva

Esta aplicação consiste de uma placa circular submetida a um carregamento

impulsivo, ou seja, com duração bem pequena, em uma tentativa de se simular

condições iniciais de velocidade. O carregamento é aplicado na face superior da placa

de duas formas: na situação 1, ele é aplicado em um círculo com centro coincidente

com o da placa; na situação 2, ele é aplicado em uma coroa circular também tendo

como centro o próprio centro da placa. Neste problema, ao invés de se utilizar a

pressão como variável básica, utiliza-se o deslocamento vertical u. Esta consideração

é feita desprezando-se o deslocamento horizontal, considerado bem menor em

comparação com o vertical. Deste modo, resolvendo-se a equação da acústica, o

resultado de pressão é igual ao deslocamento transversal da membrana.

A derivada de u em relação à normal ao contorno, multiplicada pela constante

K (módulo de compressibilidade do meio), dada por (4.2), representa carga por

unidade de área normal ao contorno. Portanto, quando neste problema prescreve-se a

derivada em relação à normal ao contorno, indiretamente se está prescrevendo carga

aplicada. Em (4.2), ρ representa a densidade do meio e c, a velocidade de propagação

da onda de compressão no mesmo.

2cK ρ= (4.2)

Ressalta-se que, para haver correspondência de amplitudes entre carga

aplicada e condições iniciais de velocidade, devem-se igualar as quantidades de

movimento produzidas em ambos os casos, chegando-se à seguinte relação:

tfcev Δ= 2

0 (4.3)

Em (4.3), e é a espessura da placa, c a velocidade de propagação da onda

acústica no meio, v0 é a velocidade inicial, suposta constante ao longo da região onde

é aplicada, f a derivada em relação à normal prescrita, constante ao longo da região de

aplicação e ao longo do tempo, e Δt o intervalo de aplicação de f.

A espessura é sucessivamente reduzida em relação ao raio da placa na

tentativa de se chegar, na prática, a uma situação de membrana, para a qual os

Aplicações 92 resultados analíticos de deslocamento são conhecidos, no caso de velocidade inicial, e

apresentados no Anexo F.

É importante ressaltar que quando se refere a membrana, está se tratando de

um meio com espessura desprezível em relação às demais dimensões. Daí a opção por

se tomar espessuras pequenas para o modelo de placa acústica adotado. A esta altura,

poderia se questionar por que então a espessura não é de imediato reduzida a um valor

muito pequeno. Isto não é feito porque assim os elementos na face vertical do

contorno axissimétrico da figura 4.13 seriam muito pequenos. Como não é

recomendável no MEC a presença de elementos de contorno adjacentes de tamanhos

muito diferentes (o maior deve sempre ser menor que duas vezes o menor), os

elementos na horizontal deveriam ser também bem reduzidos, aumentando

excessivamente o número de elementos.

f(t)

f(t)

a

b

b

a

(a) (b)

e

r

r

Figura 4.13 – esquema do modelo da membrana: (a) em perfil; (b) em planta.

f(t)

tdt

1

Figura 4.14 – gráfico no tempo da carga aplicada

Na situação 1, toma-se a = 0 m; b = 0,5 m; enquanto para a situação 2, a = 1,5

m e b = 2,0 m. O raio da placa é sempre mantido como r = 4,0 m. A carga f(t) é uma

Aplicações 93 função pulso descrita pela figura 4.14, sendo o tempo de aplicação tomado como igual

a 1 intervalo de tempo de análise.

O material da placa possui velocidade de propagação c = 1 m/s. Ela é

discretizada por elementos de contorno de comprimento igual a 0,05 m, sendo o

número de elementos variável com a espessura. Mantém-se sempre β = 0,6, sendo o

intervalo de tempo, conseqüentemente, igual a 0,03 s. Emprega-se θ = 1,0, já que,

conforme inferido das figuras que contém gráficos de resposta do problema, para os

tempos totais de análise aqui empregados não ocorrem instabilidades.

Inicialmente, são feitos testes com diversas espessuras, a fim de se encontrar

uma suficientemente pequena para que o problema possa ser equiparado ao da

vibração transversal de uma membrana. Na figura 4.15, faz-se uma comparação entre

resultados de deslocamento no centro da placa obtidos para espessuras de 0,1 m, 0,2

m, 1,0 m, e o resultado analítico para membrana, na situação 2, para um tempo total

de 20 s, suficiente para que cheguem, no ponto central da membrana, duas reflexões

no contorno.

Conforme esperado, à medida que se reduz a espessura, a resposta se aproxima

da analítica. Considera-se que os resultados para uma espessura de 0,1 m, que

corresponde a uma relação entre espessura e raio de 1/40, são suficientemente

próximos da solução analítica.

0 4 8 12 16 20tempo (s)

-0.6

-0.4

-0.2

0

0.2

0.4

desl

ocam

ento

(mm

)

e = 0,20 me = 0,10 mmembrana analítica

0 4 8 12 16 20tempo (s)

-0.08

-0.04

0

0.04

0.08

0.12

0.16

desl

ocam

ento

(mm

)

e = 1,0 m

(a) (b)

Figura 4.15 – comparação de espessuras: (a) e = 0,1m; 0,2 m e teórico; (b) e = 1,0 m

Aplicações 94

Para e = 1,0 m (figura 4.15(b)), o resultado diverge bastante do analítico. Para

tal espessura, passam a ser muito significativas, em relação às que ocorrem no

contorno externo, as reflexões que ocorrem na face inferior da placa e atingem o

ponto em questão, na face superior, invalidando qualquer consideração de

comportamento de membrana.

Na figura 4.16, é feita a comparação entre a resposta para deslocamento no

centro da placa para a situação 1, com espessura de 0,10 m, e a resposta analítica para

a membrana circular na mesma situação.

0 4 8 12 16 20tempo (s)

-0.2

-0.1

0

0.1

0.2

desl

ocam

ento

(mm

)

placa numéricamembrana analítica

Figura 4.16 – comparação entre resultado numérico e analítico – situação 1

Pode-se notar que, para a situação 1, apesar de algumas oscilações em tempos

menores de análise, o resultado, à medida que o tempo cresce, aproxima-se mais do

analítico do que para a situação 2. Uma possível causa para isso é a maior ocorrência

de reflexões na segunda situação, visto que a coroa circular onde é aplicado fluxo

emite ondas para dentro e para fora de si mesma, provocando ondas "duplas" em

relação às observadas para a situação 1. Este número maior de reflexões acaba

provocando maiores alterações na resposta à medida que o tempo de análise cresce.

Pode-se considerar, no entanto, em vista da diferença dos problemas

aplicados, que os resultados são bastante satisfatórios, já que apresentam sempre

variação bem semelhante.

Aplicações 95

Deve-se ressaltar que, para que houvesse equivalência entre resultados para

velocidade inicial e carga impulsiva aplicadas, foi empregada a equação (4.3). Deste

modo, os resultados acima apresentados foram obtidos, quando analíticos, para uma

velocidade inicial v0 = 1000 m/s e, quando numéricos, para uma derivada em relação

ao contorno de 3,3333.103. Como K = 1,0 kN/m2, a carga aplicada foi de 3,3333

kN/m2.

Nas figuras 4.17 e 4.18, são apresentados os resultados para as situações 1 e 2

relativos à derivada no contorno direito da figura 4.13(a), que corresponde,

multiplicada por K, à reação vertical de apoio na membrana.

0 4 8 12 16 20tempo (s)

-0.12

-0.08

-0.04

0

0.04

0.08

0.12

∂u/∂

n

Figura 4.17 – derivada normal ao contorno direito – situação 1

0 4 8 12 16 20tempo (s)

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

∂u/∂

n

Figura 4.18 – derivada normal ao contorno direito – situação 2

Aplicações 96

0 20 40 60 80 100tempo (s)

-2.5

-2

-1.5

-1

-0.5

0

0.5

pres

são

(Pa)

menos refinadomais refinadoanalítico

4.1.4. Barra prismática circular submetida a pressão não-nula em uma

extremidade

Este exemplo corresponde a um meio idêntico àquele apresentado no item

4.1.1, figura 4.1, situação 1. No entanto, neste caso, as condições de contorno nas

extremidades da barra são trocadas. Onde se prescrevia fluxo não-nulo, passa a se

prescrever pressão não-nula, enquanto que, onde se prescrevia pressão nula, passa a se

prescrever fluxo nulo. Procedendo-se desta maneira, pode-se verificar se a utilização

de interpolação constante no tempo para o fluxo no contorno é adequada, visto que

neste caso estes valores não são prescritos, mas calculados pelo método. São

empregadas duas malhas: primeiramente o meio é discretizado através da mesma

malha utilizada em 4.1.1; posteriormente, o tamanho dos elementos é reduzido à

metade, produzindo uma malha com 60 elementos de contorno de 0,75 m.

Emprega-se sempre θ = 1,1, pois como foi analisado em 4.1.1, para a barra

com contorno próximo ao eixo de axissimetria, o problema apresenta instabilidade,

provavelmente devida às expressões presentes nas integrais ao longo do contorno.

Para manter β = 0,6 em todos os casos, utiliza-se um intervalo de tempo (dt)

igual a 0,9 s para a malha menos refinada e igual a 0,45 s para a mais refinada.

Na figura 4.19 são apresentados os resultados para pressão em um ponto

localizado no contorno a meio comprimento da barra, para as duas discretizações,

bem como o resultado analítico.

Figura 4.19 – pressão em um ponto a meia barra

Aplicações 97

Pela figura 4.19, nota-se uma sensível melhora na resposta com o refinamento

realizado na malha de elementos de contorno, com a redução das oscilações e a

conseqüente aproximação da resposta analítica.

Na figura 4.20, está apresentado o fluxo no contorno em um ponto na

extremidade da barra na qual se aplica a pressão. Os picos ocorrem a cada passagem

da frente de onda por essa extremidade e a conseqüente reflexão. Percebe-se que em

cada pico ocorre uma mudança de sinal, que significa uma mudança de fase. Esta

mudança (passagem de onda de compressão para onda de tração) é a responsável pelo

caráter oscilatório da resposta de pressão apresentada na figura 4.19.

0 20 40 60 80 100tempo (s)

-3

-2

-1

0

1

2

3

fluxo

Figura 4.20 – fluxo na extremidade da barra

Os resultados obtidos estão dentro do padrão esperado, validando portanto a

utilização das funções de interpolação no tempo constantes para o fluxo.

Aplicações 98

4.2. Acoplamento acústico-acústico

Neste item é apresentado apenas um exemplo de aplicação relativa ao

acoplamento acústico-acústico MEC-MEF, de modo a validar o método aqui utilizado

para promover este acoplamento.

Esse exemplo consiste de uma barra prismática de seção transversal circular

vazada, semelhante àquela apresentada na situação 1 do item 4.1.1. A diferença é que

aqui esta barra possui cada metade de seu comprimento constituída por um material

distinto, sendo ambos os materiais acústicos.

A metade superior, de acordo com a figura 4.21, é modelada pelo MEF, sendo

o restante modelado pelo MEC.

Os materiais possuem as seguintes propriedades:

Meio 1: ρ = 1,0 kg/m3; c = 1,0 m/s;

Meio 2: ρ = 1,0 kg/m3; c = 2,0 m/s.

Na extremidade superior é aplicado um fluxo descrito por uma função

Heaviside unitária ao longo do tempo – q(t) = H(t-0). A extremidade inferior

apresenta como condição de contorno pressão nula.

12

61

q(t)

A

B

MEIO 1 - MEC

MEIO 2 - MEF

p = 0

Figura 4.21 – modelo para acoplamento acústico-acústico

Aplicações 99

A malha do MEF é discretizada por elementos quadrados bilineares de lado

igual a 0,75 m, sendo empregado um passo de tempo igual a 0,15 s; enquanto isso, a

malha do MEC é discretizada por elementos lineares de comprimento igual a 0,75 m,

adotando-se um passo de tempo igual a 0,60 s, o que equivale a um parâmetro β = 0,8.

Neste ponto observa-se uma vantagem do acoplamento iterativo. O emprego

do MEC é vantajoso para meios homogêneos, sendo que no caso de meios não-

homogêneos, havendo, em seu interior, regiões homogêneas distintas, pode-se fazer

uso da partição em sub-regiões. O acoplamento iterativo permite que cada região seja

discretizada independentemente, por métodos distintos ou não, procedendo-se ao

acoplamento apenas na interface, iterativamente, conforme exposto no Capítulo 3.

Para efeito de comparação é utilizado um programa que emprega o Método

das Diferenças Finitas, baseado no operador heterogêneo desenvolvido por COHEN e

JOLY [34], exaustivamente testado e de eficiência comprovada.

São tomadas como referência as pressões ao longo do tempo obtidas para os

pontos A e B da figura 4.21.

A análise das figuras 4.22 e 4.23 mostra que os resultados obtidos através do

programa aqui desenvolvido para acoplamento acústico-acústico encontram-se

satisfatoriamente próximos daqueles obtidos pelo programa de MDF, validando a

implementação computacional realizada.

0 20 40 60 80 100tempo (s)

-20000

-16000

-12000

-8000

-4000

0

pres

são

(Pa)

acoplamento MEC-MEFMDF

Figura 4.22 – pressão no ponto A – acoplamento x MDF

Aplicações 100

0 20 40 60 80 100tempo (s)

-16000

-12000

-8000

-4000

0

4000

pres

são

(Pa)

acoplamento MEC-MEFMDF

Figura 4.23 – pressão no ponto B – acoplamento x MDF

Aplicações 101

4.3. Acoplamento acústico-elástico

São apresentados na seqüência três exemplos de aplicação do método aqui

desenvolvido para análise transiente de um sistema com acoplamento acústico-

elástico. Através deles, é possível a confirmação da validade dos conceitos discutidos

no Capítulo 3.

4.3.1. Barra prismática circular vazada

Neste item, a aplicação corresponde a uma barra prismática, cuja seção

transversal é uma coroa circular, idêntica àquela apresentada na situação 1 do item

4.1.1.

A diferença é que, aqui, de uma extremidade da barra até metade de seu

comprimento, a barra é constituída de um meio acústico, enquanto que o seu restante

é composto de um meio elástico.

São abordadas para esta aplicação duas situações, conforme indica a figura

4.24, na qual o meio 1 representa o meio acústico e o meio 2 o elástico:

− Situação 1: carga descrita pela função Heaviside (f(t) = H(t-0)) aplicada na

extremidade do meio elástico e fluxo nulo prescrito na extremidade do meio acústico;

− Situação 2: pressão descrita pela função Heaviside (p(t) = H(t-0)) aplicada

na extremidade do meio acústico e deslocamento na direção do eixo da barra nulo

prescrito na extremidade do meio elástico.

Adotam-se para os materiais constituintes dos meios as seguintes propriedades

(para simplificar foram escolhidos sempre valores unitários):

Meio 1: ρ = 1,00 kg/m3; c = 1,00 m/s;

Meio 2: ρ = 1,00 kg/m3; E = 1,00 N/m2; ν = 0,0. O valor nulo para o

coeficiente de Poisson é adotado para igualar a velocidade propagação da onda

primária deste meio à do meio acústico.

A malha de elementos de contorno através da qual é discretizado o meio 1 é

composta de elementos com comprimento igual a 0,75 m. Para a malha de elementos

finitos, através da qual se discretiza o meio 2, inicialmente são empregados elementos

quadrados de lado igual a 0,75 m.

O passo de tempo da análise no MEC é de 0,45 s, o que fornece um β = 0,6.

Para o MEF, esse valor também é adotado de início. O tempo total de análise é de 360

Aplicações 102 s. Na malha modelada pelo MEC, adota-se θ = 1,1, para reduzir possíveis

instabilidades.

São apresentados na seqüência alguns testes efetuados. Estes testes dizem

respeito a quatro parâmetros da análise: o valor da tolerância empregado na análise

para a verificação da convergência; o valor de α utilizado para relaxação ao final das

iterações; o refinamento da malha de elementos finitos na direção longitudinal à barra;

e o intervalo de tempo de análise adotado nesta malha. Para a realização destes testes,

é empregado o modelo da situação 1.

Com estes testes, apesar da simplicidade da aplicação, são esclarecidos

aspectos importantes relativos ao acoplamento, que se tornam muito úteis quando da

execução de outras aplicações mais complexas.

p(t)

A

B

MEIO 1

MEIO 2

12

f(t)

z

r

C

D

A

B

C

D

MEIO 1

MEIO 2

uz = 0

q = 0

12

MEC

MEF

(a) (b) (c)

Figura 4.24 – esquema da aplicação do item 4.3.1: (a) situação 1; (b) situação 2; (c) malha empregada;

cotas em metros.

De início, adota-se o parâmetro α igual a 1, ou seja, não é efetuada relaxação

ao final das iterações para determinado passo de tempo. Deste modo, tomam-se para

as cargas nodais na interface da malha de elementos finitos os valores obtidos através

de elementos de contorno ao final da iteração anterior.

Primeiramente, adota-se uma tolerância ε igual a 10-3 para verificação da

convergência. Conforme ilustrado na figura 4.25, no entanto, ela se mostra deficiente.

Quando se reduz a tolerância para 10-6, nota-se que o resultado, que tendia a não mais

oscilar em relação a uma reta horizontal à medida que a análise se estendia, passa a se

estabilizar neste sentido.

Aplicações 103

No entanto, podem-se notar ainda grandes instabilidades ocorrendo em

instantes de tempo mais elevados. O procedimento seguinte, então, é a adoção de um

parâmetro de relaxação de final de iteração α igual a 0,5. O resultado, no entanto,

apresenta-se muito semelhante àquele obtido na figura 4.25, não sendo exposto aqui.

0 100 200 300 400tempo (s)

-60

-40

-20

0

desl

ocam

ento

(m)

tol = 10e-3tol = 10e-6

0 100 200 300 400tempo (s)

-50

-40

-30

-20

-10

0

10

desl

ocam

ento

(m)

tol = 10e-3tol = 10e-6

(a) (b) Figura 4.25 – deslocamentos uz em A (a) e B (b) para tolerâncias de 10-3 e 10-6

O passo de tempo empregado no meio elástico discretizado pelo MEF é então

reduzido, primeiramente a 0,225 s e, em seguida, a 0,1125 s, o que corresponde,

respectivamente, a razões ΔtC/ΔtF de 2 e 4.

0 100 200 300 400tempo (s)

-25

-20

-15

-10

-5

0

desl

ocam

ento

(m)

dtc/dtf = 2dtc/dtf = 4

0 100 200 300 400tempo (s)

-12

-10

-8

-6

-4

-2

0

2

desl

ocam

ento

(m)

dtc/dtf = 2dtc/dtf = 4

(a) (b) Figura 4.26 – deslocamentos uz em A (a) e B (b) variando passo de tempo do MEF

Aplicações 104

Pode-se notar que com a redução do passo de tempo de análise do meio

modelado pelo MEF, o problema estabilizou-se. No entanto, assim como quando o

meio é modelado apenas pelo MEC (item 4.1.1), ocorre ao longo do tempo um

amortecimento considerável na amplitude de deslocamentos.

Além disso, observam-se desvios de trajetória no que deveriam ser retas, nos

gráficos relativos ao ponto A, correspondentes a irregularidades nos patamares nos

gráficos relativos ao ponto B.

Apesar de não muito visível na figura 4.26(a), quando se empregou ΔtC/ΔtF =

= 4, a resposta apresentou um patamar mais definido do que quando esta razão foi

igual a 2. Portanto, adota-se como referência o valor 4.

É importante ressaltar também que, apesar de a utilização do parâmetro α não

ter feito diferença quando se empregou intervalos de tempo iguais para os meios

modelados pelo MEC e pelo MEF, aqui ela assume grande importância. Quando este

parâmetro não é utilizado, o número de iterações necessárias cresce vertiginosamente,

essencialmente nos instantes de tempo de resolução de ambos os meios.

Por último, a malha de elementos finitos é refinada na direção longitudinal à

barra, tendo sua dimensão nesta direção reduzida à metade (de 0,75 para 0,375 m). Os

resultados apresentados na figura 4.27 mostram que o amortecimento numérico da

resposta foi reduzido, indicando uma considerável melhora. Também, apesar de não

muito visível na figura 4.27(b), os desvios de trajetória são razoavelmente reduzidos.

0 100 200 300 400tempo (s)

-25

-20

-15

-10

-5

0

desl

ocam

ento

(m)

MEF malha inicialMEF malha refinada

0 100 200 300 400tempo (s)

-12

-10

-8

-6

-4

-2

0

2

desl

ocam

ento

(m)

MEF malha inicialMEF malha refinada

(a) (b) Figura 4.27 – deslocamentos uz nos pontos A (a) e B (b) variando refinamento da malha do MEF

Aplicações 105

Para o modelo final adotado (ε = 10-6; α =0,5; ΔtC/ΔtF = 4 e malha do MEF

mais refinada na direção longitudinal à barra), então, é feita uma comparação entre os

resultados obtidos com o programa implementado e a solução analítica, para um

tempo total de 96 s. Para isso, são tomados os deslocamentos uz nos pontos A e B, e

as pressões nos pontos B e C.

0 20 40 60 80 100tempo (s)

-25

-20

-15

-10

-5

0

desl

ocam

ento

(m)

analíticonumérico

0 20 40 60 80 100tempo (s)

-14

-12

-10

-8

-6

-4

-2

0

2

desl

ocam

ento

(m)

analíticonumérico

(a) deslocamento uz – ponto A (b) deslocamento uz – ponto B

0 20 40 60 80 100tempo (s)

-2.5

-2

-1.5

-1

-0.5

0

0.5

pres

são

(Pa)

analíticonumérico

0 20 40 60 80 100tempo (s)

-2.5

-2

-1.5

-1

-0.5

0

0.5

pres

são

(Pa)

analíticonumérico

(c) pressão – ponto B (d) pressão – ponto C

Figura 4.28 – comparação entre resultados analíticos e numéricos - situação 1

Por fim, são mostrados os resultados obtidos para a situação 2, na figura 4.29,

de deslocamentos uz nos pontos B e D e pressões nos pontos B e C.

Aplicações 106

0 20 40 60 80 100tempo (s)

-2

0

2

4

6

8

10

12

14

desl

ocam

ento

(m)

analíticonumérico

0 20 40 60 80 100tempo (s)

-1

0

1

2

3

4

5

6

7

desl

ocam

ento

(m)

analíticonumérico

(a) deslocamento uz – ponto B (b) deslocamento uz – ponto D

0 20 40 60 80 100tempo (s)

-2.5

-2

-1.5

-1

-0.5

0

0.5

pres

são

(Pa)

analíticonumérico

0 20 40 60 80 100tempo (s)

-2.5

-2

-1.5

-1

-0.5

0

0.5

pres

são

(Pa)

analíticonumérico

(c) pressão – ponto B (d) pressão – ponto C

Figura 4.29 – comparação entre resultados analíticos e numéricos - situação 2

Nota-se que os deslocamentos, ao invés de sofrerem um amortecimento, como

na situação 1, neste caso, atingem amplitudes maiores para o segundo ciclo da

resposta, além de não apresentarem um patamar bem definido. Estendendo-se a

análise, esse comportamento vai se acentuando.

Aplicações 107

4.3.2. Fonte emissora no interior de um duto metálico

Esta aplicação corresponde a um meio acústico, constituído de um fluido, aqui

considerado água, no qual está imerso um duto metálico, composto de aço. No eixo

do duto posiciona-se uma fonte geradora de pressão, desejando-se obter a resposta

para a pressão ao longo do tempo em um ponto externo ao duto. Este ponto, em um

caso real, pode representar um receptor de sinais acústicos, que aciona determinado

equipamento, pela transformação destes sinais em impulsos elétricos.

A configuração do modelo é indicada na figura 4.30. A sua altura é adotada de

modo a, durante o tempo total de análise, não haver influência de reflexões nos bordos

superior e inferior.

A BFONTE

AÇO

ÁGUAÁGUA

C

detalhe da malha do MEC

detalhe da malha do MEF

Figura 4.30 – modelo de um meio afetado por uma fonte geradora de pressão; cotas em metros

São adotadas as seguintes características para os meios:

Água: ρ = 1000 kg/m3; c = 1500 m/s;

Aço: ρ = 7700 kg/m3; E = 2,1.1011 N/m2; ν = 0,30.

Estes valores fornecem para o aço uma velocidade de propagação da onda de

compressão igual a 6059 m/s, de acordo com a expressão (4.4).

( )( )( )ννρ

ν211

11 −+

−=

Ec (4.4)

Aplicações 108

O sinal emitido pela fonte é constituído de três ciclos de senos, de amplitude

unitária, em MPa, com freqüência igual a 10000 Hz.

Estende-se a análise por um tempo total de 0,002 s, que já é suficiente para se

perceberem várias reflexões, nos pontos onde se mede a resposta, advindas da parede

interna do duto.

O contorno do duto é discretizado através de elementos de contorno de 0,02 m

de comprimento, que representam uma cavidade no meio infinito. O interior do duto é

discretizado por elementos finitos quadrados, de lado também igual a 0,02 m.

O passo de tempo da análise por elementos de contorno é de 0,8.10-5 s, o que

fornece para o modelo um parâmetro β igual a 0,6. Para a malha modelada pelo MEF

o intervalo de tempo é quatro vezes menor.

Os resultados obtidos através do programa implementado neste trabalho são

comparados com aqueles obtidos através de um programa, anteriormente

desenvolvido e exaustivamente testado, de acoplamento acústico-elástico entre o

Método das Diferenças Finitas (MDF), para o meio acústico, e o Método dos

Elementos Finitos (MEF), para o meio elástico. A malha de pontos de diferenças

finitas, neste caso, é estendida na horizontal até uma distância de 1,96 m em relação

ao eixo de axissimetria, de modo que, nos pontos de análise, não haja influência de

reflexões ocorridas no bordo direito do modelo. Para esta comparação, ilustrada pelas

figuras 4.31 a 4.33, são tomados os históricos de pressão nos pontos A e B, e de

deslocamento horizontal no ponto C, referentes à figura 4.30.

0 0.0004 0.0008 0.0012 0.0016 0.002tempo (s)

-6

-4

-2

0

2

4

6

pres

são

(MP

a)

MEC-MEFMDF-MEF

Figura 4.31 – pressão no ponto A – comparação entre métodos

Aplicações 109

0 0.0004 0.0008 0.0012 0.0016 0.002tempo (s)

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

pres

são

(MP

a)

MEC-MEFMDF-MEF

Figura 4.32 – pressão no ponto B – comparação entre métodos

0 0.0004 0.0008 0.0012 0.0016 0.002tempo (s)

-1e-011

-5e-012

0

5e-012

1e-011

1.5e-011

pres

são

(MP

a)

MEC-MEFMDF-MEF

Figura 4.33 – deslocamento horizontal no ponto C – comparação entre métodos

Pela análise dos gráficos, pode-se afirmar que os resultados encontram-se

bastante próximos, especialmente até um instante de tempo entre 0,0016 e 0,0018 s, a

partir do qual começa a haver uma divergência mais pronunciada, principalmente na

pressão no ponto B. Diante das múltiplas reflexões e passagens pelo interior do duto

ocorridas até esse instante, torna-se difícil interpretar este comportamento. O

problema, talvez, possa advir de pequenos erros acumulados na montagem das

matrizes de convolução de elementos de contorno, que a partir de algum instante de

tempo tornam-se relevantes.

Aplicações 110 As figuras 4.34 e 4.35 mostram a comparação entre os resultados de pressão

nos pontos A e B para as configurações com e sem o duto (meio homogêneo

modelado pelo MEC).

0 0.0004 0.0008 0.0012 0.0016 0.002tempo (s)

-6

-4

-2

0

2

4

6

pres

são

(MP

a)

sem dutocom duto

Figura 4.34 – pressão em A – comparação com e sem o duto

0 0.0004 0.0008 0.0012 0.0016 0.002tempo (s)

-1.2

-0.8

-0.4

0

0.4

0.8

1.2

pres

sào

(MP

a)

sem dutocom duto

Figura 4.35 – pressão em B – comparação com e sem o duto

O ponto A, obviamente, não é afetado pela presença do duto até que seja

atingido, em 0,0002 s, por uma reflexão que ocorre na parede interna dele. Percebe-se,

na figura 4.34, que esta influência começa antes que chegue toda a emissão original

da fonte, que dura 0,0003 s.

Aplicações 111 Em relação ao ponto B, percebe-se nitidamente a redução de amplitude em

virtude da colocação do duto, como uma barreira entre fonte e receptor. Nota-se

também que a excitação, com o duto, chega antes, pela aceleração que a onda sofre na

passagem pela parede de aço, visto que este material tem velocidade de propagação da

onda primária cerca de quatro vezes maior que a da água.

A redução da amplitude pode ser medida através da atenuação, dada por (4.5),

onde pct é a pressão na configuração com o tubo e pst a pressão na configuração sem o

tubo.

⎟⎟⎠

⎞⎜⎜⎝

⎛−=

st

ct

pp

A log20 (4.5)

Tomando-se os valores do primeiro pico da figura 4.35, chega-se a uma

atenuação de 14,80. O valor teórico, para teoria de paredes planas, é dado por (4.6).

( )⎟⎟

⎠

⎞

⎜⎜

⎝

⎛⎟⎟⎠

⎞⎜⎜⎝

⎛+= ek

rr

A 22

2

1

2 sen411log10 (4.6)

Em (4.6), r2 e r1 são as partes reais das impedâncias, dadas por (4.7), da parede

e do meio a ela adjacente, respectivamente; e é a espessura da parede e k2 o número de

onda, dado por (4.8), sendo nesta última expressão f a freqüência da onda que se

propaga e c2 a velocidade de propagação da onda primária na parede.

cr ρ= (4.7)

22

2c

fk π= (4.8)

O emprego da expressão (4.6) para uma parede de espessura igual a 0,04 m,

submetida a uma excitação de freqüência igual a 10000 Hz, fornece uma atenuação de

16,05. Pode-se considerar então que o valor obtido através da modelagem numérica

encontra-se satisfatoriamente próximo do valor teórico predito, já que este se refere à

propagação através de paredes planas, que não é o caso da aplicação aqui considerada.

Aplicações 112

4.3.3. Parede circundada por fluido submetida a carga impulsiva

Esta aplicação consiste de uma parede, de seção transversal circular vazada,

constituída de um material hipotético, oca em seu interior e com material fluido em

seu redor. As dimensões do problema são indicadas na figura 4.36.

f(t)MATERIAL 1

MATERIAL 2A B

detalhe da malha do MEF

detalhe da malha do MEC

Figura 4.36 – esquema do modelo para aplicação do item 4.3.3; cotas em metros

As restrições nos contornos inferior e superior do modelo implicam no

emprego de condições de contorno de deslocamentos nulos na malha de elementos

finitos e fluxo nulo na de elementos de contorno.

As propriedades dos materiais 1 e 2 são descritas a seguir:

Material 1: ρ = 2,0 kg/m3; E = 2,66.105 N/m2; ν = 0,30;

Material 2: ρ = 1,0 kg/m3; c = 1436 m/s.

A carga f(t) aplicada na parede equivale a uma carga impulsiva, de valor

unitário, em N/m2, e duração igual a 0,012 s. A análise é estendida até um instante de

tempo igual a 0,1 s.

O meio correspondente ao material 1 é discretizado por elementos finitos

elastodinâmicos quadrados de lado igual a 0,5 m. Já o meio relativo ao material 2 tem

seu contorno discretizado até uma distância de 100 m em relação ao eixo de

axissimetria, por elementos de comprimento igual a 0,5 m (até a 40 m do eixo) e 1,0

m (até o fim).

O passo de tempo de análise, para ambas as malhas, é de 0,3.10-3 s, o que

fornece, na malha do MEC, β igual a 0,86 para os elementos de 0,5 m, e 0,45 para os

Aplicações 113 elementos de 1,0 m. Com o emprego do parâmetro θ igual a 1,1 na análise, a resposta

obtida é estável.

São tomados para análise os deslocamentos na direção horizontal dos pontos

A e B. Comparam-se, nas figuras 4.37 e 4.38, os resultados fornecidos pelo programa

aqui desenvolvido com aqueles obtidos pelo mesmo programa a que se referiu no item

anterior, a fim de se verificar a validade da resposta.

0 0.02 0.04 0.06 0.08 0.1tempo (s)

-0.002

-0.001

0

0.001

0.002

0.003

desl

ocam

ento

(m)

MEC-MEFMDF-MEF

Figura 4.37 – deslocamento horizontal em A – comparação entre métodos

0 0.02 0.04 0.06 0.08 0.1tempo (s)

-0.002

-0.001

0

0.001

0.002

0.003

0.004

desl

ocam

ento

(m)

MEC-MEFMDF-MEF

Figura 4.38 – deslocamento horizontal em B – comparação entre métodos

Aplicações 114 Os aspectos dos gráficos são bem semelhantes. Nota-se que para o ponto A, a

partir de algum instante de tempo, os resultados começam a divergir mais

pronunciadamente. Este comportamento pode ser creditado a diversos fatores:

refinamento insuficiente, deficiência na integração numérica ao longo dos elementos

de contorno, truncamento da malha, entre outros. Apesar disso, o resultado é

considerado bem satisfatório.

Na seqüência, nas figuras 4.39 e 4.40, comparam-se os resultados obtidos

através do modelo da figura 4.36 com aqueles obtidos para a parede caso não

houvesse fluido em seu exterior.

0 0.02 0.04 0.06 0.08 0.1tempo (s)

-0.005

-0.00375

-0.0025

-0.00125

0

0.00125

0.0025

0.00375

0.005

desl

ocam

ento

(m)

com fluidosem fluido

Figura 4.39 – deslocamento horizontal em A – comparação com e sem fluido

0 0.02 0.04 0.06 0.08 0.1tempo (s)

-0.005

-0.00375

-0.0025

-0.00125

0

0.00125

0.0025

0.00375

0.005

desl

ocam

ento

(m)

com fluidosem fluido

Figura 4.40 – deslocamento horizontal em B – comparação com e sem fluido

Aplicações 115

Conforme é de se esperar, as amplitudes de deslocamento para o caso sem a

presença de fluido são bem maiores, já que este fluido representa um obstáculo à

vibração da parede (funciona, a grosso modo, como um amortecedor).

Com isso, o próprio aspecto do gráfico é bem diferente. Sem o fluido, o

“período de vibração” (tempo para que a resposta atinja dois máximos consecutivos) é

maior, como pode ser visto, por exemplo, na figura 4.40, em que, para o caso com o

fluido há praticamente três picos de deslocamento positivo, enquanto que, para o caso

sem fluido, há apenas um. Ou seja, a presença do fluido altera significativamente os

modos de vibração da parede elástica.

Conclusões 116

Conclusões

Neste trabalho, foi desenvolvido um programa computacional para análise de

problemas de propagação de ondas em domínios axissimétricos, envolvendo o

acoplamento entre meios acústicos e elásticos.

Foram brevemente apresentadas as equações básicas que governam o

problema físico, tanto para meios acústicos quanto para elásticos, bem como as

formulações do Método dos Elementos Finitos e do Método dos Elementos de

Contorno, este último para o caso acústico, apenas.

As integrais ao longo do tempo presentes na solução do meio acústico

modelado pelo MEC foram todas desenvolvidas analiticamente, sendo que em alguns

casos foi utilizado o conceito de partes finitas de integrais, simplificando a álgebra

com a qual se lida, em comparação com aquela necessária quando da consideração

das derivadas das funções Heaviside e conseqüente regularização dos integrandos. As

expressões obtidas como resultados da integração tiveram sua validade assegurada em

vista das respostas obtidas nas aplicações efetuadas. Como foi visto, a integração

analítica no tempo fornece respostas mais precisas do que aquelas obtidas quando se

integram as expressões numericamente. No exemplo da barra submetida a

carregamento normal, tomado como referência para comparação, a resposta com as

integrais analíticas apresentou menor amortecimento numérico e praticamente não

apresentou defasagem de período.

Em alguns casos, verificaram-se instabilidades presentes nas respostas em

instantes de tempo de análise considerados prematuros. Esse fenômeno apareceu

principalmente em problemas de domínio fechado, próximo ao eixo de axissimetria,

levando à conclusão de que é proveniente de um mau condicionamento das

expressões obtidas quando o raio é pequeno. O domínio fechado, através das

múltiplas reflexões no contorno, favorece a propagação da instabilidade. Essas

instabilidades foram eliminadas através da adoção do esquema θ, que provou sua

eficácia.

Foi testada a influência do refinamento do passo de tempo de análise em

função da malha de elementos de contorno utilizada. Chegou-se então a uma faixa de

valores mais adequados para o parâmetro β, que relaciona a distância percorrida pela

onda em um passo de tempo e o tamanho do elemento. Para valores muito pequenos,

Conclusões 117 o elemento fica dividido em segmentos também muito pequenos para a integração,

concentrando muito as regiões de gradiente elevado do integrando e então

prejudicando o resultado da integração numérica ao longo do segmento. Já para

valores muito altos, a integração numérica fica prejudicada pois a extensão do

segmento torna-se muito grande em relação ao número de pontos de Gauss adotado.

Pôde-se perceber que o amortecimento numérico da resposta em alguns casos

mostrou-se elevado, para instantes de tempo de análise mais avançados,

comportamento semelhante àquele observado em análises bidimensionais. Essa é uma

característica dos métodos de resolução no domínio do tempo; no entanto, tem um

aspecto mais acentuado no MEC, em comparação com o MEF, por exemplo.

O acoplamento iterativo também foi implementado com sucesso, apresentando

uma precisão considerada satisfatória nos resultados. Através dos exemplos incluídos

neste trabalho, foram tiradas diversas conclusões a respeito de alguns parâmetros

necessários para o acoplamento.

Notou-se que a flexibilidade introduzida pelo acoplamento iterativo quanto à

escolha do passo de tempo adotado em cada malha é de grande utilidade, visto que

cada método apresenta diferentes características para a definição da discretização a ser

adotada.

O parâmetro de relaxação teve seu valor testado, chegando-se à conclusão de

que o valor de referência na literatura era o mais adequado para os problemas com os

quais se lidou.

Foi introduzido um segundo parâmetro de relaxação, para a adoção das

condições de contorno a serem adotadas na primeira malha resolvida (no caso, a de

elementos finitos), ao início de cada iteração, através de uma ponderação entre os

valores obtidos da malha de elementos de contorno e os valores iniciais na iteração

anterior. Sem este parâmetro, houve aplicações em que se verificaram instabilidades.

Como sugestões para trabalhos futuros dando seguimento ao aqui apresentado,

podem-se citar:

Elaboração de um programa que comporte mais de uma interface de

acoplamento. Para isso, cada sistema deve ser resolvido separadamente, sendo que a

cada iteração todas as interfaces passam pelo teste de convergência, sendo adotadas,

na iteração seguinte, as condições de contorno adequadas nos nós de cada interface do

Conclusões 118 problema. O programa atual permite a consideração de apenas dois sistemas, com

uma interface de acoplamento.

Desenvolvimento das expressões advindas das derivadas das funções

Heaviside da solução fundamental, regularizando o integrando das integrais resolvidas

aqui diretamente por partes finitas e obtendo mais uma parcela que considera as

condições iniciais no contorno. Com isso, é possível a consideração de condições

iniciais de pressão nos nós pertencentes ao contorno da malha de elementos de

contorno.

Introdução de aproximação linear para o fluxo ao longo do tempo. Esse

procedimento implica a necessidade de integração analítica no tempo de novas

expressões, necessárias para a montagem da matriz G. Torna-se também possível, em

conjunto com o procedimento descrito no item anterior, a consideração de condições

iniciais de fluxo não-nulas nos nós do contorno

Consideração das integrais de domínio presentes na equação integral do

MEC, a fim de se considerar as condições iniciais para pressão e sua derivada em

relação ao tempo no interior do domínio.

Execução de programação em paralelo, especialmente para o acúmulo das

matrizes geradas no MEC pelo processo de convolução. Como a solução fundamental

axissimétrica tem o aspecto mostrado na figura 2.4, não é possível efetuar o

truncamento até que todos os pontos da malha axissimétrica sejam atingidos pela

excitação emitida pelo ponto mais distante a cada um, dentre todos os anéis

correspondentes aos pontos fontes. Depois desse instante, todas as matrizes anulam-se

automaticamente; entretanto, até aí, não é possível a adoção de qualquer processo de

truncamento, nem mesmo pelo cálculo das matrizes em função da interpolação de

valores obtidos para alguns instantes de tempo, já que a solução fundamental

axissimétrica apresenta um crescimento, com forte variação de gradiente, à medida

que se aproxima o instante em que o ponto mais distante do anel de fontes atinge o

ponto campo em questão.

Introdução do amortecimento na análise acústica por elementos de

contorno. Sugere-se a adoção de uma adaptação do procedimento empregado em JIN

et al. [35], de simples implementação.

Conclusões 119

Obtenção da solução fundamental axissimétrica para modelagem de meios

elásticos pelo MEC. A matemática envolvida na transformação da solução

fundamental tridimensional na axissimétrica é bem mais complexa neste caso; no

entanto, acredita-se que, através de procedimentos semelhantes aos empregados no

caso acústico, pode-se chegar às expressões adequadas, a serem consideradas na

equação integral correspondente.

Bibliografia 120

Bibliografia

[1] POISSON, S.D., Mémoire sur l'équilibre et le mouvement des corps élastiques.

Mém. Acad. Sci. Paris, 1829.

[2] BATHE, K.J., Finite element procedures. 1 ed. New Jersey, Prentice Hall Inc.,

1996.

[3] COOK, R.D., Concepts and applications of finite element analysis. 3 ed. New

York, John Wiley & Sons Inc., 2001.

[4] THOMSON, W.T., Theory of vibration with applications. 1 ed. Englewood

Cliffs, New Jersey, Prentice Hall Inc., 1973.

[5] CLOUGH, R.W., PENZIEN, J., Dynamics of structures. 2 ed. University of

California, Berkeley, McGraw-Hill, 1993.

[6] ZIENKIEWICZ, O.C., TAYLOR, R.L., The finite element method – vol. 1. 5 ed.

Oxford, Butterworth-Heinemann, 2002.

[7] GIVOLI, D., NETA. B. “High-order non-reflecting boundary condition for

dispersive waves”, Wave Motion v. 37, n. 3, pp. 257-271, Mar. 2003.

[8] SARMA, G.S., MALLICK, K., GADHINGLAJKAR, V.R. “Nonreflecting

boundary condition in finite-element formulation for an elastic wave equation”,

Geophysics v. 63, n. 3, pp. 1006-1016, May-Jun. 1998.

[9] DOMINGUEZ, J., Boundary elements in dynamics. Southampton,

Computational Mechanics Publications, 1993.

[10] BECKER, A.A., The boundary element method in engineering. Berkshire,

McGraw-Hill, 1992.

[11] MANSUR, W.J., A time-stepping technique to solve wave propagation problems

using the boundary element method. Ph.D. Thesis, University of Southampton,

Southampton, England, 1983.

[12] SOARES JR., D., Análise dinâmica de sistemas não lineares com acoplamento

do tipo solo-fluido-estrutura por intermédio do método dos elementos finitos e

do método dos elementos de contorno. Tese de D.Sc., COPPE/UFRJ, Rio de

Janeiro, RJ, Brasil, 2004.

[13] VON ESTORFF, O., ANTES, H., “On FEM-BEM coupling for fluid-structure

interaction analysis in the time domain”, International Journal for Numerical

Methods in Engineering v. 31, n. 6, pp. 1151-1168, May 1991.

Bibliografia 121 [14] BELYTSCHKO, T., LU, Y.Y., “A variationally coupled FE-BE method for

transient problems”, International Journal for Numerical Methods in

Engineering v. 37, n. 1, pp. 91-105, Jan. 1994.

[15] YU, G., MANSUR, W.J., CARRER, J.A.M., LIE, S.T., “A more stable scheme

for BEM/FEM coupling applied to two-dimensional elastodynamics”,

Computers and Structures v. 79, n. 8, pp. 811-823, Mar. 2001.

[16] SOARES JR., D., VON ESTORFF, O., “Combination of FEM and BEM by an

iterative coupling procedure”. European Congress on Computational Methods

in Applied Sciences and Engineering, 1099-170, Jyväskylä, Finland, 24-28 July

2004.

[17] SOARES JR, D., VON ESTROFF, O., MANSUR, W.J., “Iterative coupling of

BEM and FEM for nonlinear dynamic analyses”. Computational Mechanics v.

34, n.1, pp. 67-73, Jun. 2004.

[18] SOARES JR, D., VON ESTROFF, O., MANSUR, W.J., “Efficient nonlinear

solid-fluid interaction analysis by an iterative BEM/FEM coupling”.

International Journal for Numerical Methods in Engineering v. 64, n. 11, pp.

1416-1431, Nov. 2005.

[19] GRANNELL, J.J., SHIRRON, J.J., COUCHMAN, L.S., “A hierarchic p-version

boundary–element method for axisymmetric acoustic scattering”, Journal of the

Acoustical Society of America v. 95, n. 5, pp. 2320-2329, May 1994.

[20] TSINOPOULOS, S.V., AGNANTIARIS, J.P., POLYZOS, D., “An advanced

boundary element/fast Fourier transform axisymmetric formulation for acoustic

radiation and wave scattering problems”, Journal of the Acoustical Society of

America v. 105, n. 3, pp. 1517-1526, Mar. 1999.

[21] SOENARKO, B., “A boundary element formulation for radiation of acoustic

waves from axisymmetric bodies with arbitrary boundary conditions”, Journal

of the Acoustic Society of America v. 93, n. 2, pp. 631-639, Feb. 1993.

[22] ISRAIL, A.S.M., BANERJEE, P.K., WANG, H.C., “Time-domain formulations

of BEM for two-dimensional, axisymmetric and three-dimensional scalar wave

propagation”. In: Advanced Dynamic Analysis by Boundary Element Method –

v. 7, chapter 3, London, England, Elsevier Appl. Scie., pp. 75-113, 1992.

[23] BESKOS, D.E., “Boundary element methods in dynamic analysis”, Applied

Mechanics Reviews v. 50, n. 3, pp. 149-197, Mar. 1997.

Bibliografia 122 [24] CZYGAN, O., VON ESTORFF, O., “An analytical fundamental solution of the

transient scalar wave equation for axisymmetric systems”, Computer Methods

in Applied Mechanics and Engineering v. 192, pp. 3657-3671, May 2003.

[25] MORSE, P.M., FESHBACH, H., Methods of theoretical physics. New York,

McGraw-Hill, 1953.

[26] HADAMARD, J., Lecture on Cauchy's problem in linear partial differential

equations. New York, Dover Publications, 1952

[27] KUTT, H.R., Quadrature formulae for finite-part integrals, In: CSIR Special

Report WISK 178, National Research Institute for Mathematical Sciences,

Pretoria, 1975.

[28] ABRAMOWITZ, M., STEGUN, I.A., Handbook of mathematical functions with

formulas, graphs and mathematical tables. 5 ed., New York, Dover

Publications, 1968.

[29] MANSUR, W.J., CARRER, J.A.M., “Two-dimensional transient BEM analysis

for the scalar wave equation: kernels”, Engineering Analysis with Boundary

Elements v. 12, n. 4, pp. 283-288, Oct. 1993.

[30] TELLES, J.C.F., “A self-adaptive coordinate transformation for efficient

numerical evaluation of general boundary element integrals”, International

Journal for Numerical Methods in Engineering v. 24, n. 5, pp. 959-973, May

1987.

[31] BREBBIA, C.A., TELLES, J.C.F., WROBEL, L.C., Boundary element

techniques. Berlin, Springer-Verlag, 1984.

[32] MANSUR, W.J., YU, G., CARRER, J.A.M., LIE, S.T., SIQUEIRA, E.F.N.,

“The θ scheme for time-domain BEM/FEM coupling applied to the 2-D scalar

wave equation”, Communications in Numerical Methods in engineering v. 16,

n. 6, pp. 439-448, Feb. 2000.

[33] ELLEITHY, W.M., AL-GAHTANI, H.J., EL-GEBEILY, M., “Iterative coupling

of FE and BE methods in elastostatics”, Engineering Analysis with Boundary

Elements v. 25, n. 8, 685-695, Sep. 2001.

[34] COHEN, G., JOLY, P., “Fourth order schemes for the heterogeneous acoustic

equation”, Computer Methods in Applied Mechanics and Engineering v.80, pp.

397-407, Jun. 1990.

Bibliografia 123 [35] JIN, F., PEKAU, O.A., ZHANG, C.H., “A 2-D time domain boundary element

method with damping”, International Journal for Numerical Methods in

Engineering v. 51, n. 6, pp. 647-661, Jun. 2001.

[36] NOWACKI, W., Dynamic of elastic systems. New York, John Wiley & Sons.

Inc., 1963.

[37] MORSE, P.M., Vibration and sound. 2 ed. New York, McGraw-Hill, 1948.

Expressões analíticas para as integrais no tempo 124

Anexo A

Expressões analíticas para as integrais no tempo

Neste anexo, são apresentadas as expressões analíticas para as integrais ao

longo do tempo presentes na montagem das matrizes H e G. O procedimento para

obtenção destas expressões é mostrado no item 2.2.

A.1. Expressões para g

− Caso 2:

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −−−

=d

rdrdru

udF

dg i

i

22

22

22

2 ,2

1π

(A.1)

− Caso 3:

⎟⎟⎠

⎞⎜⎜⎝

⎛ −=

drdK

dg

22

221

π (A.2)

− Caso 4:

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −−

−−

−⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −−−

=

drd

rdru

udF

d

drd

rdru

udF

dg

f

f

i

i

22

22

22

2

22

22

22

2

,2

1

,2

1

π

π (A.3)

− Caso 5:

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −−

−=

drd

rdud

Fd

g f22

22

22

2 ,2

1π

(A.4)


A.2. Expressões para hI e hF

As expressões finais para hI e hF são obtidas pelas fórmulas (A.5) e (A.6).

( )432212122112 ECEuCECEuC

tch ffi −++−

Δ=

π (A.5)

( )432212122112 ECEuCECEuC

tch iif +−−

Δ=

π (A.6)

As expressões E1, E2, E3 e E4 são fornecidas a seguir. Estas expressões

relacionam-se com a separação do integrando nas Parcelas 1, 2, 3 e 4 realizada no

item 2.2.1, representando o resultado da integral ao se desconsiderar as constantes do

integrando. As demais variáveis têm o mesmo significado que aquele indicado no

Capítulo 2.

E1:

− Caso 2:

( ) ( )( ) ( )( )

( )( )

( ) ⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −−−

−

+−

−⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −−−

−+

+−−−

−−−=

drd

rdru

udE

drrdrd

drd

rdru

udF

drd

ruudurd

udruE

i

i

i

i

iii

ii

22

22

22

2222

22

22

22

22

222

2222222

2222

,

,2

1

(A.7)

− Caso 3:

( )( )

( ) ⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

−−

−

+−⎟

⎟⎠

⎞⎜⎜⎝

⎛ −

−= 22

22

2222

2222

222

21rdru

udE

drrdrd

drdK

drdE i

i

(A.8)


− Caso 4:

( ) ( )( ) ( )( )

( )( )

( )( ) ( )

( ) ( )( )

( )( )

( ) ⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −−

−

−

++

+⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −−

−

−−

−−−−

−−−−

−⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −−−

−

+−

−⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −−−

−+

+−−−

−−−=

drd

rdru

udE

drrdrd

drd

rdru

udF

drd

ruudurd

udru

drd

rdru

udE

drrdrd

drd

rdru

udF

drd

ruudurd

udruE

f

f

f

f

fff

ff

i

i

i

i

iii

ii

22

22

22

2222

22

22

22

22

222

2222222

2222

22

22

22

2222

22

22

22

22

222

2222222

2222

,

,2

,

,2

1

(A.9)

− Caso 5:

( ) ( )[ ]( ) ( )( )

( )( )

( ) ⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −−−

−

+−

−⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −−−

−+

+−−−

−−−−=

drd

rdud

Edrrd

rd

drd

rdud

Fdrd

ruudrdrd

uddruruE

i

i

ff

fff

22

22

22

2222

22

22

22

22

222

222222222

222222

,

,2

1

(A.10)


E2:

− Caso 2:

( ) ( )( ) ( )( )2222222

2222

2ruudrd

udruE

ii

ii

−−−

−−−= (A.11)

− Caso 3:

02 =E (A.12)

− Caso 4:

( ) ( )( ) ( )( )

( ) ( )( ) ( )( )2222222

2222

2222222

2222

2ruudrd

udru

ruudrd

udruE

ff

ff

ii

ii

−−−

−−−−

−−−

−−−=

(A.13)

− Caso 5:

( ) ( )[ ]( ) ( )( )2222222

2222

2ruudrd

udruE

ff

ff

−−−

−−−−= (A.14)


E3:

− Caso 2:

( ) ( )

⎥⎥⎦

⎤

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −−−

−

⎢⎢⎣

⎡−

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −−−

−+

−−

−=

drd

rdru

udE

drd

rdru

udF

rddudrdu

ruE

i

i

i

iii

i

22

22

22

22

22

22

222222

22

,

,13

(A.15)

− Caso 3:

( ) ⎥⎥⎦

⎤

⎢⎢⎣

⎡

⎟⎟⎠

⎞⎜⎜⎝

⎛ −−⎟

⎟⎠

⎞⎜⎜⎝

⎛ −−

=d

rdEd

rdKrdd

E2222

22

13 (A.16)

− Caso 4:

( ) ( )

( )

( ) ⎥⎥⎦

⎤

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −−−

⎢⎢

⎣

⎡−⎟

⎟

⎠

⎞

⎜⎜

⎝

⎛ −−

−

−−

−−−

−−

⎥⎥⎦

⎤

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −−−

−

⎢⎢⎣

⎡−

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −−−

−+

−−

−=

drd

rdru

udE

drd

rdru

udF

rdd

udrdu

ru

drd

rdru

udE

drd

rdru

udF

rddudrdu

ruE

i

i

f

f

ff

fi

i

i

iii

i

22

22

2222

22

22

22

2222

2222

22

22

22

22

22

222222

22

,,1

,

,13

(A.17)

− Caso 5:

( ) ( )

⎥⎥

⎦

⎤

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ −−

−−

⎢⎢

⎣

⎡−⎟

⎟

⎠

⎞

⎜⎜

⎝

⎛ −−

−⋅

−+

−−

−−=

drd

rdud

E

drd

rdud

Frddudrdd

ruuE

f

f

i

ff

22

22

22

22

22

22

2222222

22

,

,13

(A.18)


E4:

− Caso 2:

( ) 2222

22

4i

i

udrd

ruE

−−

−= (A.19)

− Caso 3:

04 =E (A.20)

− Caso 4:

( ) ( ) 2222

22

2222

22

4f

f

i

i

udrd

ru

udrd

ruE

−−

−−

−−

−= (A.21)

− Caso 5:

( ) 2222

22

4f

f

udrd

ruE

−−

−−= (A.22)

Integrais elípticas 130

Anexo B

Integrais elípticas

Para o desenvolvimento analítico das integrais presentes nas expressões (2.12)

a (2.14) é necessária a compreensão do conceito de integral elíptica. Neste anexo, são

apresentadas algumas definições básicas, bem como o procedimento numérico

empregado para se calcular estas integrais. Mais informações sobre este assunto

podem ser obtidas em ABRAMOWITZ e STEGUN [28].

Uma integral é classificada como elíptica se o seu integrando R(x,y) é uma

função racional de x e y, sendo y2 um polinômio de terceira ou quarta ordem em x. De

modo geral, as integrais elípticas não podem ser expressas em termos de funções

elementares, sendo representadas através de séries ou aproximadas numericamente.

B.1. Tipos de integrais elípticas

Toda integral elíptica, através de determinadas transformações, pode ser

expressa através de três formas canônicas. Neste trabalho, são tratados apenas os dois

primeiros tipos de forma canônica, já que o terceiro não tem utilidade para as integrais

com as quais se lida. O primeiro tipo é representado pela expressão (B.1), e o segundo

por (B.2).

( ) ( ) ( ) ( )( )[ ]∫∫ −−−−=−==

x

dttmtdmxFF0

2/1222

0

2/122 11sensen1,\ϕ

θθααϕ

(B.1)

( ) ( ) ( ) ( ) ( )∫∫ −−=−==−

x

dttmtdmxEE0

2/1222/12

0

2/122 11sensen1,\ϕ

θθααϕ

(B.2)

Em (B.1) e (B.2) são adotadas as definições expressas em (B.3) e (B.4).

ϕsenx = (B.3)

αsenm = (B.4)

Integrais elípticas 131

Nas expressões contidas no Anexo A, adota-se sempre a representação das

integrais elípticas correspondente às formas F(x,m) e E(x,m) presentes em (B.1) e

(B.2). Nas expressões (B.1) a (B.4), m é chamado de parâmetro da integral, α de

ângulo modular e ϕ de amplitude.

Quando ϕ = π/2, ou seja, x = 1, as integrais são ditas completas, sendo

representadas por K(m), para o primeiro tipo, e E(m), para o segundo.

B.2. Cálculo das integrais

As integrais elípticas apresentadas em B.1 são, neste trabalho, calculadas

numericamente. Para as integrais completas, adota-se a expansão através de

polinômios, que são combinados conforme as expressões (B.5) e (B.6).

1ln)( mBAmK KK −= (B.5)

1ln)( mBAmE EE −= (B.6)

Para (B.5) e (B.6) adotam-se as definições expressas em (B.7) a (B.11).

414

313

212110 mamamamaaAK ++++= (B.7)

414

313

212110 mbmbmbmbbBK ++++= (B.8)

414

313

212110 mcmcmcmccAE ++++= (B.9)

414

313

212110 mdmdmdmddBE ++++= (B.10)

21 1 mm −= (B.11)

Os coeficientes das expressões (B.7) a (B.11) são definidos na Tabela A.1.

i ai bi ci di 0 1,38629436112 0,5 1 0 1 0,09666344259 0,12498593597 0,44325141463 0,24998368310 2 0,03590092383 0,06880248576 0,06260601220 0,09200180037 3 0,03742563713 0,03328355346 0,04757383546 0,04069697526 4 0,01451196212 0,00441787012 0,01736506451 0,00526449639

Tabela A.1 – coeficientes para polinômios de aproximação de integrais elípticas completas

Integrais elípticas 132 Para as integrais elípticas de primeiro e segundo tipos incompletas, se o

parâmetro m é igual a 1, utilizam-se as expressões (B.12) e (B.13).

( ) ⎟⎟⎠

⎞⎜⎜⎝

⎛

−

+=

211ln,

xxmxF (B.12)

( ) xmxE =, (B.13)

Para o cálculo das integrais incompletas, emprega-se o chamado “Processo

Aritmético-Geométrico”.

Neste processo, definem-se inicialmente quatro parâmetros:

0;;1;0,1 01

02

00 ===−== − Gxsenmba ϕϕ (B.14)

Procedem-se então N passos, sendo que em cada passo i são calculados os

parâmetros ai, bi, ci, ϕi e Gi, expressos em (B.15) e (B.16).

( ) ( )2

;;2

1111

11 −−−−

−− −==

+= ii

iiiiii

iba

cbabba

a (B.15)

iiiiii

iii cGGtg

ba

tg ϕϕϕϕ sen; 11

1 +=⎟⎟⎠

⎞⎜⎜⎝

⎛+= −

−− (B.16)

O processo pára no N-ésimo passo, de modo que cN = 0, ou seja, aN-1 = bN-1.

Então, obtêm-se F(x,m) e E(x,m) pelas expressões (B.17) e (B.18).

( )N

NN

amxF

2,

ϕ= (B.17)

( ) ( )NG

mKmEmxFmxE +=

)()(,, (B.18)

Parte finita de integrais 133

Anexo C

Parte finita de integrais

Dependendo dos limites de integração presentes nas expressões (2.18) e

(2.19), que fornecem respectivamente hI e hF, o integrando apresenta uma

singularidade não-integrável, visto que o denominador possui fatores elevados a 3/2.

Neste caso, o cálculo desta integral passa a ser possível apenas no sentido de parte

finita.

Neste anexo, é apresentado o conceito de parte finita de integral, mostrando o

modo de se calcular a parte finita das integrais com as quais se trabalha em problemas

axissimétricos em Elementos de Contorno.

A integral (C.1), apesar de conter um integrando que apresenta singularidade

em b (limite infinito), possui resultado finito, como se pode inferir a partir do

resultado de sua integral indefinida. Portanto, diz-se que ela apresenta uma

singularidade integrável.

abxb

dxb

a

−=−∫ 2 (C.1)

Ao se derivar o primeiro membro de (C.1) em relação a b, no entanto, chega-

se à expressão (C.2), que apresenta duas parcelas, a primeira sem significado pelo fato

de o integrando possuir um termo singular em b elevado a 3/2 no denominador, e a

segunda por não possuir limite finito em b. No entanto, apesar de nenhuma das

parcelas possuir limite finito em x = b, sabe-se que a soma delas possui, visto que

produz o mesmo resultado que a derivada de uma função sabidamente derivável

(segundo membro de (C.1)).

( ) bx

b

axbxb

dx

=⎥⎦

⎤⎢⎣

⎡−

+−

− ∫ 121

2/3 (C.2)

Do mesmo modo, a integral (C.3) não possui limite finito quando x tende a b,

visto que o resultado da integral indefinida a ela correspondente também não

Parte finita de integrais 134 apresenta este limite. Esta integral possui, então, uma singularidade não-integrável em

x = b, ou seja, tanto o integrando quanto a própria integral não possuem limites finitos

quando x tende a b.

( )∫ −

x

a

dxxbxA

2/3

)( (C.3)

No entanto, pode-se provar [26] que a expressão (C.4) possui limite finito

nesta situação, em x = b, contanto que B(x) satisfaça à condição de Lipschitz

( )1212 )()( xxKxBxB −<− e B(b) = - 2A(b).

( ) xbxBdx

xbxA

x

a−

+−∫ )()(

2/3 (C.4)

O limite de (C.4) quando x tende a b é chamado de parte finita da integral

(C.3), sendo indicado da forma apresentada em (C.5).

( )∫ −

b

a

dxxbxA

2/3

)( (C.5)

Para se efetuar a integral apresentada por (C.5), primeiro deve ser definida a

parte finita para A(x) =1, conforme exposto em (C.6), onde se faz:

B(x) = - 2A(x).

Na expressão (C.6), aproveita-se o fato de que (C.2) deve ser igual à derivada

do segundo membro de (C.1).

( ) ( )

( )( )

ab

abbxbxb

dx

xbdx

xbdx

xb

bx

b

a

x

abx

b

a

−−

=

=−∂∂

−=⎟⎟

⎠

⎞

⎜⎜

⎝

⎛⎥⎦

⎤⎢⎣

⎡−

+−

−−=

=⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

−−

−=

−

=

→

∫

∫∫

2

221212

21lim1

2/3

2/32/3

(C.6)

=

=

Parte finita de integrais 135 Obtida a expressão (C.6), substitui-se em (C.5) A(x) por [A(x) – A(b)] + A(b),

chegando-se finalmente à expressão (C.7).

( ) ( ) ( )

( ) abbAdx

xbbAxA

dxxbbAdx

xbbAxAdx

xbxA

b

a

b

a

b

a

b

a

−−

−−

=

=−

+−

−=

−

∫

∫∫∫)(2)()(

)()()()(

2/3

2/32/32/3

(C.7)

A primeira parcela do último membro em (C.7) é uma integral imprópria, que

apresenta um resultado finito. Deste modo, fica completamente definido o cálculo da

parte finita de uma integral cujo integrando contém uma singularidade elevada a 3/2

no denominador.

= =

Expressões para integrais das singularidades na frente de onda 136

Anexo D

Expressões para integrais das singularidades na frente de onda

No item 2.3.2 foram enunciados todos os casos onde pode haver singularidade,

nas expressões provenientes das integrais ao longo do tempo, para a frente de onda.

Foram apresentados também os resultados do processo de integração semi-analítico

no espaço para dois deles. Neste anexo, são apresentados os resultados para os demais

casos, através das expressões A(x) e A(ui) ou A(uf), conforme a singularidade, para

singularidade em r, e B(x) e B(ui) ou B(uf), conforme a singularidade, para

singularidade em d, equivalentes àqueles resultados obtidos em 2.3.2.

D.1. Singularidade em r

− Caso 2 (r = ui):

Da expressão (A.7):

( )( )

)()()(4

)()()()( '11

2/31

11222

22

iij

ii

ij

i

i uCuuxu

uAxCxxurd

udxA η

ξη

−

−=⇒−

−−=

(D.1)


( )( ) )()()(4)()()()( '

112/3

111222

22

iij

iiji uCuuxuAxCxx

rd

udxA ηξη −−=⇒

−

−−=

(D.2)

− Caso 4 (r = uf):


( )( )

)()()(4

)()()()( '11

2/31

11222

22

ffj

ff

fj

f

i uCuuxu

uAxCxxurd

udxA η

ξη

−

−=⇒−

−=

(D.3)



( )( ) )()()(4)()()()( '

112/3

111222

22

ifj

fijf uCuuxuAxCxx

rd

udxA ηξη −−=⇒

−

−=

(D.4)

− Caso 5 (r = uf):


( )( )

)()()(4

)()()()( '11

2/31

1122222

222

iij

ff

fjff uCuux

uuAxCxx

rdrd

udduxA η

ξη

−

=⇒−

−=

(D.5)


( )( ) )()()(4)()()()( '

112/3

111222

22

iij

ffjf uCuuxuAxCxx

rd

udxA ηξη −=⇒

−

−=

(D.6)

D.2. Singularidade em d

− Caso 2 (d = ui):


( )( )

)()()(4

)()()()( '11

2/31

11222

22

iij

ii

ij

i

i uCuuxu

uBxCxxurd

ruxB η

ξη

−

=⇒−

−=

(D.7)



( )( ) )()()(4)()()()( '

112/3

111222

22

iij

iiji uCuuxuBxCxx

rd

ruxB ηξη −=⇒

−

−=

(D.8)


( )( )

)()()(4)()()()( '11

2/11

1122

22

iij

ii

ij

i

i uCuuxu

uBxCxxurd

ruxB η

ξη

−

=⇒−

−=

(D.9)


( ) ( ) )()()(4)()()()( '11

2/111122

22

iij

iiji uCuuxuBxCxx

rdru

xB ηξη −=⇒−

−=

(D.10)

− Caso 4 (d = ui):

Da expressão (A.9): igual a (D.7)





− Caso 5 (d = uf):


( )( )

)()()(4)()()()( 21

2/31

21222

22

ffj

ff

fj

f

f uCuuxu

uBxCxxurd

ruxB η

ξη

−

−=⇒−

−−=

(D.11)


( )( ) )()()(4)()()()( 21

2/3121222

22

ffj

ffjf uCuuxuBxCxx

rd

ruxB ηξη −−=⇒

−

−−=

(D.12)


( )( )

)()()(4)()()()( 21

2/11

21222

22

ffj

ff

fjff uCuux

uuBxCxx

rdd

ruuxB η

ξη

−

−=⇒−

−−=

(D.13)


( ) ( ) )()()(4)()()()( 212/1

12122

22

ffj

ffjf uCuuxuBxCxx

rd

ruxB ηξη −−=⇒

−

−−=

(D.14)

Expressões analíticas para os termos relativos a fontes 140

Anexo E

Expressões analíticas para os termos relativos a fontes

Neste anexo, são apresentadas as expressões analíticas para as integrais ao

longo do tempo presentes na montagem do vetor S. O procedimento para obtenção

destas expressões é mostrado no item 2.4.

As expressões finais para wI e wF são obtidas pelas fórmulas (E.1) e (E.2).

( )212

12 WWu

tcw fI +−

Δ=

π (E.1)

( )212

12 WWu

tcw iF −

Δ=

π (E.2)

As expressões W1 são idênticas às expressões (A.1) a (A.4), multiplicadas por

2π2; já as expressões W2 são fornecidas a seguir, de acordo com o caso.

W2:

− Caso 2:

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

−−

−= 22

22

2rdru

arcsenW i (E.3)

− Caso 3:

22 π

−=W (E.4)

− Caso 4:

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

−

−+

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

−−

−= 22

22

22

22

arcsenarcsen2rdru

rdru

W fi (E.5)

Expressões analíticas para os termos relativos a fontes 141

− Caso 5:

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

−

−+−= 22

22

22

rdru

arcsenW fπ (E.6)

Soluções analíticas de interesse 142

Anexo F

Soluções analíticas de interesse

Neste anexo são apresentadas as soluções analíticas para os problemas de que

tratam as aplicações apresentadas nos itens 4.1.1 a 4.1.4.

F.1. Barra engastada em uma extremidade e livre na outra, submetida a

carregamento axial de impacto

Esta solução é fornecida em NOWACKI [36], sendo válida para barra

homogênea. Originalmente deduzida para uma barra elástica, submetida a um

carregamento concentrado P(t) em sua extremidade, pode ser equiparada ao caso de

uma barra acústica, na qual, no lugar de carga, aplica-se fluxo, e no lugar de

deslocamentos, trata-se de pressões. Na aplicação do item 4.1.1, aplicou-se um fluxo

por unidade de área uniformemente distribuído, que, multiplicado pela área da barra,

fornece a carga concentrada.

P(t)

Figura F.1 – barra engastada submetida a carga normal

Carregamento: P(t) = P0.H(t – t0)

Deslocamento axial (ou pressão, no caso acústico):

( )( )∑

∞

=

−

⎭⎬⎫

⎩⎨⎧

⎟⎟⎠

⎞⎜⎜⎝

⎛⎟⎠⎞

⎜⎝⎛ −

−⎟⎠⎞

⎜⎝⎛ −

−−

=1

2

10

2 212cos1

212sen

1218),(

n

n

Lctn

Lxn

nEALP

txu πππ

(F.1)

Em (F.1):

A – área da seção transversal;

E – módulo de elasticidade, ou compressibilidade, no caso acústico;

c – velocidade de propagação da onda primária no meio.


F.2. Cavidade esférica em meio infinito

a

Figura F.2 – cavidade esférica em meio infinito

Para uma cavidade esférica em meio infinito (figura F.2), submetida a um

fluxo uniformemente distribuído em sua superfície, descrito ao longo do tempo por

uma função Heaviside (q(t) = q0.H(t – t0)), a pressão em qualquer ponto do meio, a

uma distância r do centro da cavidade, em qualquer instante de tempo t, é dada por

(F.2), conforme deduzido em CZYGAN e VON ESTORFF [24].

( )

( )[ ])1),(2

0 arctHeraq

trp aarct

−−⎟⎟⎠

⎞⎜⎜⎝

⎛−=

−−−

(F.2)

F.3. Membrana circular submetida a um campo de velocidades iniciais

A solução contida na equação (F.3) é deduzida em MORSE [37], para o

deslocamento vertical η em um ponto de raio r em relação ao centro e ângulo φ em

relação ao eixo horizontal em planta (figura F.3), contido em uma membrana circular

de raio a, cujo material possui velocidade de propagação da onda primária c,

submetida a um campo qualquer de velocidades iniciais v0(r,ϕ). A função J representa

função de Bessel, com seu respectivo argumento.


Figura F.3 – membrana circular

( ) ( ) ( )∑ ∑∑∞

=

∞

=

∞

= ⎪⎭

⎪⎬⎫

⎪⎩

⎪⎨⎧

+=1 0

000

2sen2sen,,n m

mnmnmnm

mnemnemn tCtCtr πνψπνψϕη (F.3)

Em (F.3):

( )∫ ∫Λ=

a

emnmnmn

emn drdrva

C0

2

0

022 ,2

1π

ϕψϕνπ

(F.4)

( )[ ]

( )[ ]2010

212

1

nn

mnmmn

J

J

πβ

πβ

=Λ

=Λ − (F.5)

mnmn ac βν

2= (F.6)

41

2−+≅

mnmnβ (F.7)

( )

( ) ( )

( ) ( ) ⎟⎠⎞

⎜⎝⎛=

⎟⎠⎞

⎜⎝⎛=

⎟⎠⎞

⎜⎝⎛=

ar

Jmr

ar

Jmr

ar

Jr

mnmmn

mnmemn

nne

πβϕϕψ

πβϕϕψ

πβϕψ

sen,

cos,

,

0

000

(F.8)


F.4. Barra acústica fechada em uma extremidade e livre na outra,

submetida a pressão de impacto

A solução apresentada em (F.9) refere-se à aplicação apresentada no item

4.1.4, que consiste de uma barra com fluxo nulo em uma extremidade e pressão

aplicada na outra. Esta pressão se comporta ao longo do tempo como uma função

Heaviside (p(t) = p0.H(t – t0)). O esquema é semelhante àquele dado pela figura F.1.

( )∑∞

=

−

⎭⎬⎫

⎩⎨⎧

⎟⎟⎠

⎞⎜⎜⎝

⎛⎟⎠⎞

⎜⎝⎛ −

−⎟⎠⎞

⎜⎝⎛ −

−−

=1

1

212cos1

212sen

1214),(

n

n

Lctn

Lxn

ntxp ππ

π (F.9)

Documents

ACOPLAMENTO MEC-MEF PARA PROBLEMAS ACÚSTICO