ads:
ads:
Livros Grátis
http://www.livrosgratis.com.br
Milhares de livros grátis para download.
ads:
ii
Aos meus pais
Waldomiro e Guaraciaba
Aos meus filhos
Nícolas, Yuri e Mikhail.
ix
Pai e Mãe
"De sua existência resta comigo o exemplo, a saudade imensa, eterno agradecimento,
além do pesar por não poder abraçá-los agora e partilhar juntos a alegria da tarefa
cumprida”.
xi
AGRADECIMENTOS
A Deus, por me ter dado serenidade e paciência nos momentos difíceis da minha vida.
Ao Prof. Dr. Isaías Vizotto, pela orientação, pela amizade e pela disponibilidade sempre
dispensada.
Ao Prof. Dr. Leandro Palermo Júnior pelos esclarecimentos na fase inicial do trabalho.
Ao Prof. Dr. Wilson Sérgio Venturini pelos ensinamentos e pela colaboração em
diversas partes do trabalho.
Ao Prof. Dr. Humberto Breves Coda pela atenção e pelo auxílio prestado.
Ao Prof. Dr. Renato Soliani pelo apoio e incentivo.
A todos os professores e funcionários, que de alguma forma, contribuíram para o
desenvolvimento desse trabalho.
A minha família pelo apoio e estímulo, sem os quais este trabalho não teria sido
possível.
xiii
______________________________________________________________________
ÍNDICE
CAPÍTULO 1 INTRODUÇÃO ...................................................................................... 4
1.1 C
ONSIDERAÇÕES
G
ERAIS
......................................................................................5
1.2 O
RGANIZAÇÃO DO
T
RABALHO
............................................................................10
CAPÍTULO 2 ELASTOSTÁTICA............................................................................. 12
2.1 I
NTRODUÇÃO
.......................................................................................................13
2.2 E
STADO
P
LANO DE
T
ENSÃO
................................................................................13
2.3 E
STADO
P
LANO DE
D
EFORMAÇÃO
.......................................................................14
2.4 E
QUAÇÕES DE
E
QUILÍBRIO
..................................................................................16
2.5 R
ELAÇÃO
D
EFORMAÇÃO
- D
ESLOCAMENTO
.......................................................21
2.6 R
ELAÇÃO
T
ENSÃO
- D
EFORMAÇÃO
.....................................................................22
2.7 E
QUAÇÕES DE
N
AVIER
........................................................................................23
2.8 C
ONDIÇÕES DE
C
ONTORNO
.................................................................................24
2.9 E
QUAÇÃO
I
NTEGRAL PARA
P
ONTOS DE
D
OMÍNIO
................................................26
2.10 S
OLUÇÃO
F
UNDAMENTAL
...................................................................................31
2.11 E
QUAÇÃO
I
NTEGRAL PARA
P
ONTOS DO
C
ONTORNO
............................................34
2.12 R
EGIÕES
I
NFINITAS
..............................................................................................40
2.13 T
ENSÕES NOS
P
ONTOS
I
NTERNOS
........................................................................43
2.14 T
ENSÕES NOS
P
ONTOS DO
C
ONTORNO
.................................................................44
CAPÍTULO 3 ELASTODINÂMICA ......................................................................... 46
3.1 I
NTRODUÇÃO
.......................................................................................................47
3.2 E
QUAÇÕES
G
ERAIS DA
E
LASTODINÂMICA
...........................................................47
3.3 I
DENTIDADE DE
S
OMIGLIANA
..............................................................................51
3.4 E
QUAÇÃO
I
NTEGRAL DE
C
ONTORNO
...................................................................52
3.5 E
QUAÇÃO
I
NTEGRAL PARA
P
ONTOS
E
XTERNOS
..................................................53
3.6 E
QUAÇÃO
I
NTEGRAL DAS
C
OMPONENTES DE
T
ENSÃO
........................................53
3.7 V
IBRAÇÕES
L
IVRES DE
C
ORPOS
F
INITOS
.............................................................54
2
CAPÍTULO 4 MÉTODO DOS ELEMENTOS DE CONTORNO (MEC) ............. 56
4.1 F
ORMULAÇÃO
E
STÁTICA
.....................................................................................57
4.1.1 Introdução ..................................................................................................... 57
4.1.2 Discretização Geométrica e das Variáveis ................................................... 59
4.1.2.1 Elemento Linear ..................................................................................... 59
4.1.2.2 Elemento Quadrático .............................................................................. 65
4.1.3 Equações Integrais na Forma Matricial ....................................................... 68
4.1.4 Matriz dos Coeficientes de Influência do Elemento ...................................... 72
4.1.4.1 Submatriz g............................................................................................ 75
4.1.4.2 Submatriz h........................................................................................... 77
4.1.5 Condensação de Deslocamentos em Nós Duplos.......................................... 81
4.2 F
ORMULAÇÃO
D
INÂMICA
....................................................................................84
4.2.1 Introdução ..................................................................................................... 84
4.2.2 Formulação com Reciprocidade Dual .......................................................... 86
4.2.3 Formulação com Integração Direta.............................................................. 95
4.2.4 Obtenção da Solução Numérica.................................................................. 109
CAPÍTULO 5 MÉTODO DOS ELEMENTOS FINITOS (MEF) ......................... 115
5.1 F
ORMULAÇÃO
E
STÁTICA
...................................................................................116
5.2 F
ORMULAÇÃO
D
INÂMICA
..................................................................................121
5.3 E
LEMENTO
F
INITO
U
NIDIMENSIONAL
................................................................122
CAPÍTULO 6 MÉTODOS DE INTEGRAÇÃO NUMÉRICA .............................. 128
6.1 I
NTRODUÇÃO
.....................................................................................................129
6.2 M
ÉTODO DE
N
EWMARK
.....................................................................................129
6.3 M
ÉTODO DE
H
OUBOLT
......................................................................................133
6.4 A
PLICAÇÃO DOS
I
NTEGRADORES AO
M
ÉTODO DOS
E
LEMENTOS DE
C
ONTORNO
............................................................................................................................... 136
6.4.1 Método de Newmark Aplicado ao MEC...................................................... 137
6.4.2 Método de Houbolt Aplicado ao MEC ........................................................ 138
CAPÍTULO 7 COMBINAÇÃO DOS MÉTODOS ELEMENTOS FINITOS -
ELEMENTOS DE CONTORNO............................................................................................. 140
3
7.1 I
NTRODUÇÃO
.....................................................................................................141
7.2 E
QUAÇÕES
B
ÁSICAS DA
F
ORMULAÇÃO
MEC
E
MEF........................................142
7.3 T
ÉCNICA DE
A
COPLAMENTO
MEC – MEF........................................................144
CAPÍTULO 8 APLICAÇÕES NUMÉRICAS ......................................................... 148
8.1 A
NÁLISE
E
STÁTICA
...........................................................................................149
8.1.1 Exemplo 1: Viga Solicitada por Força Normal Uniformemente Distribuída
Constante na Extremidade Livre ......................................................................................... 149
8.1.2 Exemplo 2: Viga Solicitada por Força Cortante na Extremidade Livre..... 153
8.1.3 EXEMPLO 3: Viga Solicitada por Momento Concentrado na Extremidade
Livre..................................................................................................................................... 161
8.2 C
ONDENSAÇÃO
E
STÁTICA
.................................................................................167
8.3 MEC – V
IBRAÇÃO
L
IVRE DE UMA
V
IGA EM
B
ALANÇO
.....................................171
8.4 MEC – V
IBRAÇÃO
L
IVRE
: A
RCO
......................................................................176
8.5 A
NÁLISE
T
RANSIENTE DE
P
ROPAGAÇÃO DE
O
NDAS
E
LÁSTICAS
: V
IGA
E
NGASTADA
S
UBMETIDA A
C
ARREGAMENTO
T
RANSVERSAL
S
ÚBITO EM SUA
E
XTREMIDADE
. 179
8.6 MEC - A
NÁLISE
T
RANSIENTE DE
P
ROPAGAÇÃO DE
O
NDAS
E
LÁSTICAS
: V
IGA
E
NGASTADA
S
UBMETIDA A
C
ARREGAMENTO
L
ONGITUDINAL
S
ÚBITO EM SUA
E
XTREMIDADE
182
8.7 A
COPLAMENTO
MEF-MEC–A
NÁLISE
T
RANSIENTE
: C
ARREGAMENTO
T
RANSVERSAL
..........................................................................................................................187
CAPÍTULO 9 CONCLUSÕES ................................................................................. 190
REFERÊNCIAS BIBLIOGRÁFICAS ........................................................................ 194
ANEXO 1: NOTAÇÃO INDICIAL CARTESIANA.................................................. 211
ANEXO 2: CÁLCULO DO JACOBIANO DA TRANSFORMAÇÃO PARA O
ELEMENTO LINEAR ............................................................................................................. 214
ANEXO 3: INTEGRAÇÃO NUMÉRICA DE GAUSS.............................................. 218
ANEXO 4: OBTENÇÃO DA SOLUÇÃO DESLOCAMENTO ............................... 222
ANEXO 5: ASPECTOS DA IMPLEMENTAÇÃO DO SOFTWARE......................227
4
CAPÍTULO 1 INTRODUÇÃO
5
1.1 CONSIDERAÇÕES GERAIS
A crescente exigência da engenharia moderna, aliada à necessidade de se obter
resultados úteis para os problemas em estudo, faz com que a maior parte das
pesquisas na área de mecânica dos sólidos esteja voltada para os métodos numéricos.
Com a utilização de computadores cada vez mais potentes, tornou-se possível a
análise detalhada de problemas mais gerais que envolvem sistemas com milhares de
equações.
O método mais difundido e ainda o mais utilizado é o Método dos Elementos
Finitos (MEF), pela facilidade em modelar casos com geometria complexa, problemas
não lineares e com variação das propriedades no domínio.
Apesar da eficiência, pela versatilidade e pela possibilidade de modelar a maior
parte dos problemas da engenharia de maneira real, o MEF exige a exaustiva
discretização de todo o domínio para a resolução de um problema. Surge, entre outras,
a dificuldade de se modelar domínio infinito, como é necessário em problemas da
elastodinâmica, onde existe a propagação de ondas para o infinito. Como a malha é
finita aparecem implicações relacionadas com a reflexão das ondas provocadas pela
formulação de um contorno fictício.
No final da década de 70, destacou-se uma técnica conhecida como Método dos
Elementos de Contorno (MEC), ganhando notoriedade graças a uma série de
características vantajosas. Uma delas é diminuição considerável nos dados de entrada,
bem como do esforço computacional, uma vez que o sistema gerado tem dimensões
menores que o obtido nos métodos que discretizam o domínio. Outra vantagem é a
habilidade em tratar domínios complicados, como por exemplo, os meios infinitos, pois
neste método apenas a discretização do contorno do domínio é necessária. A idéia
básica do método consiste em relacionar as variáveis dos diferentes pontos do contorno
através de funções analíticas, resultando numa série de coeficientes de influência, os
6
quais são colocados em forma matricial. As condições de contorno são aplicadas de
maneira semelhante ao MEF. Uma característica do sistema de equações gerado é ter
a matriz final cheia e não simétrica, o que impede a utilização das técnicas largamente
empregadas na resolução de sistemas com matrizes simétricas e em banda.
O MEC tem apresentado excelentes resultados comprovados em numerosos
trabalhos, mostrando ser uma técnica bastante eficiente na solução dos problemas
analisados.
Nota-se, porém, que não existe uma superioridade de determinada técnica sobre
a outra. Isso depende, entre outras coisas, do tipo de problema a ser resolvido.
Assim, combinações dos métodos tem sido objeto de pesquisa, visando o melhor
aproveitamento de cada técnica.
As equações integrais para problemas potenciais e elasticidade apareceram na
literatura no século passado. Em 1872 Betti demonstrou o Teorema da Reciprocidade.
A Identidade de Somigliana, que é a base dos métodos diretos das equações integrais
(BREBBIA,1984), foi apresentada em 1886. Essa identidade é a equação integral que
estabelece a relação entre as forças e deslocamentos no contorno de um corpo e seus
deslocamentos internos.
A formulação de Somigliana é chamada direta. Nela, as variáveis, que são os
deslocamentos e as forças no contorno, têm significado físico. Rizzo, em 1967,
apresentou para a elastostática a formulação que até hoje é usada no MEC,
relacionando deslocamentos e forças. Ele trabalhou com a solução fundamental da
equação diferencial que governa o problema, a qual corresponde à força concentrada.
A outra formulação, chamada indireta, emprega funções fictícias sem qualquer
significado físico. Ela está relacionada ao trabalho publicado por Fredholm em 1903,
que discutiu soluções baseadas no processo da discretização.
7
O estudo da elastodinâmica começou com as teorias de movimento de corpos
elásticos de Navier, em 1827 e Cauchy, em 1828 e continuou com as contribuições de
Green (1838), Stokes (1849) e outros pesquisadores nos problemas de propagação de
ondas.
A solução de problemas que envolvem a variável tempo tem gerado muitas
pesquisas no estudo da formulação adequada. Sendo uma área de atuação
relativamente recente, ainda são muitos os aspectos discutíveis. O número de trabalhos
publicados nos últimos anos, nos quais diferentes formulações são empregadas, mostra
o interesse que o assunto tem despertado, como também, as potencialidades do MEC
nessa área.
Cruse (1968) e Cruse e Rizzo (1968) apresentaram os primeiros trabalhos para
solução de problemas elastodinâmicos, utilizando o MEC. A partir do trabalho de Rizzo
(1967) em elastostática, utilizaram a formulação direta e a Transformada de Laplace
(BESKOS, 1987) para resolver um problema transiente de propagação de onda. Nesse
procedimento, o problema torna-se independente da variável tempo e as equações
integrais são estabelecidas em função dos parâmetros de transformação.
Numa etapa posterior, Manolis e Beskos (1983) aperfeiçoaram essa metodologia
para resolver problemas transientes de dispersão de ondas.
A partir de 1983 a maior parte das pesquisas se desenvolveu usando na sua
formulação soluções fundamentais dependentes do tempo. Destacam-se os trabalhos
de Mansur (1982), que abordou os problemas de propagação de ondas elásticas, e
Mansur e Brebbia (1986) que desenvolveram uma formulação geral no domínio tempo
que usa núcleos bidimensionais. Essa formulação, embora bastante precisa em certos
tipos de problemas, implica em desenvolvimento matemático elaborado e em alguns
casos, considerável esforço computacional.
Nos problemas com cargas aplicadas em todo o corpo, ou em parte dele, pode-
se utilizar discretização do domínio em células (CODA, 1993). Os resultados obtidos
8
com essa técnica são excelentes, mas tal recurso foge um pouco à filosofia do MEC,
que limita a discretização ao contorno.
Para o cálculo da vibração livre de corpos finitos, as poucas opções disponíveis
até 1981 não eram adequadas. Nardini e Brebbia, em 1982, apresentaram uma técnica
que recebeu o nome de Dupla Reciprocidade. Esse procedimento consiste do emprego
de soluções fundamentais independentes do tempo, juntamente com um procedimento
que transforma as integrais de domínio em integrais de contorno. Os problemas de
vibrações livres foram reduzidos à solução algébrica de problemas de autovalores,
como apresentado por Coda (1990). A maior vantagem dessa técnica é que as integrais
do contorno são calculadas apenas uma vez, pois são independentes da freqüência.
Motivados pela propriedade da formulação e pelos ótimos resultados, Nardini e
Brebbia fizeram uma aplicação para obtenção da resposta em problemas de
elastodinâmica transiente, publicada em 1983. A formulação permite resolver esses
problemas empregando a solução fundamental da elastostática na transformação da
integral de domínio em integral de contorno.
O método da dupla reciprocidade é uma maneira geral de se construir soluções
particulares que podem ser usadas para resolver problemas dependentes do tempo,
problemas não-lineares, bem como para representar qualquer distribuição interna de
forças. Podem ser citados entre os inúmeros trabalhos publicados utilizando esse
conceito, os de Loeffler (1988, 1994), Partridge, Brebbia e Wrobel (1992), Coda e
Venturini (1990), Zhu e Zhang (1993), Partridge e Sensale (1997), Chen, Brebbia e
Power (1999), Almeida e Coda (2001) e Barbirato e Venturini (2005).
É interessante notar que, apesar das diversas características vantajosas que
essa formulação oferece, ainda há muitos pontos a serem abordados, visando a melhor
utilização da sua potencialidade como algoritmos auxiliares adequados, recursos que
melhorem a precisão da resposta, entre outros.
9
O Método dos Elementos Finitos (MEF) tem-se constituído numa ferramenta
poderosa e amplamente utilizada para construção de soluções aproximadas para os
mais variados problemas de engenharia. Sua utilização não abrange apenas a
engenharia estrutural, mas também muitas outras áreas, proporcionando significativos
avanços nos mais diversos campos.
Na década de 1940, Hrenikoff (1941), McHenry (1943) e Newmark (1949)
apresentaram trabalhos com a utilização do MEF no tratamento de problemas de
mecânica dos sólidos. Em 1956 foi estabelecida a formulação da forma como é hoje
conhecida, com a publicação de um trabalho de Turner, Clough, Martin e Topp. Os
primeiros livros textos surgiram com Holand e Bell, em 1969, e Zienkiewicz em 1971.
Atualmente existem dezenas de monografias, revistas e livros voltados para o estudo
desse método e esse número está crescendo exponencialmente com revelações
adicionais do poder e da versatilidade de suas aplicações.
Por causa das vantagens e desvantagens que cada método apresenta, muitos
autores têm se interessado no desenvolvimento de algoritmos onde o acoplamento
MEC/MEF é considerado. O primeiro trabalho, que surgiu em 1977, foi de Zienkiewicz,
Kelly e Bess, que utilizaram uma aproximação de energia para o domínio do MEC.
Em se tratando de elastodinâmica, podem ser citados Kobayashi e Kawakami,
em 1985 e Kobayashi e Mori, em 1986 na análise de problemas no domínio freqüência.
Entre os primeiros trabalhos, nos quais formulações de acoplamento para a completa
análise no domínio tempo de problemas transientes são apresentadas, encontram-se
as contribuições de Karabalis e Beskos e Spirakos e Beskos, em 1986, que
investigaram a resposta dinâmica em fundações flexíveis.
Podem-se citar ainda os trabalhos de Von Estorf e Kausel (1989), Coda (1993),
Araújo (1994), Lei e Qinghua (1997), Siqueira (1999), Coda e Venturini (2000), Almeida
(2003) e Almeida, Coda e Mesquita (2004).
10
O presente trabalho apresenta o estudo de problemas elastodinâmicos,
governados por uma equação de campo vetorial, através do MEC. Foram analisados os
problemas de vibrações livres e transientes, com o conceito de matriz de massa,
através de duas técnicas.
A primeira, que desenvolve a integração no domínio por células, e a segunda
através da reciprocidade dual. A discretização do domínio por células foi implementada
para se obter resultados confiáveis nos problemas a serem analisados. Essa técnica
leva a ótimos resultados nesse tipo de análise, embora seja preciso discretizar todo o
domínio.
Na reciprocidade dual, no estudo de regiões onde o contorno apresenta cantos
ou angulosidade, houve a necessidade de se buscar um caminho alternativo na
geração das matrizes de transformação, uma vez que se chegava a uma matriz não
inversível. Isso foi feito através de um algoritmo auxiliar, que modifica o sistema de
equações. Depois de implementada a parte estática, o algoritmo foi adaptado para os
problemas elastodinâmicos de vibrações livres. Com a obtenção de bons resultados,
essa técnica foi levada aos problemas transientes e ao acoplamento MEC/MEF.
1.2 ORGANIZAÇÃO DO TRABALHO
O trabalho está organizado em nove capítulos, sendo o primeiro constituído pela
presente introdução.
No Capítulo 2 apresenta-se uma revisão dos conceitos básicos da teoria da
elasticidade para problemas no regime elástico linear. Demonstra-se a Identidade de
Somigliana para corpo elástico utilizando a técnica dos resíduos ponderados, onde a
função ponderadora é a solução fundamental. A equação integral do contorno é obtida
através do acréscimo de um setor circular e são mostradas as expressões para a
obtenção das tensões.
11
No Capítulo 3 são apresentadas as equações gerais para a elastodinâmica e
para a vibração livre de corpos finitos.
No Capítulo 4 demonstra-se a formulação do Método dos Elementos de
Contorno (MEC). Descreve-se a formulação estática, com as funções de interpolação
utilizadas para a discretização da geometria e das variáveis. Apresentam-se os
elementos disponíveis para a análise e as matrizes dos coeficientes de influência
desses elementos. Demonstra-se a condensação da equação integral de
deslocamentos, proposta para a eliminação das linhas e colunas dependentes. Em
seguida apresenta-se a formulação dinâmica, com a matriz de massa vinda da
integração dos termos inerciais pelos processos da reciprocidade dual e integração
direta. A implementação numérica é descrita com a eliminação dos deslocamentos e
acelerações dependentes.
O Capítulo 5 é dedicado à apresentação da formulação do Método dos
Elementos Finitos (MEF) aplicada às estruturas reticuladas nesse trabalho.
Os métodos de integração numérica para problemas dinâmicos são descritos no
Capítulo 6. Foram utilizados os esquemas Houbolt e Newmark de avanço no tempo.
No Capítulo 7 desenvolveu-se o acoplamento entre os métodos dos elementos
finitos e dos elementos de contorno (MEF e MEC), através da técnica das sub-regiões.
No Capítulo 8 apresentam-se alguns exemplos e uma comparação dos
resultados obtidos com soluções analíticas ou de outras modelagens existentes na
bibliografia para comprovar o programa desenvolvido.
As conclusões e as sugestões para novos estudos e continuidade do trabalho
estão no Capítulo 9.
Nos anexos estão as deduções teóricas complementares, conceitos da notação
indicial cartesiana e alguns detalhes da implementação computacional.
12
CAPÍTULO 2 ELASTOSTÁTICA
13
2.1 INTRODUÇÃO
Apresentam-se, neste capítulo, os conceitos da teoria da elasticidade que serão
utilizados nos capítulos subseqüentes. Estes conceitos serão importantes na aplicação
do Método dos Elementos de Contorno para problemas bidimensionais, no regime
elástico linear.
Admitem-se as hipóteses de pequenos deslocamentos, pequenas deformações e
material obedecendo a Lei de Hooke, isto é, linearidade das relações entre
deformações e tensões, bem como a não influência da mudança de configuração da
estrutura na formulação das equações de equilíbrio.
Para a representação concisa das equações gerais da Teoria da Elasticidade e
dos teoremas dela conseqüentes será utilizada, além da notação usual, a notação
indicial que é uma forma compacta usual de se representar longas expressões. Torna-
se uma alternativa vantajosa que permite uma melhor compreensão das grandezas
envolvidas, bem como facilita a manipulação de somatórios e sistemas de equações.
Tal notação é representada no sistema cartesiano convencional e usa índices 1, 2 e 3
para representar os eixos x, y e z. Algumas regras básicas estão descritas no Anexo 1.
2.2 ESTADO PLANO DE TENSÃO
Seja uma chapa delgada solicitada por forças paralelas ao plano
1 2
x x
(Figura
2.1) e distribuídas uniformemente ao longo da espessura t. Desde que as faces
perpendiculares a
3
x
estejam descarregadas, as tensões
33 13 23
,
e
σ σ σ
são nulas e
pode-se admitir, em princípio, que sejam nulas também no interior da chapa.
Sendo pequena a espessura t, admite-se, aproximadamente, que
11 22 12
,
e
σ σ σ
sejam independentes de
3
x
. (TIMOSHENKO e GOODIER,1980).
14
Figura 2.1 - Chapa delgada submetida a carregamentos no seu plano.
As deformações
13 23
e
ε ε
são nulas e a deformação específica
33
ε
pode ser
calculada em função de
11 22
e
ε ε
conforme equações das componentes de deformação
apresentadas no item 2.5:
11 22
33
( )
1
ν ε ε
ε
ν
+
=
−
(2.1)
onde
ν
é o coeficiente de Poisson.
2.3 ESTADO PLANO DE DEFORMAÇÃO
Seja um corpo prismático, localizado entre dois planos indeslocáveis e sem
atrito, solicitado por forças agindo em planos perpendiculares ao eixo longitudinal
3
x
2
x
2
x
1
x
3
x
1
Q
2
Q
3
Q
4
Q
t
15
(Figura 2.2) e que estas, bem como as condições de contorno, não variam ao longo do
comprimento.
O deslocamento axial
3
u
é nulo nas extremidades e, por simetria, na seção do
meio, podendo admitir que o mesmo ocorra em todas as seções transversais. As
componentes
1 2
u e u
são funções de
1 2
x e x
, portanto
33 13 23
,
e
ε ε ε
são nulas e as
demais independentes de
3
x
.
As tensões
13 23
e
σ σ
são nulas e a tensão longitudinal
33
σ
pode ser
determinada em função de
11 22
e
σ σ
conforme as equações constitutivas apresentadas
no item 2.6:
33 11 22
( )
σ ν σ σ
= +
(2.2)
Figura 2.2 - Sólido prismático submetido a forças paralelas ao plano
1 2
x x
.
2
x
2
x
3
x
1
x
1
Q
2
Q
3
Q
4
Q
16
2.4 EQUAÇÕES DE EQUILÍBRIO
Considere um corpo em equilíbrio sob a ação de um sistema de forças externas
1 2
, , ....,
n
Q Q Q
e, um plano fictício
π
passando através do ponto O no interior desse
corpo, que fica dividido em duas partes, denominadas A e B (Figura 2.3).
Figura 2.3 - Tensões no ponto interno.
A parte A está em equilíbrio com as forças
1 2 3
, ,
Q Q Q
e o efeito da parte B.
Assume-se que esse efeito é continuamente distribuído sobre a seção
π
. Em torno do
ponto O considera-se uma pequena superfície
dA
, e um vetor unitário
ˆ
n
saindo da
superfície e normal a ela. O efeito de B em
dA
pode ser reduzido a um diferencial de
força. Define-se o vetor tensão como (LOVE, 1944):
0
lim
dA
dQ
dA
σ
→
§ ·
=
¨ ¸
© ¹
(2.3)
Pode-se decompor a tensão em duas componentes com relação ao plano: uma
normal, que é a projeção de
σ
na direção da normal
ˆ
n
, e outra tangencial, que é a
ˆ
n
17
projeção de
σ
no plano
π
(Figura 2.4). Essa última ainda pode ser projetada nesse
plano nas duas direções ortogonais.
Figura 2.4 - Decomposição do vetor tensão
σ
.
Forças como
1 2
, ,...,
n
Q Q Q
são chamadas forças de superfície. Forças
distribuídas sobre o volume do corpo, tais como forças gravitacionais e forças
magnéticas são chamadas forças de volume.
Por facilidade, pode-se admitir um sistema de coordenadas cartesianas e as
componentes de tensão
ij
σ
como indicado na Figura 2.5: uma componente normal e
duas tangenciais, todas nas direções dos eixos coordenados. O primeiro índice (i)
indica a direção do eixo perpendicular ao plano em questão, e o segundo (j) indica a
direção da componente de tensão.
Um plano cuja normal externa tem o sentido positivo do eixo é um plano positivo.
A tensão normal na direção dessa normal externa é considerada positiva. A tensão
tangencial de um plano positivo e que tem o sentido positivo do eixo é uma tensão
tangencial positiva.
Já foi visto que, num plano
π
passando por O atua um vetor tensão dado pela
Equação (2.3). Em outro plano passando por O um vetor de tensão diferente irá agir.
Mostra-se que, o vetor de tensão em qualquer plano passando por um ponto pode ser
π
σ
n
σ
t
σ
O
18
obtido quando são conhecidos os vetores de tensão nos três planos normais aos eixos
coordenados, passando por esse ponto, ou seja:
1 11 1 12 2 13 3
2 21 1 22 2 23 3
3 31 1 32 2 33 3
n n n
n n n
n n n
ρ σ σ σ
ρ σ σ σ
ρ σ σ σ
= + +
= + +
= + +
em notação indicial:
i ij j
n
ρ σ
=
(2.4)
onde
i
ρ
é a componente do vetor tensão num plano qualquer e
j
n
são os cossenos
diretores da normal
n
ao plano em questão, com relação ao sistema de coordenadas
definido.
Aplicando a fórmula (2.4) num plano coincidente com a superfície, obtém-se a
condição de equilíbrio no contorno do corpo (Figura 2.5):
1 11 1 12 2 13 3
2 21 1 22 2 23 3
3 31 1 32 2 33 3
p n n n
p n n n
p n n n
σ σ σ
σ σ σ
σ σ σ
= + +
= + +
= + +
em notação indicial:
i ij j
p n
σ
=
(2.5)
19
Figura 2.5 - Representação das forças de superfície no contorno do corpo.
O equilíbrio estático das forças que agem no paralelepípedo mostrado na Figura
2.6 requer que (SAADA, 1994):
13
11 12
1
1 2 3
23
21 22
2
1 2 3
31 32 33
3
1 2 3
0
0
0
b
x x x
b
x x x
b
x x x
σ
σ σ
σ
σ σ
σ σ σ
∂
∂ ∂
+ + + =
∂ ∂ ∂
∂
∂ ∂
+ + + =
∂ ∂ ∂
∂ ∂ ∂
+ + + =
∂ ∂ ∂
em notação indicial:
,
0
ij j i
b
σ
+ =
(2.6)
20
onde
i
b
são as componentes das forças por unidade de volume.
Além disso, se não existe nenhum momento aplicado no corpo, o equilíbrio de
momentos nos leva a igualdade:
12 21
13 31
23 32
σ σ
σ σ
σ σ
=
=
=
ou :
ij ji
σ σ
=
(2.7)
Figura 2.6 – Sólido elementar submetido a esforços internos.
1
dx
2
dx
3
dx
1
x
2
x
3
x
1
b
2
b
3
b
21
2.5 RELAÇÃO DEFORMAÇÃO - DESLOCAMENTO
Um corpo, submetido a um sistema de forças, terá sua configuração inicial
modificada, ou seja, sofrerá uma deformação. Se as componentes de deslocamento
i
u
são tais que sua primeira derivada é tão pequena que quadrados e produtos das
derivadas parciais podem ser desprezados, então a seguinte relação, entre as
componentes do tensor de deformação (
ij
ε
) e do deslocamento (
i
u
) pode ser escrita em
notação indicial (TIMOSHENKO,GOODIER,1980) :
1
, ,
2
( )
ij i j j i
u u
ε
= +
(2.8)
Explicitando os termos:
1
11
1
2
22
2
u
x
u
x
ε
ε
∂
=
∂
∂
=
∂
3
33
3
1 2
12
2 1
3
1
13
3 1
1
2
1
2
u
x
u u
x x
u
u
x x
ε
ε
ε
∂
=
∂
§ ·
∂ ∂
= +
¨ ¸
∂ ∂
© ¹
§ ·
∂
∂
= +
¨ ¸
∂ ∂
© ¹
22
32
23
3 2
1
2
u
u
x x
ε
§ ·
∂
∂
= +
¨ ¸
∂ ∂
© ¹
2.6 RELAÇÃO TENSÃO - DEFORMAÇÃO
O tensor de tensão relaciona-se com o tensor de deformação, para sólidos
elásticos, através da Lei de Hooke. Se o material é isótropo, ou seja, tem as mesmas
propriedades em qualquer direção, essa relação se expressa como (BREBBIA,
TELLES, WROBEL, 1984):
(
)
( )
( )
11 11 22 33 11
22 11 22 33 22
33 11 22 33 33
2
2
2
G
G
G
σ λ ε ε ε ε
σ λ ε ε ε ε
σ λ ε ε ε ε
= + + +
= + + +
= + + +
12 12
13 13
23 23
2
2
2
G
G
G
σ ε
σ ε
σ ε
=
=
=
em notação indicial:
2
ij kk ij ij
G
σ λ ε δ ε
= +
(2.9)
onde
ij
δ
é o delta de Kronecker.
23
As constantes de Lamé,
λ
e
G
, podem ser escritas em função do Módulo de
Young ou de Elasticidade (E) e do Coeficiente de Poisson (
ν
) pelas relações:
2 ( 1 )
E
G
ν
=
+
2
1 2
G
ν
λ
ν
=
−
(2.10)
As equações (2.6), (2.8) e (2.9) formam um conjunto de 15 equações para as 15
incógnitas
,
ij ij
σ ε
e
i
u
e definem o problema de elasticidade linear para estado plano de
deformação. Para o caso de estado plano de tensão deve-se substituir
ν
por
(1 )
ν ν
+
.
2.7 EQUAÇÕES DE NAVIER
Substituindo a equação (2.8) em (2.9), obtém-se a expressão das tensões em
termos das derivadas dos deslocamentos que, substituída na equação de equilíbrio
(2.6) leva à equação de Navier-Cauchy (LOVE, 1944):
( )
( )
( )
2
2 2 2 2 2
31 1 1 1 2
1
2 2 2 2
1 2 3 1 2 1 3 1
2
2 2 2 2 2
3
2 2 2 1 2
2
2 2 2 2
1 2 3 1 2 2 2 3
2 2 2
2
3 3 3
1
2 2 2
1 2 3 1
0
0
uu u u u u
G G b
x x x x x x x x
uu u u u u
G G b
x x x x x x x x
u u u
u
G G
x x x x
λ
λ
λ
§ · ª º
∂∂ ∂ ∂ ∂ ∂
+ + + + + + + =
¨ ¸
« »
∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂
© ¹ ¬ ¼
§ · ª º
∂∂ ∂ ∂ ∂ ∂
+ + + + + + + =
¨ ¸
« »
∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂
© ¹ ¬ ¼
§ ·
∂ ∂ ∂
∂
+ + + +
¨ ¸
∂ ∂ ∂ ∂
© ¹
2
2
3
2
3
2
3 2 3 3
0
u
u
b
x x x x
ª º
∂
∂
+ + + =
« »
∂ ∂ ∂ ∂
¬ ¼
24
em notação indicial:
, ,
( ) 0
j kk k kj j
G u G u b
λ
+ + + =
(2.11)
Tomando-se as equações (2.8) e (2.9), substituindo na equação (2.5) para pontos
no contorno, as forças de superfície são obtidas:
3 3
1 2 1 1 2 1
1 1 1 2 3
1 2 3 1 2 1 3 1
2
u u
u u u u u u
p n G n G n G n
x x x x x x x x
λ
§ · § ·
§ ·
∂ ∂
∂ ∂ ∂ ∂ ∂ ∂
= + + + + + + +
¨ ¸ ¨ ¸
¨ ¸
∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂
© ¹
© ¹ © ¹
3 3
1 2 1 2 2 2
2 1 2 2 3
2 1 1 2 3 2 3 2
3 3 3 3
1 2 1 2
3 1 2 3 3
3 1 3 2 1 2 3 3
2
2
u u
u u u u u u
p G n n G n G n
x x x x x x x x
u u u u
u u u u
p G n G n n G n
x x x x x x x x
λ
λ
§ · § ·
§ ·
∂ ∂
∂ ∂ ∂ ∂ ∂ ∂
= + + + + + + +
¨ ¸ ¨ ¸
¨ ¸
∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂
© ¹
© ¹ © ¹
§ · § · § ·
∂ ∂ ∂ ∂
∂ ∂ ∂ ∂
= + + + + + + +
¨ ¸ ¨ ¸ ¨ ¸
∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂
© ¹ © ¹ © ¹
em notação indicial:
, , ,
( )
i k k i i j j i j
p u n G u u n
λ
= + +
(2.12)
2.8 CONDIÇÕES DE CONTORNO
A equação (2.11) governa o comportamento linear de um corpo elástico
homogêneo e isótropo submetido a ações estáticas. Para esse corpo,
Ω
define o
domínio e
Γ
o contorno. É necessário conhecer as condições de contorno, as quais o
corpo está submetido:
25
a) Condições essenciais ou deslocamentos prescritos:
i i
u u
=
1
m
e
Γ
(2.13)
b) Condições naturais ou forças de superfícies prescritas:
i i
p p
=
2
em
Γ
(2.14)
Nas equações (2.13) e (2.14)
i
u
e
i
p
são os valores conhecidos no contorno. A
superfície externa total do corpo é
1 2
Γ = Γ ∪ Γ
( Figura 2.7).
Essa subdivisão do contorno em duas partes torna possível a resolução dos
casos em que, num mesmo ponto existem os dois tipos de condições de contorno em
diferentes direções ou mesmo uma combinação delas, como no caso de apoios
elásticos.
Figura 2.7 - Domínio
Ω
com as condições de contorno em
Γ
.
a
Q
b
Q
c
Q
26
2.9 EQUAÇÃO INTEGRAL PARA PONTOS DE DOMÍNIO
A equação de equilíbrio da estática (2.6) pode ser aproximada numericamente
através da técnica dos resíduos ponderados. Utilizando-se uma solução fundamental
qualquer como função ponderadora, obtém-se:
(
)
*
,
. 0
ij i j j
b u
σ
+ =
(2.15)
Efetuando uma integral em todo o domínio do corpo, a equação (2.15) passa a
ser representada por:
(
)
*
,
. 0
ij i j j
b u d
σ
Ω
+ Ω =
³
(2.16)
A expressão (2.16) pode ser escrita como:
* *
,
. - .
ij i j j j
u d b u d
σ
Ω Ω
Ω = Ω
³ ³
(2.17)
Integrando-se por partes a primeira integral de (2.17), tem-se:
* * *
, ,
. . . - .
ij i j ij i j ij j i
u d n u d u d
σ σ σ
Ω Γ Ω
Ω = Γ Ω
³ ³ ³
(2.18)
Sendo
*
ij
ε
o campo de deformação correspondente à solução fundamental, a
equação (2.8) e a simetria dos tensores envolvidos, permite escrever:
* * *
,
. . .
ij j i ij ji ij ij
u d d d
σ σ ε σ ε
Ω Ω Ω
Ω = Ω = Ω
³ ³ ³
(2.19)
Substituindo a relação (2.19) na segunda integral do lado direito de (2.18), tem-
se:
27
* * *
,
. . . - .
ij i j ij i j ij ij
u d n u d d
σ σ σ ε
Ω Γ Ω
Ω = Γ Ω
³ ³ ³
(2.20)
Substituindo-se a condição de equilíbrio no contorno, dada pela equação (2.5), em
(2.20):
* * *
,
. . - .
ij i j j j ij ij
u d p u d d
σ σ ε
Ω Γ Ω
Ω = Γ Ω
³ ³ ³
(2.21)
As equações (2.17) e (2.21) têm a mesma integral de domínio no lado esquerdo.
Logo:
* * *
. - . - .
j j ij ij j j
p u d d b u d
σ ε
Γ Ω Ω
Γ Ω = Ω
³ ³ ³
(2.22)
Na expressão (2.22), a segunda integral é de domínio. É preciso transformá-la em
uma integral de contorno.
Sendo
*
ij
σ
o campo de tensão correspondente à solução fundamental, a lei de
Hooke para materiais isótropos faz verdadeira a relação:
* *
. .
ij ij kl kl
σ ε ε σ
=
(2.23)
De forma análoga à equação (2.19):
* *
,
. .
ij ij j i ij
u
ε σ σ
=
(2.24)
Portanto a integral de domínio a ser transformada pode ser escrita como:
* *
,
. .
ij ij j i ij
d u d
σ ε σ
Ω Ω
Ω = Ω
³ ³
(2.25)
Aplicando a integração por partes em (2.25):
28
* * *
, ij ,
. . . - .
j i j i ij j ij i
u d u n d u d
σ σ σ
Ω Γ Ω
Ω = Γ Ω
³ ³ ³
(2.26)
Substituindo-se a condição de equilíbrio no contorno, dada pela equação (2.5) em
(2.26):
* * *
, ,
. . - .
j i ij j j j ij i
u d u p d u d
σ σ
Ω Γ Ω
Ω = Γ Ω
³ ³ ³
(2.27)
Levando-se essa integral para a equação (2.22), obtém-se:
* * * *
,
. - . . . 0
j j j j j ij i j j
p u d u p d u d b u d
σ
Γ Γ Ω Ω
Γ Γ + Ω + Ω =
³ ³ ³ ³
(2.28)
O campo de tensão correspondente à solução fundamental satisfaz a equação de
equilíbrio:
* *
,
0
ij i
j
b
σ
+ =
(2.29)
E a equação (2.28) fica escrita como:
* * * . *
. - . - . . 0
j j j j j j j j
p u d u p d b u d b u d
Γ Γ Ω Ω
Γ Γ Ω + Ω =
³ ³ ³ ³
(2.30)
Torna-se necessário modificar a equação (2.30) para uma forma integral a fim de
que se possa realizar uma análise numérica sobre ela. Será analisada a terceira
integral da equação, que é a equação de domínio:
* .
.
j j
b u d
Ω
Ω
³
(2.31)
29
Adotando-se as componentes das forças de volume
*
j
b
como forças concentradas
unitárias aplicadas no ponto
*
s
∈ Ω
em cada uma das três direções ortogonais
definidas pelo vetor de componentes
j
p
, tem-se:
*
( , ) .
j
j
b s q p
δ
=
(2.32)
onde
1
j
p
=
e
( , )
s q
δ
é a função Delta de Dirac.
A distribuição Delta de Dirac é muito importante para a formulação do Método
dos Elementos de Contorno. Seja
s
o ponto onde as forças unitárias serão aplicadas,
q
o ponto onde as respostas a essas forças serão avaliadas e
f
uma função contínua
qualquer. A função tem as seguintes propriedades:
( , )
s q
δ
= ∞
se q = s
( , ) 0
s q
δ
=
se q
≠
s (2.33)
{ }
( ) se
( ) ( , )
0 se
f s s
f q s q d
s
δ
∞
Ω
∈ Ω
°
Ω =
®
∈ Ω − Ω − Γ
°
¯
³
Substituindo (2.33) em (2.31) , tem-se:
* .
. ( ) .
j j j j
b u d u s p
Ω
Ω =
³
(2.34)
A equação (2.33) considera as cargas unitárias aplicadas ao mesmo tempo. Se
cada uma delas atuar independente da outra, os deslocamentos e forças de superfície
podem ser escritos como:
30
* *
* *
( , )
( , )
j ij i
j ij i
u u s q p
p p s q p
=
=
(2.35)
onde
* *
( , ) ( , )
ij ij
u s q e p s q
representam os deslocamentos e forças de superfície na
direção
j no ponto q, devido à uma força unitária aplicada na direção i no ponto s.
De maneira geral, sendo
ij
δ
o Delta de Kronecker, a equação (2.32) fica:
*
. ( , )
ij
ij
b s q
δ δ
=
(2.36)
Substituindo o estado de deslocamento
*
ij
u
e seu carregamento correspondente
*
ij
b
na equação (2.30), tem-se a equação integral do contorno para deslocamentos:
( )
* *
*
( ) ( , ) - ( ) ( , )
( ) ( , )
i j ij j ij
j ij
u s p Q u s Q d u Q p s Q d
b q u s q d
Γ Γ
Ω
= Γ Γ +
+ Ω
³ ³
³
(2.37)
onde
s
é o ponto fonte,
Q
é ponto do contorno e
q
é ponto do domínio.
Essa equação é conhecida como a Identidade de Somigliana para
deslocamentos e representa o deslocamento
i
u
em qualquer ponto
q
do domínio. Ela
pode ser escrita na forma de equilíbrio da equação de Navier - Cauchy (2.11), como:
* *
, ,
( ) ( , ) 0
j kk k kj j
G u G u s q p
λ δ
+ + + =
(2.38)
As soluções para essa equação são chamadas soluções fundamentais.
31
2.10 SOLUÇÃO FUNDAMENTAL
Utilizam-se as soluções fundamentais, obtidas de problemas específicos, para
desenvolver formulações do Método dos Elementos de Contorno (LOVE, 1944;
BREBBIA, WALKER, 1980; BREBBIA, TELLES, WROBEL, 1984). Seja
*
Ω
um meio
elástico infinito e, por conseqüência,
*
Γ
o contorno infinito. Esse caso corresponde à
solução fundamental de Kelvin, que representa fisicamente o efeito de uma carga
estática unitária concentrada atuando no ponto
s
∈
*
Ω
(Figura 2.8). Essa é a solução
para a equação de equilíbrio de Navier - Cauchy (Equação(2.38)).
Figura 2.8 - Carga concentrada unitária aplicada na região infinita.
As expressões das soluções fundamentais são:
a) Deslocamentos
2
x
1
x
*
Γ
* *
,
ij ij
u p
1
1
Γ
s
q
1
r
2
r
*
Ω
Ω
32
*
, ,
-1
( , ) (3 - 4 ) ln( )
8 (1- )
ij ij i j
u s q r r r
G
ν δ
π ν
ª º
= −
¬ ¼
(2.39)
O primeiro índice corresponde à direção de aplicação da carga unitária, o
segundo à direção do deslocamento e r = r (s, q) à distância entre os pontos
q
e
s
(Figura 2.8).
b) Forças de Superfície
Introduzindo a equação (2.39) em (2.17), obtêm-se as forças de superfície nas
direções j, devidas à carga unitária em
s
na direção i. Essas forças atuam em uma
superfície definida pelos cossenos diretores (
j
n
) de um vetor normal a ela:
*
, , , ,
-1
( , ) (1 - 2 ) 2 - (1 - 2 )( - )
4 (1- )
ij ij i j i j j i
r
p s q r r r n r n
r n
ν δ ν
π ν
½
∂
ª º
= +
® ¾
¬ ¼
∂
¯
¿
(2.40)
c) Deformações Específicas
Considerando
*
kji
ε
(Figura 2.9) como a deformação
*
ji
ε
em um ponto q, devido a
uma carga unitária aplicada no ponto
s
na direção k, a relação deformação -
deslocamento (equação (2.8)) fornece:
*
, , , , , ,
- 1
( , ) (1- 2 )( ) - 2
8 (1- )
kji i kj j ki k ji i j k
s q r r r r r r
G r
ε ν δ δ δ
π ν
ª º
= + +
¬ ¼
(2.41)
33
Figura 2.9 - Deformações específicas fundamentais.
d) Tensões
Seja
*
kji
σ
(Figura 2.10) a tensão
*
ji
σ
em um ponto
q
, devido a uma carga unitária
aplicada no ponto
s
na direção k. A equação (2.41) em (2.9) fornece:
*
, , , , ,
- 1
( , ) (1- 2 ) 2
4 (1- )
kji i kj j ki i j k
s q r r r r r
r
σ ν δ δ
π ν
ª º
= + +
¬ ¼
(2.42)
q
q
s
s
1
p
=
1
p
=
r
r
*
122
ε
*
121
ε
*
111
ε
*
111
ε
*
112
ε
*
122
ε
ponto
q
*
222
ε
*
221
ε
*
211
ε
*
211
ε
*
212
ε
*
222
ε
ponto
q
34
Figura 2.10 - Tensões específicas fundamentais.
Para todas as expressões, as derivadas de
( , )
r r s q
=
são tomadas com
referência às coordenadas do ponto
q
, ou seja:
,
( )
j
j
j
r
r
r
x q r
∂
= =
∂
(2.43)
As soluções fundamentais acima são válidas para o estado plano de
deformação. No caso do estado plano de tensão, deve-se substituir
ν
por
1 (1 )
ν ν
= +
.
2.11 EQUAÇÃO INTEGRAL PARA PONTOS DO CONTORNO
Considerando a solução fundamental de Kelvin, a Identidade de Somigliana
(equação (2.37)) só é satisfatória para a obtenção de deslocamentos em qualquer
q
q
s
s
1
p
=
1
p
=
r
r
*
122
σ
*
121
σ
*
111
σ
*
111
σ
*
112
σ
*
122
σ
ponto
q
*
222
σ
*
221
σ
*
211
σ
*
211
σ
*
212
σ
*
222
σ
ponto
q
35
ponto interno
s
(Figura 2.11), quando os deslocamentos e forças no contorno são
conhecidos.
Figura 2.11 - Ponto fonte interior.
É preciso determinar os valores incógnitos de deslocamentos e forças de
superfície de todos os pontos do contorno. Para isso escreve-se a equação integral
para pontos do contorno utilizando um artifício que transforma, inicialmente, o ponto de
contorno em um ponto do domínio, sobre o qual se pode aplicar a Identidade de
Somigliana (CODA, 1990).
Nessa técnica o corpo representado como na Figura 2.12 tem o domínio
Ω
acrescido de um setor de círculo com o centro no ponto fonte
S
de raio
ε
. O ponto
S
do
contorno agora pertence ao domínio
ε
Ω + Ω
.
n
Γ
Ω
s
q
36
Figura 2.12 – Contorno expandido por superfície esférica.
O novo domínio passa a ser
ε
Ω + Ω
e o seu contorno
-
ε
Γ Γ + Γ
, para os quais a
equação (2.37) pode ser aplicada:
( )
* *
*
( ) ( , ) - ( ) ( , )
( ) ( , )
i j ij j ij
j ij
u S p Q u S Q d u Q p S Q d
b q u S q d
ε ε
ε
Γ−Γ+Γ Γ−Γ+Γ
Ω+Ω
= Γ Γ +
+ Ω
³ ³
³
(2.44)
onde :
Γ
= contorno da superfície que foi expandida
ε
Γ
= contorno da superfície esférica acrescida
ε
Ω
= domínio da parte esférica acrescida
ε
Ω
ε
Γ
Γ
S
Q
q
Ω
Γ
37
Separando-se as integrais para cada trecho do domínio e do contorno, escreve-
se:
* *
* *
* *
( ) ( ) ( , ) - ( ) ( , )
( ) ( , ) ( ) ( , ) -
- ( ) ( , ) ( ) ( , )
i j ij j ij
j ij j ij
j ij j ij
u S p Q u S Q d u Q p S Q d
b q u S q d p Q u S Q d
u Q p S Q d b q u S q d
ε
ε ε
Γ−Γ Γ−Γ
Ω Γ
Γ Ω
= Γ Γ +
+ Ω + Γ
Γ + Ω
³ ³
³ ³
³ ³
(2.45)
Para que o ponto
S volte a pertencer ao contorno é necessário fazer o limite de
ε
tender a zero, para que
εε
ΩΓ e
tendam a zero também.
As integrais com núcleos semelhantes a
*
ij
u
(Equação (2.45)), são chamadas
integrais singulares de camada simples e demonstra-se que para
0
ε
→
o seu limite é
nulo. Desta forma, tem-se:
*
0
*
0
lim ( ) . ( , ) 0
lim ( ) . ( , ) 0
j ij
j ij
b q u S q d
p Q u S Q d
ε
ε
ε
ε
→
Ω
→
Γ
Ω =
Γ =
³
³
(2.46)
No trecho do contorno
-
Γ Γ
, desde que a integral do lado direito da igualdade
seja calculada de acordo com o conceito de valor principal de Cauchy, pode-se
escrever que:
* *
0
lim ( ) . ( , ) ( ) . ( , )
j ij j ij
p Q u S Q d p Q u S Q d
ε
→
Γ−Γ Γ
Γ = Γ
³ ³
(2.47)
38
Já as integrais com núcleo
*
ij
p
são chamadas de camada dupla e espera-se uma
descontinuidade no limite
0
ε
→
. O desenvolvimento do limite sobre
ε
Γ
para esse
termo pode ser feito:
* *
0 0
*
0
lim ( ) ( , ) lim ( ) - ( ) ( , )
lim ( ) ( , )
j ij j j ij
j ij
u Q p S Q d u Q u S p S Q d
u S p S Q d
ε ε
ε ε
ε ε
ε
ε
ε
→ →
Γ Γ
→
Γ
ª º
Γ = Γ +
¬ ¼
+ Γ
³ ³
³
(2.48)
O primeiro termo do lado direito da equação (2.48) é identicamente nulo, uma vez
que o campo de deslocamento
( )
j
u Q
satisfaz a condição de H
o
lder :
( ) - ( ) ( , )
k k
u Q u S B r S Q
α
≤
(2.49)
sendo B e
α
constantes positivas.
Assim, na expressão (2.48) o limite sobre
ε
Γ
é igual a:
* *
0 0
lim ( ) ( , ) lim ( ) ( , ) ( ) ( )
j ij j ij ij j
u Q p S Q d u S p S Q d C S u S
ε ε
ε ε
ε ε
→ →
Γ Γ
Γ = Γ =
³ ³
(2.50)
O coeficiente
ij
C
é dado pela seguinte expressão:
( )
ij ij ij
C S I
δ
= +
(2.51)
onde:
*
0
lim ( , )
ij ij
I p S Q d
ε
ε
→
Γ
= Γ
³
(2.52)
39
com
ij
δ
, que é o Delta de Kronecker.
Demonstra-se que
( )
2
ij
ij
C S
δ
= para contornos suaves.
Se o ponto
S
pertence ao domínio, tem-se
( ) 1
ij
C S
=
, e para pontos externos ao
corpo, tem-se
( ) 0
ij
C S
=
.
Para pontos definidos em contornos com angulosidade, o limite da equação (2.52)
pode ser expresso como função dos ângulos
1 2
e
θ θ
definidos na Figura 2.13:
1 2
-1
8 ( 1 - ) 2 3
ij
A A
I
A A
π ν
ª º
=
« »
¬ ¼
com:
2 1 1 2
2 1
2 1 2 1
1 4 ( 1 - ) ( - ) 2 - 2
2 cos 2 - cos 2
3 4 ( 1 - ) ( - ) 2 - 2
A sen sen
A
A sen sen
ν π θ θ θ θ
θ θ
ν π θ θ θ θ
= + +
=
= + +
(2.53)
Figura 2.53 - Definição dos ângulos para cálculo dos termos da matriz C.
1
n
2
n
1
θ
2
θ
40
A equação integral do contorno fica dada por:
* *
*
( ) ( ) ( ) ( , ) ( , ) ( )
( ) ( , )
ij j j ij ij j
j ij
C S u S u Q p S Q d u S Q p Q d
b q u S q d
Γ Γ
Ω
+ Γ = Γ +
+ Ω
³ ³
³
(2.54)
Essa equação fornece uma relação entre os deslocamentos e as forças de
superfície que deve ser satisfeita. Introduzindo as condições de contorno, obtêm-se as
incógnitas apenas no contorno.
2.12 REGIÕES INFINITAS
As equações integrais analisadas até agora levam em conta apenas corpos finitos.
A extensão da equação (2.54) para regiões infinitas com uma ou mais cavidades
internas, requer uma análise cuidadosa do comportamento das funções envolvidas.
Essa análise está relacionada ao comportamento das funções sobre uma superfície de
contorno infinitamente distante das cavidades.
Seja
r
o raio de uma esfera de superfície
r
Γ
, centrada em
S
, que envolve as
cavidades do problema representado na Figura 2.14. A equação (2.54) pode ser escrita
para esse corpo com
r
finito, como:
* *
* *
( ) ( ) ( ). ( , ) ( ). ( , )
( , ) . ( ) ( , ) . ( )
ij j j ij j ij
r
ij j ij j
r
C S u S u Q p S Q d u Q p S Q d
u S Q p Q d u S Q p Q d
Γ Γ
Γ Γ
+ Γ + Γ =
= Γ + Γ
³ ³
³ ³
(2.55)
41
Figura 2.64 - Definição da região infinita com cavidade interna
No limite, quando
r
→ ∞
, a equação (2.55) apresenta o significado desejado se
for satisfeita a condição de regularidade:
* *
lim ( ) . ( , ) - ( , ) ( ) 0
j ij ij j
r
r
u Q p S Q u S Q p Q d
→∞
Γ
ª º
Γ =
¬ ¼
³
(2.56)
Para problemas bidimensionais, considerando o contorno infinito (
r
Q
∈ Γ
), tem-
se:
*
*
onde
( ) ;
( ln 1)
( , )
(1)
1
( , )
ij
ij
d G d G O r
O r i j
u S Q
O i j
p S Q O
r
ϕ
Γ = =
+ =
=
®
≠
¯
§ ·
=
¨ ¸
© ¹
(2.57)
onde O ( ) é o comportamento assintótico de uma função para
0
r ou r
→ ∞ →
.
r
S
Γ
Ω
Q
infinito
Γ
42
Se a carga total aplicada sobre a superfície
Γ
não for auto-equilibrada, o princípio
de Saint-Venant mostra que
( )
j
u Q
e
( )
j
p Q
terão o mesmo comportamento da solução
fundamental correspondente a uma carga concentrada na direção resultante. Logo,
( ) 0 (ln 1) ( ) 0 (1 )
j j
u Q r e p Q r
= + =
são obtidos, o que não garante, em geral, a
anulação de cada termo separadamente. Mas, pode-se substituir
( ) ( )
j j
u Q e p Q
pelos
tensores correspondentes à solução fundamental e verificar que a equação (2.55) é
satisfeita, pois os termos se cancelam quando
r
→ ∞
.
Pode-se afirmar que as condições de regularidade são sempre satisfeitas se
( ) ( )
j j
u Q e p Q
se comportam, na pior das hipóteses, como a solução fundamental no
infinito.
Neste caso, os problemas de cavidade em meio infinito podem ser representados
pela equação (2.54), com a normal apontando para dentro da cavidade.
Figura 2.15 - Definição da normal.
n
n
Γ
43
2.13 TENSÕES NOS PONTOS INTERNOS
É necessário conhecer as componentes de tensão em qualquer ponto do sólido
em estudo, para que se possa considerar resolvido o problema elastostático.
A equação (2.37), que é uma representação contínua dos deslocamentos nos
pontos internos, pode ser derivada em relação às coordenadas do ponto
S
.
Substituindo as derivadas nas equações (2.8) e (2.9), obtém-se a expressão das
tensões para pontos interiores (CODA, 1990):
* *
*
( ) ( , ) ( ) - ( , ) ( )
( , ) ( )
ij kij k kij k
kij k
s D s Q p Q d S s Q u Q d
D s q b q d
σ
Γ Γ
Ω
= Γ Γ +
+ Ω
³ ³
³
(2.58)
Para o estado plano de deformação, tem - se:
( )
*
, , , , , ,
*
, , ,
2
, , , , , ,
1
1- 2 - 2
4 (1- )
2 (1- 2 ) ( ) -
2 (1- )
- 4 2 (
kij ki j kj i ij k i j k
kij ij k ik j jk i
i j k i j k j
D r r r r r r
r
G r
S r r r
r n
r r r n r r n r
ν δ δ δ
π ν
ν δ ν δ δ
π ν
ν
ª º
= + +
¬ ¼
∂
ª
= + +
®
¬
∂
¯
º
+ +
¼
,
, ,
) (1- 2 )
(2 ) - (1- 4 ) )
i k
k i j j ik i jk k ij
r
n r r n n n
ν
δ δ ν δ
+
½
+ +
¾
¿
(2.59)
Para o estado plano de tensão substitui–se
por 1 (1 )
ν ν
+
.
44
2.14 TENSÕES NOS PONTOS DO CONTORNO
As tensões em pontos do contorno podem ser calculadas com bastante precisão
através dos valores conhecidos das forças de superfície e dos deslocamentos, após a
resolução do problema. Esses valores são obtidos no sistema global de coordenadas.
Utilizando a matriz de transformação para tensor de primeira ordem, determinam-
se os deslocamentos e forças de superfície em relação a um sistema de referência local
ao elemento considerado (Figura 2.16).
Figura 2.16 - Sistema de referência local ao elemento.
A equação de Cauchy (2.4) para o sistema local, fornece:
22 2
12 1
p
p
σ
σ
=
=
(2.60)
Sendo
1
u
o deslocamento interpolado ao longo do elemento em função dos
valores nodais, a componente de deformação para esse sistema é dada por:
1
11
1
u
x
ε
∂
=
∂
(2.61)
2
x n
=
1
x
1
x
2
x
Ω
45
A lei de Hooke (Equação (2.9)) permite escrever:
]
22
22 11
1
(1- 2 ) -
(1- 2 ) 2G
σ
ε ν ν ε
ν
ª
=
«
¬
(2.62)
Com os valores de
11 22
e
ε ε
pode-se obter as tensões em um ponto qualquer do
contorno em relação ao sistema local do elemento, através da lei de Hooke
generalizada.
Assim:
11 11 22
2 -
1-
G
ν
σ ε ε
ν
ª º
=
« »
¬ ¼
(2.63)
As equações (2.63) e (2.60) fornecem as tensões em um ponto qualquer do
contorno em relação ao sistema local de coordenadas.
46
CAPÍTULO 3 ELASTODINÂMICA
47
3.1 INTRODUÇÃO
Nesse capítulo as expressões estão também representadas sob a forma de
notação indicial e no sistema de coordenadas cartesianas convencionais. Serão
apresentados os conceitos necessários para o desenvolvimento da formulação da
vibração livre e elastodinâmica transiente bidimensional.
3.2 EQUAÇÕES GERAIS DA ELASTODINÂMICA
Admitem-se novamente as hipóteses de pequenos deslocamentos, pequenas
deformações e material obedecendo à Lei de Hooke. O corpo em estudo é elástico,
homogêneo e isótropo.
A condição de equilíbrio no contorno de um corpo é dada por:
i ij j
p n
σ
=
(3.1)
A relação entre as componentes do tensor de deformação e do deslocamento é
escrita como:
( )
, ,
1
2
ij i j j i
u u
ε
= + (3.2)
A Lei de Hooke relaciona tensor tensão com tensor deformação:
2
ij kk ij ij
G
σ λ ε δ ε
= +
(3.3)
com G e
λ
dados pela equação (2.10).
O equilíbrio dinâmico das forças agindo no paralelepípedo mostrado na Figura
(3.1) requer que (MANOLIS, 1988):
48
,
ij j i i
b u
σ ρ
+ =
(3.4)
onde:
i
b
são as forças de volume;
ρ
é a densidade do corpo;
2
2
i
i
u
u
t
∂
=
∂
é a aceleração com relação à direção
i
x
.
Figura 3-1 - Sólido elementar submetido a esforços internos.
As equações de equilíbrio de forças resultantes de tensões (3.4), relações
cinemáticas (3.2) e Lei de Hooke (3.3) formam um conjunto de 15 equações para 15
i i i
b b u
ρ
= −
1
dx
2
dx
3
dx
1
x
2
x
3
x
49
incógnitas (
, ,
ijj ij i
u
σ ε
). Substituindo a expressão (3.2) na expressão (3.3), a Lei de
Hooke é escrita em termos dos deslocamentos:
, , ,
( )
ij k k ij i j j i
u G u u
σ λ δ
= + + (3.5)
As equações (3.5) e (3.4) nos levam à equação geral de Navier-Cauchy:
, ,
( )
i jj j ji i i
G u G u b u
λ ρ
+ + + =
(3.6)
A equação (3.6) constitui um sistema linear de equações diferenciais hiperbólicas
para a variável
i
u
.
As ondas que ocorrem no estado elastodinâmico, para o qual as rotações são
nulas, são chamadas ondas longitudinais e se propagam com uma velocidade dada
por:
( 2 )
d
c G
λ ρ
= +
(3.7)
Quando a variação do volume é nula (
,
i i
u
= 0), as ondas num corpo elástico são
chamadas distorcionais e se propagam com uma velocidade
s
c
, dada por:
ρ
G = c
s
(3.8)
As constantes
d s
c e c
podem ser usadas para escrever a equação de
equilíbrio em termos dos deslocamentos, isto é, a equação de Navier-Cauchy:
2 2 2
, .
( - )
s i jj d s j ji i i
c u c c u b u
ρ
+ + =
(3.9)
Para resolver o problema elastodinâmico isotrópo é necessário determinar as
componentes de deslocamentos
( , )
i
u S t
que satisfaçam:
50
a) a equação (3.9) para
0
t t
≥
em todos os pontos do domínio
Ω
;
b) as condições iniciais:
0 0
0 0
0
( , ) ( )
( , )
( , ) ( )
i i
i
i i
t t
u s t u s
u s t
u s t v s
t
=
°
=
°
°
®
°
∂
ª º
°
= =
« »
∂
°
¬ ¼
¯
(3.10)
em todos os pontos s pertencentes ao domínio. O índice ‘0’ significa que os valores
prescritos correspondem ao tempo
0
0
t
=
. O ponto sobre a variável representa derivada
parcial com relação ao tempo e v é a velocidade.
Figura 3-2- Deslocamentos e velocidades prescritas em t = 0.
c) Condições de Contorno:
1
2
( , ) ( , )
( , ) ( , )
i i
i ij j i
u S t u S t S
p S t n p S t S
σ
= ∈ Γ
°
®
°
= = ∈ Γ
¯
(3.11)
onde
i
p
são os valores conhecidos das forças em
2
Γ
,
i
u
são os valores conhecidos
dos deslocamentos em
1
Γ
, S são pontos do contorno e t é tempo.
(
)
( )
v
io io
io io
u s
u u s
=
=
Ω
Γ
51
Figura 3-3- Condições de contorno
i i
p e u
no tempo t em
1 2
e
Γ Γ
.
Reordenando as equações (3.7) e (3.8) para obter
e G
λ
em função das
velocidades
d
c
e
s
c
das ondas e substituindo na equação (3.5), as tensões podem ser
também escritas como:
= + +
2 2 2
, , ,
( - 2 ) ( )
ij d s k k ij s i j j i
c c u c u u
σ ρ δ ρ
(3.12)
A segunda das condições de contorno dada pela equação (3.11) é escrita em
função dos deslocamentos usando a equação (3.1):
= + +
2 2 2
, , ,
( - 2 ) ( )
i d s k k i s i j j i j
p c c u n c u u n
ρ ρ
(3.13)
3.3 IDENTIDADE DE SOMIGLIANA
O procedimento para se obter a Identidade de Somigliana é idêntico ao
apresentado para problemas elastostáticos, utilizando a solução fundamental de Kelvin,
n
d
Γ
s
S
1
Γ
2
Γ
Ω
52
que representa fisicamente o efeito de uma carga unitária estática atuando em um
domínio infinito. A versão correspondente ao problema tratado aqui pode ser obtida
diretamente, pela substituição da componente da força de volume
i
b
por
i
u
ρ
na
equação (2.37).
Assim:
* *
( , ) ( , ) ( , ) ( ) - ( , ) ( , ) ( )
k i ki i ki
u s t p Q t u s Q d Q u Q t p s Q d Q
Γ Γ
= Γ Γ +
³ ³
*
( , ) ( , )
ki i
u s q u Q t
ρ
Ω
+
³
(3.14)
onde, com exceção dos tensores fundamentais, todas as variáveis envolvidas na
análise são também função da variável tempo.
A formulação baseada na equação (3.14) é recomendada devida a sua
simplicidade e universalidade. O problema está no fato de que, para essa formulação, a
integral inercial é estendida a todo o domínio, o que acarreta um esforço adicional de
discretização.
3.4 EQUAÇÃO INTEGRAL DE CONTORNO
A Identidade de Somigliana (3.14) é válida para pontos
s
contidos no interior do
corpo (
Ω
). Quando o ponto fonte está no contorno (
Γ
) a obtenção de uma expressão
similar é elaborada através do processo limite apresentado na seção 2.11 . A equação
integral fica:
53
* *
( ) ( ) ( , ) ( , ) - ( , ) ( , )
ij j ij j ij j
C S u S u S Q p Q t d p S Q u Q t d
Γ Γ
= Γ Γ +
³ ³
*
( , ) ( , )
ij j
u S q u q t d
ρ
Ω
+ Ω
³
(3.15)
3.5 EQUAÇÃO INTEGRAL PARA PONTOS EXTERNOS
A equação integral para pontos externos é muito útil em problemas com
descontinuidade de reações de contorno, ocasionada por deslocamentos prescritos.
Considera-se o campo virtual de deslocamentos como causado por uma carga unitária
que não esteja contida na região ocupada pelo corpo estudado. A equação será:
* *
*
0 ( , ) ( , ) - ( , ) ( , )
( , ) ( , )
ij j j ij
ij j
u s Q p Q t d u Q t p s Q d
u S q u q t d
ρ
Γ Γ
Ω
= Γ Γ +
+ Ω
³ ³
³
(3.16)
3.6 EQUAÇÃO INTEGRAL DAS COMPONENTES DE TENSÃO
A expressão para se calcular as tensões nos pontos internos do sólido em
estudo é obtida substituindo o valor de (equação (3.14) ) na relação tensão deformação
dada pelas equações (2.9) e (2.8):
54
*
*
*
( , ) ( , ) ( , ) ( ) -
- ( , ) ( , ) ( )
( , ) ( , ) ( )
ij ijk k
ijk k
ijk k
s t u s q p q t d q
p s q u q t d q
u s q u q t d q
σ
ρ
Γ
Γ
Ω
= Γ
Γ +
+ Γ
³
³
³
(3.17)
onde os valores de são dados pela equação (2.54).
3.7 VIBRAÇÕES LIVRES DE CORPOS FINITOS
A equação diferencial que governa a vibração livre de um corpo elástico, linear,
homogêneo e isótropo pode ser escrita como:
,
ij j i i
b u
σ ρ
+ =
(3.18)
É um caso de movimento governado pela equação de Navier-Cauchy, com
condições de contorno homogêneas. A equação geral em termo dos deslocamentos,
suprimindo as forças de volume, fica:
, ,
( )
i jj j ji i
G u G u u
λ ρ
+ + =
(3.19)
Uma alternativa para se obter as freqüências naturais e os modos de vibração de
uma estrutura, (BREBBIA e NARDINI, 1982) é a utilização da solução fundamental de
Kelvin, apresentada no item 2.10.
O problema se reduz ao estudo de autovalor e autofunção a ele associado, e, a
equação do movimento (3.19) se simplifica para:
2
, ,
( )
i jj j ji i
G u G u u
λ ρ ω
+ + =
(3.20)
55
onde
ω
é a freqüência natural.
As condições de contorno são:
1
2
( , ) 0
( , ) 0
i
i ij j
u S t S
p S t n S
σ
= ∈ Γ
= = ∈ Γ
(3.21)
e as condições iniciais:
0 0
0 0 0
0
( , )
( , ) ( , ) , 0
i i
i i i
t t
u s t u
u s t u s t v s t
t
=
=
∂
ª º
= = ∈ Ω =
« »
∂
¬ ¼
(3.22)
A equação integral fica:
* *
( ) ( ) ( , ) ( , ) - ( , ) ( , ) -
ki i ki i ki i
C S u S u S Q p Q t d p S Q u Q t d
Γ Γ
= Γ Γ
³ ³
2 *
( , ) ( )
ki
u S q d q
ω ρ
Ω
− Ω
³
(3.23)
56
CAPÍTULO 4 MÉTODO DOS ELEMENTOS DE
CONTORNO (MEC)
57
4.1 FORMULAÇÃO ESTÁTICA
4.1.1 INTRODUÇÃO
A equação integral de contorno para problemas estáticos (2.54) pode ser escrita
como:
* *
*
( ) ( ) ( , ) ( ) ( ) - ( , ) ( ) ( )
( ) ( , ) ( )
ki i ki i ki i
i ki
C S u S u S Q p Q d Q p S Q u Q d Q
b q u S q d q
Γ Γ
Ω
= Γ Γ +
+ Ω
³ ³
³
(4.1)
A obtenção de soluções analíticas fechadas, para a maioria dos problemas no
campo da engenharia, apresenta grandes dificuldades. Para evitar que os problemas
possíveis de serem resolvidos se restrinjam aos que se consegue a solução analítica,
transforma-se a equação (4.1) num conjunto de equações algébricas que possam ser
resolvidas por procedimentos numéricos.
Essa transformação envolve a discretização do contorno em uma série de
elementos sobre os quais os deslocamentos e forças de superfície são interpolados por
funções polinomiais associadas a certo número de nós, ou pontos nodais, do elemento.
Os valores ligados aos nós são chamados valores nodais.
A equação (4.1) é escrita, em forma discretizada, para cada ponto fonte aplicado
nos pontos funcionais, e assim as integrais sobre cada elemento são determinadas.
Conseqüentemente o sistema de equações algébricas com valores nodais de
deslocamentos e forças de superfície é obtido para a estrutura.
58
Figura 4-1- Discretização do contorno.
As condições de contorno e as equações de compatibilidade e equilíbrio nos
pontos nodais são impostas e o sistema pode então ser resolvido. Com a solução do
sistema são obtidos os deslocamentos e forças de superfície no contorno. Isso
possibilita os cálculos dos valores dos deslocamentos e tensões em quaisquer pontos
do domínio.
A Figura 4-2 ilustra os elementos que serão utilizados nesse trabalho.
1
x
2
x
Ω
nós
n
Γ
j
Γ
representação
geométrica do
contorno
59
Figura 4-2 - Tipos de elementos de contorno.
Para facilitar as integrações numéricas sobre os elementos, foram
parametrizadas as coordenadas de cada ponto do elemento com relação aos nós
extremos (elemento linear) ou nós extremos e nó central (elemento quadrático) através
de coordenadas locais homogêneas.
4.1.2 DISCRETIZAÇÃO GEOMÉTRICA E DAS VARIÁVEIS
4.1.2.1 ELEMENTO LINEAR
As coordenadas cartesianas x dos pontos do contorno situados ao longo do
elemento
Γ
j
são expressas em termos das coordenadas dos nós e das funções de
interpolação
ϕ
:
elemento quadratico
elemento linear
60
( )
T N
x x
ϕ η
=
(4.2)
onde:
1
2
x
x
x
½
=
® ¾
¯ ¿
1 2
1 2
0 0
0 0
T
ϕ ϕ
ϕ
ϕ ϕ
ª º
=
« »
¬ ¼
(4.3)
{
}
1 2 1 1
1 1 2 2
T
N
x x x x x=
Figura 4-3 - Descrição geométrica do elemento linear.
As funções aproximadoras
)ij(Ș
são dadas por:
elemento
1
x
1
1
x
2
1
x
1
2
x
2
2
x
2
x
1
nó
2
nó
1
η
=
1
η
= −
η
l
61
[ ]
1
2
1
(1- )
2
-1 , 1
1
(1 )
2
ϕ η
η
ϕ η
=
∈
= +
(4.4)
A aproximação discretizada das variáveis, que são os deslocamentos e forças de
superfície dos pontos do contorno situados ao longo do elemento, é feita da mesma
forma que a discretização geométrica. Eles são expressos em termos dos
deslocamentos e forças de superfície dos pontos nodais, relacionados pelas funções de
interpolação.
Figura 4-4 - Deslocamentos e forças de superfície.
Assim:
T N
u u
ϕ
=
T N
p p
ϕ
=
(4.5)
1
x
2
x
l
1
η
=
1
η
= −
2 2
valores nodais de u ou p
1 1
valores nodais de u ou p
1 1 1 1
u ou p
ϕ ϕ
2 2 2 2
u ou p
ϕ ϕ
62
com
ϕ
definido pela equação (4.4) e:
1
2
u
u
u
½
=
® ¾
¯ ¿
{
}
=
1 2 1 2
1 1 2 2
T
N
u u u u u
1
2
p
p
p
½
=
® ¾
¯ ¿
{
}
=
1 2 1 2
1 1 2 2
T
N
p p p p p
(4.6)
Surge a necessidade da utilização de tipos diferentes de elementos para
discretização do contorno, uma vez que pode ocorrer descontinuidade de carregamento
ou reações.
Figura 4-5 - Descontinuidade de carregamentos e reações
Serão definidos os elementos lineares utilizados: contínuo, descontínuo e de
transição.
j
Γ
0
u
=
p
pontos com
descontinuidade
63
O elemento contínuo é utilizado quando não há descontinuidade de
carregamento ou reações nas ligações com o elemento vizinho. Nesse elemento os nós
geométricos coincidem com os pontos fontes S.
Figura 4-6 - Elemento linear contínuo.
O elemento descontínuo é aquele em que, para um mesmo ponto, os valores
de carregamentos e/ou reações são diferentes daqueles representados no elemento
vizinho.
Figura 4-7 - Elemento linear descontínuo.
Uma maneira de se tratar esse problema é deslocar os pontos fonte S, que antes
coincidiam com os nós geométricos, para o exterior do corpo. Assim, para esses
pontos, utiliza-se a equação integral de contorno para pontos externos (equação (4.1))
elemento j
1
elemento j
+
p ou u
1
elemento j
+
elemento j
−
1
elemento j
p ou u
64
com
( ) = 0
ki
C S
. Define-se ‘nós duplos’ na interface dos elementos para garantir a
descontinuidade de carregamentos e/ou reações.
Figura 4-8 - Nó duplo em pontos com descontinuidade.
Para o elemento linear descontínuo, a aproximação da geometria é feita também
pela expressão (4.2), pois os pontos nodais estão localizados nos extremos dos
elementos. A aproximação das variáveis é dada pela expressão (4.5). Os pontos fonte
são deslocados para o exterior, perpendicularmente ao elemento e a uma distância que
depende do comprimento do elemento ao qual pertence.
Figura 4-9 - Novo posicionamento do ponto fonte.
O elemento de transição é aquele em que um de seus nós é tratado como
contínuo e o outro como nó duplo.
1
elemento j
+
elemento j
p ou u
1
elemento j
−
j
a L
j
a L
1
j
aL
+
1
j
a L
−
1
elemento j
+
elemento j
−
1
elemento j
p ou u
nós duplos
nós duplos
65
Figura 4-10 - Elemento de transição.
4.1.2.2 ELEMENTO QUADRÁTICO
Seja um trecho curvo do contorno definido por três pontos nodais, dois nas
extremidades e um na metade do seu comprimento (Figura 4-2 b). As coordenadas
dos pontos pertencentes ao elemento
j
Γ
podem ser calculadas pela equação (4.2)
com:
1 2 3
1 2 3
0 0 0
0 0 0
T
ϕ ϕ ϕ
ϕ
ϕ ϕ ϕ
ª º
=
« »
« »
¬ ¼
{
}
=
1 2 3 1 2 3
1 1 1 2 2 2
T
N
x x x x x x x
(4.7)
1
elemento j
+
elemento j
p ou u
1
elemento j
−
Nós Duplos
66
Figura 4-11 - Descrição geométrica do elemento quadrático.
As funções aproximadoras
( )
ϕ η
são polinômios quadráticos dependentes da
variável homogênea
η
, com características de funções de forma, isto é, cada uma tem
valor unitário para um nó e zero para os outros dois (Figura 4-12):
1
2
3
1
( - 1)
2
1
(1 )
2
(1 - )(1 )
ϕ η η
ϕ η η
ϕ η η
=
= +
= +
(4.8)
1
=
η
0
=
η
1
η
= −
2
x
1
x
2
2
x
3
2
x
1
2
x
2
1
x
1
1
x
3
1
x
67
Figura 4-12 - Funções de interpolação quadráticas.
Sendo o contorno curvo (Figura 4-13), a avaliação das integrais da equação (4.1)
requer o uso do Jacobiano, pois
i
ϕ
é função de
Ș
e as integrais são relativas ao
contorno
Γ
.
Figura 4-13 - Notação para contorno curvo.
Assim, o Jacobiano da transformação é dado por:
2 2
1 2
x x
d
J
d
η η η
∂ ∂
§ · § ·
Γ
= + =
¨ ¸ ¨ ¸
∂ ∂
© ¹ © ¹
(4.9)
2
2
r
dx
x
∂
∂
G
1
1
r
dx
x
∂
∂
G
1
x
2
x
Γ
r
G
68
Para o elemento linear
2
L
J
=
, onde L é o comprimento do elemento. O cálculo
do Jacobiano para esse elemento está no Anexo 2.
Portanto:
d J d
η
Γ =
O elemento quadrático, da mesma forma que o linear, pode ser contínuo,
descontínuo ou de transição. Para o elemento descontínuo ou de transição faz-se o
mesmo procedimento adotado para o linear, com o ponto fonte deslocado para o
exterior e a utilização de nós duplos na solução dos problemas de descontinuidade.
4.1.3 EQUAÇÕES INTEGRAIS NA FORMA MATRICIAL
A equação integral de contorno (4.1) com a eliminação dos termos inerciais
torna-se:
* *
( ) ( ) ( , ) ( ) ( ) - ( , ) ( ) ( )
ki i ki i ki i
C S u S u S Q p Q d Q p S Q u Q d Q
Γ Γ
= Γ Γ
³ ³
(4.10)
A etapa seguinte baseia-se na utilização do MEC como ferramenta numérica
propriamente dita. Isto consiste na discretização da equação integral e na formação de
um sistema matricial preparado para sua posterior solução computacional.
Tendo discretizado o contorno
Γ
em elementos
j
Γ
definidos por seus nós,
simples ou duplos, que permitem uma aproximação das variáveis parametrizadas por
seus valores nodais, a equação integral (4.10) será substituída por um somatório de
todas as integrais desenvolvidas em cada elemento
j
Γ
, chegando-se a:
69
* *
1 1
( ) ( )
NE NE
N N
T T
j j
ki i
j j
j j
C u p d u u d p
ϕ η ϕ η
= =
Γ Γ
ª º ª º
« » « »
+ Γ = Γ
« » « »
« » « »
¬ ¼ ¬ ¼
¦ ¦
³ ³
(4.11)
sendo NE o número de elementos e
N
j
u
ou
N
j
p
os vetores
N
u
ou
N
p
do j-ésimo
elemento.
A equação (4.9) permite a transformação da integral para a coordenada
adimensional
η
, através do Jacobiano. Desse modo a equação (4.11) fica:
1
*
1
-1
1
*
1
-1
( )
= ( )
NE
N
T
j
ki i
j
NE
N
T
j
j
C u p J d u
u J d p
ϕ η η
ϕ η η
+
=
+
=
ª º
« »
+ =
« »
¬ ¼
ª º
« »
« »
¬ ¼
¦
³
¦
³
(4.12)
Desenvolvendo-se as integrais da equação (4.12) pode-se escrever para cada
ponto S:
1
2
1 2
ˆ ˆ ˆ ˆ
( ) ( ) ...
:
S S SS SN
S
N
U
U
C S U S h h h h
U
U
½
° °
° °
° °
ª º
+ =
® ¾
¬ ¼
° °
° °
° °
¯ ¿
70
1
2
S1 S2 SS SN
S
N
P
P
g g ... g g
:
P
P
½
° °
° °
° °
ª º
=
® ¾
¬ ¼
° °
° °
° °
¯ ¿
(4.13)
onde N é o número total de nós. As submatrizes
SK SK
h e g
contêm os coeficientes de
influência no ponto nodal S do nó K, pertencente ao elemento j. Sua dimensão é 2x2 no
caso bidimensional e seu desenvolvimento será detalhado no item 4.1.4.
Aplicando a equação (4.13) aos N nós do contorno, obtém-se um sistema de 2N
equações lineares, contendo 2N incógnitas, entre elas deslocamentos e forças de
superfície:
11 12 1
1
11 12 1 1
21 21 2
2
21 22 2 2
1 2
1 2
..
..
..
..
:: : : : :
: : : :
..
..
(
N
N
N
N
NN N NN N
N N NN
g g g
U
h h h P
g g g
Uh h h P
C
U
h h h P
g g g
ª º
§ ·
½ª º ½
« »
¨ ¸
° °
« » ° °
« »
¨ ¸
° °
« » ° °
+ = « »
® ¾ ® ¾
¨ ¸
« »
« »
° ° ° °
¨ ¸
« »
« »
° ° ° °
¨ ¸
« »
« »
¬ ¼ ¯ ¿
¯ ¿
© ¹
¬ ¼
ˆ
)
C H U G P
+ =
(4.14)
onde
U e P
contém respectivamente os deslocamentos e as forças de superfície dos
pontos nodais,
ˆ
ˆ
e
H G
são as matrizes de coeficientes de influência da estrutura e
C
é
a matriz cujos elementos dependem da geometria, dada pela expressão (3.37).
Incorporando a matriz
ˆ
em
C H
, resultando
H
, o sistema de equações pode ser
escrito como:
71
H U G P
=
(4.15)
Num problema com N nós, o sistema (4.15) representa 2N equações lineares
algébricas, no qual
1
N
deslocamentos e
2
N
forças são conhecidos. O número total de
deslocamentos (
1
N
) e forças (
2
N
) prescritas deve ser igual a N. Com essas condições
de contorno, o sistema é agora composto de N equações correspondentes aos N
pontos nodais definidos no contorno.
Reorganizando as matrizes
e
H G
, para obter uma matriz
A
, contendo os
coeficientes referentes às incógnitas ao problema, e um vetor
B
, contendo a soma das
contribuições de
H e G
referentes aos valores prescritos, obtém-se um sistema de
equações:
A X B
=
(4.16)
onde
X
é o vetor das incógnitas.
Esse sistema pode ser resolvido através das técnicas usuais do cálculo numérico,
obtendo-se assim os deslocamentos e forças de superfícies desconhecidos.
A partir dos valores do contorno, podem-se calcular os deslocamentos e as
tensões em qualquer ponto do domínio. Para isto, utiliza-se a equação (3.23), escrita na
forma matricial:
1
*
1
-1
1
*
1
-1
( ) ( , ) ( ) -
- ( , ) ( )
NE
N
j
j
NE
N
j
j
u s u s Q J d p
p s Q J d u
ϕ η η
ϕ η η
+
=
+
=
½
° °
=
® ¾
° °
¯ ¿
½
° °
® ¾
° °
¯ ¿
¦
³
¦
³
(4.17)
72
sendo s e e
N N
j j
p u
∈ Ω
valores já determinados.
A equação (4.17) pode ser escrita como:
'
U H U G P
= − +
(4.18)
De maneira análoga, a equação integral (3.44) para cálculo de tensões será
substituída por:
1
1
1
1
1
1
( ) ( , ) ( ) -
- ( , ) ( )
NE
N
j
j
NE
N
j
j
s D s Q J d p
S s Q J d u
σ ϕ η η
ϕ η η
+
=
−
+
=
−
½
° °
=
® ¾
° °
¯ ¿
½
° °
® ¾
° °
¯ ¿
¦
³
¦
³
(4.19)
ou:
' '
H U G P
σ
= − +
(4.20)
sendo
' '
e
H G
as matrizes que se obtém ao realizar as integrais da equação (4.19).
4.1.4 MATRIZ DOS COEFICIENTES DE INFLUÊNCIA DO ELEMENTO
As submatrizes
~
~
e
SK SK
h g
obtidas a partir da equação (4.12) dependem somente
da solução fundamental utilizada e da geometria do contorno, pois:
73
1
*
-1
1
*
-1
( ) ( , ) ( , )
( ) ( , ) ( , )
m m
ki ki
m m
ki ki
h S p S Q Q J d
g S u S Q Q J d
ϕ η η
ϕ η η
+
+
=
=
³
³
(4.21)
onde:
m = indica a função aproximadora utilizada;
i = direção da carga unitária no ponto fonte;
k = direção da força(ou deslocamento) sobre o elemento em estudo;
Ș
= coordenada adimensional sobre o contorno;
S = ponto fonte, onde a carga unitária é aplicada, para o qual se escreve a
equação. Pode estar no contorno (S), no domínio(s) ou externo (f);
Q = ponto do elemento integrado.
O cálculo das integrais da expressão (4.21) é realizado levando-se em conta a
coincidência ou não do ponto fonte com um dos nós do elemento do contorno.
Quando o ponto fonte não está contido no elemento considerado, utiliza-se
integração numérica de Gauss, dada por:
( )
1
1
1
( )
n
i i
i
f d w f
η η η
+
=
−
=
¦
³
(4.22)
onde:
74
i
w
= fator peso associado a cada ponto;
n = total de pontos de integração;
i
η
= coordenada do i-ésimo ponto de integração.
Para que se obtenham resultados suficientemente precisos, deve-se escolher o
número de pontos de integração, levando-se em conta:
a) a função a ser integrada;
b) o comprimento do elemento;
c) a distância do ponto fonte ao elemento;
d) o ângulo formado pelo vetor posição.
No Anexo 3 estão listados os pontos de integração, as coordenadas
adimensionais desses pontos e o fator peso.
Quando o ponto fonte pertence ao elemento, a integração é realizada
analiticamente devido à singularidade apresentada quando r = 0. As integrais analíticas
serão aplicadas para os elementos contínuos.
Para o elemento linear as funções aproximadoras são:
1
-
1 0
1 - > 0
r a
L
r a
L
ϕ
§ ·
+ Γ ≤
¨ ¸
°
© ¹
°
°
=
®
°
+
§ ·
°
Γ
¨ ¸
°
© ¹
¯
(4.23)
75
2
-
0
> 0
a r
L
a r
L
ϕ
Γ ≤
°
°
=
®
°
+
°
Γ
¯
(4.24)
Figura 4-14- Elemento a ser integrado.
4.1.4.1 SUBMATRIZ G
Os termos da submatriz
g
são dados pela expressão:
*
m m
ki ki
g u d
ϕ
= Γ
³
(4.25)
onde m indica a função aproximadora utilizada.
A expressão para o cálculo da função de Interpolação
1
ϕ
é dada pela integral:
1
2
r
r
a
L a
−
L
Γ
ponto de integração
76
1 *
1 , , 1
, ,
0
-
, ,
0
-1
(3 - 4 ) ln -
8 (1- )
1 ( - )
(3 - 4 ) ln - 1
8 (1- )
-1 ( )
(3 4 ) ln - 1-
8 (1- )
ki ki k i
a
ki k i
L a
ki k i
u d r r r d
G
r a
r r r dr
G L
r a
r r r dr
G L
ϕ ν δ ϕ
π ν
ν δ
π ν
ν δ
π ν
½
ª º
Γ = Γ
® ¾
¬ ¼
¯ ¿
½
−
§ ·
ª º
= + +
® ¾
¨ ¸
¬ ¼
© ¹
¯ ¿
½
+
§ ·
ª º
+ −
® ¾
¨ ¸
¬ ¼
© ¹
¯ ¿
³ ³
³
³
(4.26)
Resolvendo as integrais e rearranjando os termos chega-se a:
[ ]
1 *
2 2
, ,
-1 ( - )
( - ) (3 4 ) ln ln( - ) -
8 (1- ) L 2
( - )
-L (2ln -1) (3 - 4 )
4 4 2
ki ki
ki k i
L a
u d L a a a L a
G
L a a L
a r r
ϕ δ ν
π ν
δ ν
ª
ª
Γ = − +
®
«
«
¬
¬
¯
ª º ½
º º
°
+ + −
« »
¾
» »
« »
¼ ¼
°
¬ ¼ ¿
³
(4.27)
Para o ponto fonte aplicado no nó inicial (a = 0):
1 *
, ,
(3 4 ) (1,5 - ln L) -
16 (1- )
ki ki k i
L
u d r r
G
ϕ ν δ
π ν
ª º
Γ = −
¬ ¼
³
(4.28)
Para o ponto fonte aplicado no nó final (a = L):
1 *
, ,
(3 4 ) (0,5 - ln L)
16 (1- )
ki ki k i
L
u d r r
G
ϕ ν δ
π ν
ª º
Γ = − +
¬ ¼
³
(4.29)
A expressão para o cálculo da função de Interpolação
2
ϕ
é dada pela integral:
77
2 *
, , 2
, ,
0
-
, ,
0
-1
(3 4 ) ln -
8 (1- )
-1 -
(3 4 ) ln -
8 (1- )
-1
(3 4 ) ln -
8 (1- )
ki ki k i
a
ki k i
L a
ki k i
u d r r r d
G
a r
r r r dr
G L
r r r
G
ϕ ν δ ϕ
π ν
ν δ
π ν
ν δ
π ν
½
ª º
Γ = − Γ =
® ¾
¬ ¼
¯ ¿
½
§ ·
ª º
= − +
® ¾
¨ ¸
¬ ¼
© ¹
¯ ¿
½
ª º
+ −
® ¾
¬ ¼
¯ ¿
³ ³
³
³
a r
dr
L
+
§ ·
¨ ¸
© ¹
(4.30)
Calculando as integrais e rearranjando os termos chega-se a:
2 2 2
2 *
2
, ,
-1
(3 - 4 ) ln( - ) ln -
8 (1- ) L 2 2
-L -
2 4 2
ki ki
k i
L a a
u d L a a
G
a L L
r r
ϕ δ ν
π ν
ª
−
Γ = +
«
®
«
¯
¬
º ½
§ ·
+
¾
¨ ¸
»
© ¹
¼ ¿
³
(4.31)
Para o ponto fonte aplicado no nó inicial (a = 0):
2 *
, ,
(0,5 - ln L)(3 - 4 )
16 (1- )
ki ki k i
L
u d r r
G
ϕ ν δ
π ν
ª º
Γ = +
¬ ¼
³
(4.32)
Para o ponto fonte aplicado no nó final (a = L):
2 *
, ,
(3 - 4 ) (1,5 - ln L)
16 (1- )
ki ki k i
L
u d r r
G
ϕ ν δ
π ν
ª º
Γ = +
¬ ¼
³
(4.33)
4.1.4.2 SUBMATRIZ H
Os termos da submatriz
h
são dados por:
78
*
m m
ki ki
h p d
ϕ
= Γ
³
(4.34)
O valor de
*
ki
p
vem da expressão (3.26):
( )
*
, , , ,
-1
1- 2 2 - (1- 2 )( - )
4 (1- )
ki ki k i k i i k
r
p r r r n r n
r n
ν δ ν
π ν
∂
½
ª º
= +
® ¾
¬ ¼
∂
¯ ¿
Nesse caso, em que o ponto fonte pertence ao elemento a ser integrado, r é
ortogonal a n, logo
r
n
∂
∂
= 0 e:
{ }
*
, ,
1
(1 2 )( -
4 (1- )
ki k i i k
p r n r n
r
ν
π ν
= − (4.35)
A expressão para o cálculo da função de Interpolação
1
ϕ
é dada pela integral:
( )
( )
{ }
( )
( )
( )
( )
1 *
, , 1
, ,
0
-
0
1
1- 2 -
4 (1- )
-
1 1
1 2 - 1
4 (1- )
1
1-
ki k i i k
a
k i i k
L a
p d r n r n d
r
r a
r n r n dr
r L
r a
dr
r L
ϕ ν ϕ
π ν
ν
π ν
Γ = Γ
§ ·
½
= − + +
¨ ¸
® ¾ ®
¯ ¿ ¯
© ¹
§ ·
+
½
+
¨ ¸
¾
¿
© ¹
³ ³
³
³
(4.36)
Para o ponto fonte aplicado no nó final ( a = L ):
1 *
1- 2
-
4 (1- )
ki
p d K
ν
ϕ
π ν
§ ·
Γ =
¨ ¸
© ¹
³
(4.37)
Onde:
79
0 se
1 se >
-1 se <
K i
K K i
K i
=
°
°
°
=
®
°
°
°
¯
(4.38)
A expressão para o cálculo da função de Interpolação
2
ϕ
é dada pela integral:
( )
( )
( )
( )
( )
( )
2 *
, 2
, ,
0
-
0
1
1 2 -
4 1-
1 1 -
1 2 - -
4 1-
1
-
ki k i i k
a
k i i k
L a
p d r n r n d
r
a r
r n r n dr
r L
a r
dr
r L
ϕ ν ϕ
π ν
ν
π ν
½
Γ = − Γ
® ¾
¯ ¿
½
§ ·
= −
® ¾ ®
¨ ¸
© ¹
¯ ¿ ¯
½
+
§ ·
¾
¨ ¸
© ¹
¿
³ ³
³
³
(4.39)
Para o ponto fonte aplicado no nó inicial (a = 0):
2 *
1- 2
4 (1 )
ki
p d K
ν
ϕ
π ν
§ ·
Γ =
¨ ¸
−
© ¹
³
(4.40)
Nos casos em que o ponto fonte é aplicado no nó inicial do elemento com a
função de interpolação
1
ϕ
e o ponto fonte é aplicado no nó final com a função de
interpolação
2
ϕ
, a integral deve ser calculada sob o conceito do valor principal de
Cauchy. A integral sobre os dois elementos que contém o nó deve ser feita
simultaneamente para eliminar a singularidade. Não é possível expressar o valor de
m
ki
h
individualmente, mas o valor total de cada componente da submatriz
ki
H
correspondente ao nó em estudo pode ser calculado como:
80
2
1
1- 2
ln
4 (1 )
ki
L
H K
L
ν
π ν
§ ·
§ ·
=
¨ ¸
¨ ¸
−
© ¹
© ¹
(4.41)
Figura 4-15- Integrais simultâneas para
~
h
quando o ponto fonte está na interface.
Nas expressões(4.40) e (4.41) o valor de K é dado por (4.38).
É possível verificar a precisão das integrais realizadas sobre elementos de
contorno através de uma propriedade da matriz (equação (4.15))
Impondo a condição de que translações de corpo rígido correspondem a forças
de superfície nulas, aplicam-se duas translações independentes
1 2
e
k k k k
u u
δ δ
= =
e chega-se a seguinte relação:
1
0
N
nm mn
m
h u
=
=
¦
( n = 1,........, N) (4.42)
sendo
nm
h
submatrizes (2 X 2) de
H
e
m
u
as translações escritas como:
m
u I
=
(4.43)
onde
I
representa a matriz identidade de ordem 2.
1
ϕ
2
ϕ
1
J
=
2
J
=
1
1
2
2
81
Assim:
11 12 1
21 22 2
1 2
...
0
:
:
...
:
:
: : : :
0
...
N
N
N N NN
h h h
I
h h h
I
h h h
ª º
½
½
« »
° °
° °
« »
° ° ° °
=
® ¾ ® ¾
« »
° ° ° °
« »
° ° ° °
« »
¯ ¿
¯ ¿
¬ ¼
(4.44)
Onde:
1
0
N
nm
ki
m
h
=
=
¦
( n , i , k = 1,......, N) (4.45)
4.1.5 CONDENSAÇÃO DE DESLOCAMENTOS EM NÓS DUPLOS
Quando se utiliza os nós duplos para melhor representar as descontinuidades de
força ou deslocamento num ponto do contorno, criam-se linhas e colunas iguais nas
matrizes
H e G
, uma vez que as ordenadas são as mesmas para os dois nós.
Na formulação estática, quando são impostas as condições de contorno, esse
problema desaparece, pois sempre um dos nós tem deslocamento ou força conhecidos.
Na formulação dinâmica, através da reciprocidade dual (3.2.2) há a necessidade
de se inverter a matriz
F
, que relaciona aceleração
u
com o vetor
α
, dependente do
tempo. Com os nós duplos, as linhas e colunas de
F
, correspondentes a esses nós,
são linearmente dependentes e a matriz não pode ser invertida.
82
Uma solução para esse problema, proposta nesse trabalho, é um arranjo de
matrizes, eliminando linhas e colunas de um dos nós duplos. Essa proposta foi testada
nos problemas estáticos, com soluções conhecidas, para comprovar sua validade. O
mesmo raciocínio, de maneira ampla, será utilizado na resolução dos problemas
dinâmicos.
Na formulação estática, o sistema de equações tem a forma (4.15):
H U G P
=
onde
e
H G
são matrizes quadradas de ordem 2N , sendo N o número de nós do
problema.
Faz-se uma distinção entre os dois tipos de condição de contorno:
1
Γ será a
região do contorno em que os deslocamentos são prescritos, e
2
Γ
onde as forças são
prescritas.
Impondo-se as condições de contorno do problema, chega-se o sistema:
ˆ
ˆ
H U G P
=
(4.46)
O novo vetor
U
contém apenas forças e deslocamentos incógnitos e o vetor
P
contém os valores conhecidos.
Os vetores
U
(
P
) podem ser particionados definindo-se:
(
)
M M
U P
= deslocamentos (forças) em nós duplos considerados “master”;
(
)
D D
U P
= deslocamentos (forças) em nós duplos considerados “dependentes”;
83
(
)
I I
U P
= deslocamentos (forças) em nós não duplos chamados
“independentes”.
A equação (4.46) é rescrita na forma:
ˆ ˆ ˆ
ˆ ˆ ˆ
ˆ ˆ ˆ
ˆ ˆ ˆ
ˆ ˆ ˆ ˆ ˆ ˆ
DD DM DI
DD DM DI
D D
MD MM MI M MD MM MI M
I I
ID IM II
ID IM II
G G G
H H H
U P
H H H U G G G P
U P
H H H
G G G
ª º
ª º
½ ½
« »
« »
° ° ° °
=
« »
« »
® ¾ ® ¾
« »
« »
° ° ° °
¯ ¿ ¯ ¿
« »
« »
¬ ¼
¬ ¼
(4.47)
com
ˆ
ˆ
e
ij ij
H G
s submatrizes da partição.
Sabe-se que:
M D
U U
=
(4.48)
Da imposição (4.48) escreve-se:
0
0
0
D
M M
M
I I
I
U I
U U
U I A
U U
U I
½ ª º
½ ½
° °
« »
= =
® ¾ ® ¾ ® ¾
« »
¯ ¿ ¯ ¿
° °
« »
¯ ¿ ¬ ¼
(4.49)
Substituindo a equação (4.49) em (4.47) e pré-multiplicando o sistema resultante
por
T
A
:
84
ª º
ª º
« »
½
ª º
« »
=
« »
® ¾
« »
« »
¬ ¼
« »
¯ ¿
« »
¬ ¼
« »
¬ ¼
ª º
« »
ª º
« »
« »
« »
¬ ¼
« »
¬ ¼
ˆ ˆ ˆ
0
0
ˆ ˆ ˆ
0
0 0
ˆ ˆ ˆ
0
ˆ ˆ ˆ
0
ˆ ˆ ˆ
=
0 0
ˆ ˆ ˆ
DD DM DI
M
MD MM MI
I
ID IM II
DD DM DI
D
MD MM MI M
ID IM II
H H H
I
U
I I
H H H I
U
I
I
H H H
G G G
P
I I
G G G P
I
G G G
½
° °
® ¾
° °
¯ ¿
I
P
(4.50)
Efetuando o produto das matrizes, chega-se a:
ˆ ˆ ˆ ˆ ˆ ˆ
ˆ ˆ ˆ
ˆ ˆ ˆ ˆ ˆ ˆ
ˆ ˆ ˆ
DD MD DM MM DI MI
M
I
ID IM II
DD MD DM MM DI MI D
M
I
ID IM II
H H H H H H
U
U
H H H
G G G G G G P
P
P
G G G
ª º
+ + + +
« » ½
=
® ¾
« »
¯ ¿
« »
+
¬ ¼
ª º
+ + +
½
« »
° °
=
® ¾
« »
° °
« »
¯ ¿
¬ ¼
(4.51)
Na equação (4.51) foram eliminados os deslocamentos incógnitos nos nós
duplos dependentes, no lado esquerdo da equação.
A comparação dos resultados obtidos através da metodologia de condensação
de matrizes com aqueles que utilizam a maneira tradicional estão no Quadro (6.6).
4.2 FORMULAÇÃO DINÂMICA
4.2.1 INTRODUÇÃO
A equação integral de contorno que governa o comportamento dinâmico dos
corpos (5.2) é escrita como:
85
* *
( ) ( ) ( , ) ( , ) ( ) - ( , ) ( , ) ( )
ki i ki i ki i
C S u S u S Q p Q t d Q p S Q u Q t d Q
Γ Γ
= Γ Γ +
³ ³
*
( , ) ( , ) ( )
ki i
u S q u q t d q
ρ
Ω
+ Ω
³
(4.52)
No caso de vibrações livres, analisadas no item (5.6), a equação integral tem a
forma:
* *
( ) ( ) ( , ) ( , ) ( ) - ( , ) ( , ) ( ) -
ki i ki i ki i
C S u S u S Q p Q t d Q p S Q u Q t d Q
Γ Γ
= Γ Γ
³ ³
2 *
( , ) ( )
ki
u S q d q
ω ρ
Ω
− Ω
³
(4.53)
Do mesmo modo que para a formulação estática, será realizada a transformação
das equações (4.52) e (4.53) num conjunto de equações algébricas que possam ser
resolvidas pelo Método dos Elementos de Contorno.
O contorno
Γ
do corpo analisado vai ser discretizado em elementos
j
Γ
, exatos
ou aproximados, onde se interpolam as variáveis físicas do problema. A discretização
geométrica (4.2) é válida também para essa análise.
As integrais de domínio que aparecem no último termo da direita nas equações
(4.52) e (4.53) foram resolvidas de duas maneiras, descritas a seguir.
Inicialmente as integrais de domínio foram resolvidas pela Técnica da
Reciprocidade Dual, idealizada por Nardini e Brebbia em 1982. Ela emprega
diretamente a solução fundamental da estática, associada a um procedimento de
substituição de variáveis, que visa eliminar as integrais de domínio existentes mediante
a aplicação do Teorema da Reciprocidade de Betti uma outra vez. Na análise de
86
problemas de vibração livre e transiente será utilizada a condensação das matrizes,
proposta no item anterior, na formação do sistema de equações.
Em seguida, apresenta-se na seção (4.2.3) o método da integração direta, que
transforma os termos da integral de domínio em algébricos, com o emprego da técnica
da discretização em células.
4.2.2 FORMULAÇÃO COM RECIPROCIDADE DUAL
O Método da Reciprocidade Dual consiste num modo geral de se construir
soluções particulares que podem ser usadas para resolver problemas não lineares ou
dependentes do tempo, como também para representar qualquer distribuição de carga
em pontos internos, ou seja, pontos do domínio. É uma formulação que, mediante o uso
de funções auxiliares, permite a transformação das integrais de domínio em integrais de
contorno, de acordo com os procedimentos clássicos do Método dos Elementos de
Contorno.
Para a formulação da equação integral com a solução fundamental estática,
podem-se restringir as influências externas somente para o contorno e, nesse caso, a
equação (4.4) fica:
,
- 0
ij j i
u
σ ρ
=
(4.54)
O Teorema da Reciprocidade de Maxwell-Betti, que relaciona dois modos de
deslocamentos independentes, pode ser escrito como:
(
)
(
)
* * * *
, ,
- -
ki ij j i kij j ki i i ki
u u d u p u p d
σ σ
Ω Γ
Ω = Γ
³ ³
(4.55)
Na equação (4.55)
i
u
é o primeiro campo de deslocamentos, com as tensões
ij
σ
e as forças
i
p
, e
*
i
u
o segundo campo de deslocamentos com tensões
*
ij
σ
e forças
*
i
p
.
87
Esse campo de deslocamentos é virtual e pode ser associado à solução da equação
estática com carga pontual, ou seja, à solução de:
*
,
( - ) 0
kij j ki
x
σ δ ξ
+ ∆ =
(4.56)
onde:
(
)
( )
0 para
1
x x
x d
ξ ξ
ξ
Ω
∆ − = ≠
∆ − Ω =
³
(4.57)
A solução da equação (4.57) é
*
ki
u
e a força correspondente
*
ki
p
. Com as
equações (4.56), (4.57), (4.55), obtém-se:
* * *
,
- -
ij j ki ki i ki i ki i
u d c u u p d p u d
σ
Ω Γ Γ
Ω = + Γ Γ
³ ³ ³
(4.58)
Substituindo (4.54) em (4.58), com
ρ
constante:
* * *
=
ki i ki i ki i i ki
c u p u d u p d u u d
ρ
Γ Γ Ω
+ Γ Γ + Ω
³ ³ ³
(4.59)
Para a solução dessa equação integral, é possível utilizar a solução fundamental
simples evitando assim a solução fundamental dependente do tempo utilizada no
Método dos Elementos de Contorno convencional. Entretanto, a presença do integral no
volume indica que uma discretização adicional no domínio é necessária. Nardini e
Brebbia [1985] transformaram essa integral de volume numa integral no contorno,
utilizando a solução fundamental da elastostática e chegando ao Método da
Reciprocidade Dual (DRM).
88
De maneira análoga ao caso mais simples, que é a equação de Poisson, o
Método da Reciprocidade Dual propõe o uso de uma série de soluções particulares
i
u
para a função deslocamento, ao invés de uma única solução
u
. O número de soluções
i
u
utilizados é igual ao número de pontos nodais total do problema analisado (m), ou
seja, a soma dos nós do contorno com os pontos internos.
O deslocamento
( , )
i
u t
ξ
é expresso como a soma de m funções
( )
j
f
ξ
, que
dependem da geometria, multiplicadas pelas funções incógnitas
( )
j
i
t
α
, que dependem
do tempo:
( , ) ( ) ( )
j j
i
i
u t t f
ξ α ξ
=
(4.60)
com somatório de j = 1 até m, sendo o ponto
ξ
um ponto qualquer do domínio.
O índice j da função de interpolação
j
f
mostra sua associação com o ponto j, o
qual é também um dos m pontos de colocação.
A aceleração fica:
( , ) ( ) ( )
j j
i
i
u t t f
ξ α ξ
=
(4.61)
A expansão (4.61) é válida sobre todo o domínio, como no caso de um super-
elemento.
A integral de domínio do lado esquerdo da equação (4.59) pode ser escrita
como:
* * *
j j j j
i ki ki li ki
i l
u u d f u d f u d
α α δ
Ω Ω Ω
Ω = Ω = Ω
³ ³ ³
(4.62)
89
Essa integral de domínio pode ser transformada numa equação integral de
contorno equivalente. O primeiro passo é procurar a solução do problema estático com
contorno infinito, ou seja, a solução de:
lim,
0
j
m li
f
σ δ
+ =
(4.63)
O índice j representa a base da função carregamento, os índices l e i são
respectivamente, a direção em que o carregamento
( )
j
f
ξ
atua e a direção do
deslocamento ou força de superfície resultante.
A solução da equação (4.63) em termos de deslocamentos será chamada de
j
li
ψ
e as suas correspondentes forças de superfície serão
Ș
j
li
. Essas duas funções
dependem exclusivamente das funções
j
f
, que são arbitrárias.
Pode-se reduzir a equação de domínio (4.62) a uma equação de contorno,
através da transformação (4.58):
* * * *
lim,
- -
j j j j j
li ki ki ki ki ki
m li li li
f u d u d c u d p d
δ σ ψ η ψ
Ω Ω Γ Γ
Ω = Ω = Γ + Γ
³ ³ ³ ³
(4.64)
Substituindo esse resultado em (4.62) e rearranjando (4.59), tem-se:
* *
* *
-
- 0
ki i ki i ki i
j j j j
ki ki ki
li li li l
c u p u d u p d
c p d u d
ρ ψ ψ η α
Γ Γ
Γ Γ
+ Γ Γ +
§ ·
¨ ¸
+ + Γ Γ =
¨ ¸
© ¹
³ ³
³ ³
(4.65)
Essa é a equação de contorno cuja solução numérica vai ser calculada, através
do Método dos Elementos de Contorno.
90
No procedimento descrito acima foi utilizada a aplicação do Teorema da
Reciprocidade de Betti uma segunda vez, na transformação da integral de domínio em
uma soma de integrais de contorno. Daí o nome de Reciprocidade Dual.
Aplicando-se uma força unitária em cada nó, a qual produz um campo de
deslocamento
*
ki
u
, a equação (4.65) se torna um sistema de equações da forma:
* *
* *
( ) ( ) ( , ) ( ) - ( , ) ( )
( ) ( ) ( , ) ( ) - ( , ) ( ) 0
ki i ki i ki i
j j j j
ki ki ki
li li li l
c A u A p A B u B d u A B p B d
c A A p A B B d u A B B d
ρ ψ ψ η α
Γ Γ
Γ Γ
+ Γ Γ +
§ ·
¨ ¸
+ + Γ Γ =
¨ ¸
© ¹
³ ³
³ ³
(4.66)
onde A é o n-ésimo ponto nodal e B é o ponto de integração no contorno.
Os deslocamentos u e as forças de superfície p dos pontos do contorno são
expressos em termos dos deslocamentos e forças de superfície nodais, relacionados
pelas funções de interpolação
n
ϕ
,de acordo com a equação (4.6):
T n
i
T n
u u
p u
ϕ
ϕ
=
=
Para uma classe de funções
j
f
escolhidas, as funções
Ș
e
ȥ
são conhecidas no
contorno e as integrais que as contém podem ser calculadas como na equação (4.66).
Porém, se for utilizado na aproximação o mesmo conjunto de funções de interpolação
ϕ
, os coeficientes serão os mesmos, e o trabalho computacional torna-se menor.
Assim:
91
j n nj
li li
j n nj
li li
ψ φ ψ
η φ η
=
=
(4.67)
onde n é a numeração local dos nós de cada elemento, e
nj nj
li li
e
ψ η
são os valores
nos nós das funções
j j
li li
e
ψ η
correspondentes a uma base genérica.
Substituindo-se as equações (4.6) e (4.67) em (4.66), resulta uma expressão
para a equação integral que pode ser escrita na forma de um somatório de integrais
sobre todos os elementos
e
Γ
que compõem o contorno
Γ
. Assim:
* *
1
* li
j
1
* li
j
( ) ( , ) ( ) - ( , ) ( ) ( )
( ) ( , ) ( ) ( ) -
- ( , ) ( ) ( )
NE
T T
ki ki i ki i
e
e e
NE
j T
ki ki
li
e
e
T
ki
e
c A p A B d B U u A B B d B P
c A p A B B d B
u A B B d B
ϕ ϕ
ρ ψ ϕ ψ
ϕ η
=
Γ Γ
=
Γ
Γ
ª
º
«
+ Γ Γ +
»
«
¼
¬
ª
+ + Γ
®
«
¯ ¬
º ½
Γ
»
¼
¦
³ ³
¦
³
³
0
j
l
α
=
¾
¿
(4.68)
Seguindo o procedimento do item (4.1.3) chega-se ao sistema de equações
lineares, expresso em termos das matrizes globais:
- ( - ) 0
H U G P H G
ρ ψ η α
+ =
(4.69)
A matriz global
H
é formada pelas matrizes
n
h
dos elementos e das matrizes
diagonais
n
c
. A matriz global
G
vem das matrizes
n
g
dos elementos. As matrizes
ψ
e
η
contêm os valores das funções nos pontos nodais.
92
A definição de pontos internos faz com que as soluções obtidas sejam mais
precisas. Como esses nós não fazem parte dos elementos, suas coordenadas devem
ser dadas de maneira independente. As matrizes da equação (4.69) são particionadas
em submatrizes NC (nós do contorno) e NI (pontos internos). Sendo
I
a matriz
identidade, escreve-se essa equação como:
-
0 0
. - .
0 0
0
. -
NC NC NC NC
NI NI NI
NC NI
NC
NI NC NI
H U G P
H I U G
H
H I
ψ
ψ
+
+
ª º ½ ª º ½
° ° ° °
« » « »
… … +
® ¾ ® ¾
« » « »
° ° ° °
« » « »
¬ ¼ ¯ ¿ ¬ ¼ ¯ ¿
½
ª º
° °
° °
« »
+ …
® ¾
« »
° °
« »
¬ ¼
° °
¯ ¿
# #
" " " "
# #
#
" "
#
-
0
. 0
0 0
NC NI
NC NC
NI NI
G
G
η
α
α
+
ª º
½
ª º ½
« »
° °
° ° ° °
« »
« »
… + =
® ¾ ® ¾
« »
« »
° ° ° °
« »
« »
¬ ¼ ¯ ¿
° °
« »
¯ ¿
¬ ¼
#
" "
#
(4.70)
Como as condições de contorno são expressas em termos dos deslocamentos,
toma-se
U
como base. A transformação de
U
para
α
está na equação (4.60).
Aplicando-se essa equação para cada ponto nodal, obtém-se:
U F
α
=
(4.71)
onde cada coluna da matriz
F
consiste de um vetor
j
f
, que contém os valores das
funções
j
f
nos m pontos nodais.
A escolha de um conjunto de funções linearmente independente, de número
igual ao número de pontos nodais, garante que a matriz
F
possui inversa. Logo:
93
1
1
F U
F U
α
α
−
−
=
=
(4.72)
Substituindo
α
na equação (4.71):
H U M U G P
+ =
(4.73)
com a matriz de massa dada por:
-1
( - )
M H G F
ρ ψ η
=
(4.74)
A escolha de uma classe de funções
j
f
, usadas para aproximação da integral de
domínio vai influenciar os resultados. Para o DRM a função aproximadora tem sido uma
limitação. A função tem que representar uma série independente, quando referida a
vários polos deve levar às soluções particulares da equação diferencial que rege o
equilíbrio de um corpo elástico e ser diferenciável até segunda ordem, para que a
representação integral seja válida. Alguns estudos foram feitos por Loeffler (1988),
Coda (1990), Coda e Venturini (1990) e Golberg et al (1999) com relação à adequação
de alguns tipos de função a certos problemas específicos.
Desde os primeiros trabalhos publicados sobre o assunto tem sido utilizada uma
família de funções escritas em termo da distância Euclidiana entre dois pontos (r).
Funções desse tipo escritas para polos diferentes serão sempre independentes. As
principais conclusões sobre os resultados obtidos com essas funções são que elas não
apenas levam à convergência das soluções, como também reduzem o tempo
computacional necessário na solução de problemas através da utilização do MEC e
DRM.
Neste trabalho é adotada a função:
94
( ) 1 ( , )
j
j
f r A
ξ ξ
= + (4.75)
onde
( , )
j
r A
ξ
é a distância entre o ponto
j
A
, no qual é definida a função, e o ponto
ξ
.
Com base na expressão (4.75) e na equação de Navier para o caso estático, têm-
se as funções
η
ψ
e
:
3
2
, , , ,
(1- 2 ) 10
3 -
(5 - 4 ) 30 (1- ) 3
li l i li l i
r
r r r r r
G G
ν ν
ψ δ
ν ν
ª º
§ ·
= + −
¨ ¸
« »
© ¹
¬ ¼
e
( )
( )
, , ,
2
, , , , , ,
2 (1- 2 ) 1 1 1
5 - 4 1 2 2 2
4 - 5 - - (1- 5 )
15 (1- )
li l i i l li j j
li j j l i l i j j i l
r n r n n r r
r
n r n r r r r n n r
G
ν ν
η δ
ν ν
ν δ ν
ν
+
ª º
= + + +
« »
−
¬ ¼
ª º
+ +
¬ ¼
(4.76)
O cálculo detalhado para a obtenção da expressão da matriz
ψ
está no Anexo 4.
Torna-se necessária uma análise do sinal da distância
( , )
j
r A
ξ
e de suas
derivadas.
No cálculo dos elementos das matrizes
H
e
G
é definido um ponto fonte k e a
integração é feita sobre todo o contorno. Então
ki
r
é o módulo do vetor
ki
r
G
, onde i é o
ponto que varia ao longo do contorno (Figura 4-16):
2 2
= ( - ) ( - )
ki i k i k
r x x y y+
(4.77)
95
No cálculo dos elementos das matrizes
ψ
e
η
é definido um ponto fonte i e a
integração é feita sobre todo o contorno. Então
il
r
é o módulo do vetor
il
r
G
, onde l é o
ponto que varia ao longo do contorno(Figura 4-16):
2 2
( - ) ( - )
il l i l i
r x x y y= +
(4.78)
Figura 4-16 - Vetores
ki
r
e
il
r
.
Como aparecem as derivadas de r nas equações (4.76), um cuidado especial
deve ser tomado para o uso dos sinais corretos, para não acarretar erros na solução do
problema.
4.2.3 FORMULAÇÃO COM INTEGRAÇÃO DIRETA
Uma maneira mais simples para se resolver o problema transiente bidimensional é
aproximar os termos inerciais incógnitos por polinômios interpoladores parametrizados
em células. Essas células são subdomínios que formam a superfície do corpo.
ki
r
il
r
i
k
l
(
)
Ponto Fonte
(
)
Ponto Campo
C
Ponto Variável
no ontorno
§ ·
¨ ¸
© ¹
96
Essa técnica leva à integração de todo o domínio em forma discretizada e os
resultados obtidos são bastante precisos, motivo pelo qual está sendo apresentada.
O problema a ser resolvido, idêntico ao item 4.2.1, com as influências externas
restritas ao contorno, é dado pela equação integral (4.52):
* *
( ) ( ) ( , ) ( , ) ( ) - ( , ) ( , ) ( )
ki i ki i ki i
C S u S u S Q p Q t d Q p S Q u Q t d Q
Γ Γ
= Γ Γ +
³ ³
*
( , ) ( , ) ( )
ki i
u S q u q t d q
ρ
Ω
+ Ω
³
(4.79)
A discretização do domínio por células permite transformar em termos
algébricos, os termos que aparecem na integral de domínio da equação (4.79). A
discretização consiste na divisão do domínio
Ω
em células triangulares, definidas pelos
seus nós, que são os vértices do triângulo.
Figura 4-17 - Divisão do domínio em células.
97
A posição de cada ponto da célula será dada em função da posição dos vértices,
através das coordenadas homogêneas
ξ
. Para um ponto q qualquer, pertencente à
célula I, tem-se:
1
~ ~
2
q
N
q
x
K X
x
½
=
® ¾
¯ ¿
(4.80)
Onde:
1 2 3
2
~
1 2 3
0 0 0
0 0 0
q q q
q q
q
K
ξ ξ ξ
ξ ξ ξ
ª º
=
« »
¬ ¼
{ }
1 2 3 1 2 3
1 1 1 2 2 2
~
T
N
X x x x x x x=
A equação (4.80) pela qual são calculadas as coordenadas
1 2
e
x x
do ponto q
como função de
ξ
(Figura 4-18), pode ser escrita:
1 2 3
1 1 1 2 1 3 1
1 2 3
2 1 2 2 2 3 2
q q q q
q q q q
x x x x
x x x x
ξ ξ ξ
ξ ξ ξ
= + +
= + +
(4.81)
Figura 4-18 - Coordenadas homogêneas.
1
0
ξ
=
2
0
ξ
=
3
0
ξ
=
1 2 3
( , , )
q
ξ ξ ξ
1
q
2
q
3
q
1
ξ
2
ξ
1
x
2
x
98
As variáveis
1 2 3
,
q q q
e
ξ ξ ξ
são normalmente expressas em função do sistema de
coordenadas
1 2
e
x x
, como no Método dos Elementos Finitos. Mas, como a integral de
domínio a ser calculada (equação (4.79)) contém a solução fundamental
*
ik
u
, é mais
conveniente a utilização de coordenadas polares. Isso porque
*
ik
u
(equação 3.25) é
função do raio r, que é a distância do ponto q ao ponto fonte s, e do ângulo
θ
entre a
direção
1
x
e a direção do ponto fonte (Figura 4-19).
Figura 4-19 - Coordenadas polares.
Desse modo, tem-se:
1 1
1 1
2 2
2 2
3 3
3 3
cos
2
cos
2
cos
2
q s
q s
q s
r
b a sen
A
r
b a sen
A
r
b a sen
A
ξ ξ θ θ
ξ ξ θ θ
ξ ξ θ θ
§ ·
= + +
¨ ¸
© ¹
§ ·
= + +
¨ ¸
© ¹
§ ·
= + +
¨ ¸
© ¹
(4.82)
2
ξ
1
0
ξ
=
2
0
ξ
=
3
0
ξ
=
1 2 3
( , , )
q
ξ ξ ξ
1
q
2
q
3
q
1
ξ
1
θ
3
θ
2
θ
θ
99
Com:
1 2 2 2 1 2 3 2 1
1 1 1 1 1 1
1 2 3 2 3 1 3 1 2
2 2 2 2 2 2
1 2 2 1
- - -
- - -
1
( - )
2
a x x a x x a x x
b x x b x x b x x
A b a b a
= = =
= = =
=
A aproximação discretizada dos deslocamentos é feita da mesma maneira que a
discretização geométrica, ou seja, os valores dos deslocamentos incógnitos do ponto q
serão aproximados linearmente em função dos valores nodais, usando a coordenada
homogênea
ξ
(Figura 4-20).
Figura 4-20 - Aproximação linear dos deslocamentos.
Assim:
1
~
~
2
q
N
q
u
K U
u
½
=
® ¾
¯ ¿
(4.83)
1
x
2
x
j
u
1
j
u
3
j
u
2
j
u
1
q
2
q
3
q
i
Ω
100
Onde:
1 2 3
~
1 2 3
0 0 0
0 0 0
q q q
q q q
q
K
ξ ξ ξ
ξ ξ ξ
ª º
=
« »
¬ ¼
{ }
1 2 3 1 2 3
1 1 1 2 2 2
~
T
N
U u u u u u u=
A equação (4.83) pode ser escrita como:
1 2 3
1 2 3
q q q q
j j j j
u u u u
ξ ξ ξ
= + +
(4.84)
com as coordenadas
ξ
dadas por (4.82).
De modo análogo, obtém-se as expressões para as acelerações
1 2
e
u u
.
A equação integral (4.79) será escrita como uma soma de integrais sobre os
elementos do contorno e sobre as células do domínio, após a discretização (Figura
4-21):
*
1
*
* ( )
1
( ) ( ) = ( , ) ( ) ( ) -
- ( , ) ( ) ( )
( , ) ( ) ( ) ( )
NE
T
ki i ki i
e
e
ki i
e
NC
m
ki j
m
m
C S u S p S Q Q d Q U
u S Q Q d Q P
u S q K q d q U q
φ
φ
ρ
=
Γ
Γ
=
Ω
ª
«
Γ
«
¬
Γ +
ª
º
«
+ Ω
»
«
¼
¬
¦
³
³
¦
³
(4.85)
101
Figura 4-21 - Discretização do domínio e contorno.
É preciso ressaltar que os valores dos deslocamentos internos também são
incógnitos. Escreve-se então a equação para os pontos internos, considerando o ponto
fonte no interior do corpo. Na equação (4.85) isso significa substituir S (ponto do
contorno), por s (ponto do domínio).
Adotando um procedimento análogo ao item (4.1.3), a equação (4.85) pode ser
reescrita como:
1
2
1 2
:
ˆ ˆ ˆ ˆ
... ...
:
S S SS SN
S S
S
N
U
U
C U h h h h
U
U
½
° °
° °
° °
° °
ª º
+ =
® ¾
¬ ¼
° °
° °
° °
° °
¯ ¿
102
1
2
1 2
1
2
1 2
:
... ... -
:
:
- ... ... ....
:
:
S S SS SN
S
N
S
N
S S SS SN Ss SM
s
M
P
P
g g g g
P
P
U
U
U
U
m m m m m m
U
U
ρ
½
° °
° °
° °
° °
ª º
=
® ¾
¬ ¼
° °
° °
° °
° °
¯ ¿
½
° °
° °
° °
° °
° °
° °
ª º
® ¾
¬ ¼
° °
° °
° °
° °
° °
° °
¯ ¿
(4.86)
onde
ˆ
ˆ
,
SN SN
h g
são dadas pela equação (4.21);
ˆ
SS ss
S
C h h+ =
se
S
∈ Γ
;
S
C I
=
se
s
∈ Ω
;
SN
m
= contém os coeficientes de influência de massa do nó N no ponto nodal S;
N = número de nós do contorno;
M = número de nós total (contorno + internos).
Lembrando que os termos em
h
e
g
estão relacionados apenas com
deslocamentos e forças do contorno, pode-se escrever de maneira global:
103
11 1 1( 1) 1 1
1 ( 1)
1 ( 1)
.. 0 .. 0
: : : : : : :
.. 0 .. 0
: : : : : : :
: : : : : : :
.. 0 ..
N N M
N NN N N NM N
M MM M N MM M
h h U
h h U
h h I U
+
+
+
ª º ½
« » ° °
« » ° °
« » ° °
° °
=
« »
® ¾
« »
° °
« »
° °
« »
° °
« »
° °
¬ ¼ ¯ ¿
1
11 1 1( 1)
1
1 ( 1)
1 ( 1)
11 1
.. 0 .. 0
: : : : : :
:
.. 0 : 0
-
:
: : : : : :
:
: : : : : :
.. 0 : 0
- .. .. .. ..
: .. .. .. .. :
: .. .. .. .. :
-
: .
M
N N
N NN N N NM
N
M
M MM M N MM
M
g g
P
g g
P
P
g g
m m
ρ
+
+
+
ª º
½
« »
° °
« »
° °
« »
° °
° °
« »
=
® ¾
« »
° °
« »
° °
« »
° °
« »
° °
¯ ¿
« »
¬ ¼
−
1
1
:
. .. .. .. : :
: .. .. .. .. : :
.. .. .. ..
N
M MM
M
U
U
m m U
½
ª º
° °
« »
° °
« »
° °
« »
° °
« »
® ¾
« »
° °
« »
° °
« »
° °
− −
« »
° °
¬ ¼
¯ ¿
(4.87)
onde as submatrizes
0
representam, para índices maiores que N, o não relacionamento
direto entre as variáveis do domínio e os termos do contorno. As submatrizes
I
representam a existência de um ponto fonte em s e, por conseqüência, um valor
incógnito a considerar.
A equação (4.87) pode ser escrita:
H U M U G P
+ =
(4.88)
Esse problema pode ser resolvido do mesmo modo descrito no item (4.2.1), com o
reagrupamento das matrizes.
104
Com a discretização geométrica do domínio e a aproximação dos deslocamentos
e acelerações, pode-se fazer a integração sobre as células.
O termo que representa a integração sobre uma célula é expresso por:
* ( )
1
1
2
1
* * 3
1 2 311 12 1
* * 1
1 2 3
21 22 2
2
2
3
2
( , ) ( ) ( )
0 0 0
.
0 0 0
k
ki j
k
R
M M
q q q
q q q
R
m m
u S q K q d k U
u
u
u u u
r dr d
u u u
u
u
θ
θ
ρ
ξ ξ ξ
θ
ξ ξ ξ
Ω
ª º
« »
Ω =
« »
¬ ¼
½
° °
° °
ª º
° °
ª º
ª º
° °
« »
=
® ¾
« »« »
« »
¬ ¼ ¬ ¼
° °
¬ ¼
° °
° °
° °
¯ ¿
³
³ ³
(4.89)
onde:
i
j
u
é a aceleração na direção ‘j’ do nó ‘i’;
e
m M
θ θ
são os ângulos extremos de varredura do vetor r (S, q);
e
m M
R R
são os valores extremos que o vetor r (S, q) assume.
Com a representação de
ξ
em coordenadas polares (equação (4.82)), pode-se
escrever a equação (4.89) como:
*
( ) ( )
R R
M m M M
ki kij s ki
R R
k m m m m
u K d Z K dr d Y dr d
θ θ
θ
θ θ
θ ζ θ
Ω
Ω = +
³ ³ ³ ³ ³
(4.90)
onde:
105
*
* 2
1 2 3
~
( )
1 2 3
1 2 3
~
( )
1 2 3
.
.
0 0 0
0 0 0
0 0 0
0 0 0
ki ki
ki ki
s s s
s s s
s
Z u r
Y u r
K
θ θ θ
θ θ θ
θε
ξ ξ ξ
ξ ξ ξ
ζ ζ ζ
ζ
ζ ζ ζ
=
=
ª º
=
« »
¬ ¼
ª º
=
« »
¬ ¼
( )
( )
( )
1 1
1
2 2
2
3 3
3
1
cos
2
1
cos
2
1
cos
2
b a sen
A
b a sen
A
b a sen
A
θ
θ
θ
ζ θ θ
ζ θ θ
ζ θ θ
= +
= +
= +
Os valores de
e
h h
b a
são dados por (4.82).
As primitivas em r das funções
e
ki ki
Z Y
, determinadas de forma analítica, são
dadas por:
( )
( )
2 2
, ,
3
3
, ,
1 1
3 - 4 ln -
8 (1- ) 2 2 2
1 1
3 - 4 ln -
8 (1- ) 3 3 3
ki ki ki k i
r
ki ki ki k i
r
r r
Z dr M r r r
G
r
r
Y dr N r r r
G
ν δ
π ν
ν δ
π ν
ª º
½
§ ·
= = +
® ¾
« »
¨ ¸
© ¹
¯ ¿
¬ ¼
½
ª º
§ ·
= = +
® ¾
¨ ¸
« »
© ¹
¬ ¼
¯ ¿
³
³
(4.91)
A equação de cada aresta que limita a célula I em estudo é dada por:
106
- 2
( , )
cos
s
k
ml
k k
A
R s
b a sen
ξ
θ
θ θ
=
+
(4.92)
A distância
( , )
ml
R s
θ
é a distância do ponto s a um ponto q pertencente ao lado
definido pelos pontos
e
m l
q q
, quando o triângulo é definido pelos pontos
, e
m l k
q q q
como na Figura 4-22.
Os limites de integração (
e
m M
R R
) são dados pela expressão (4.92). Sempre
haverá um dos extremos representados por duas equações diferentes, devendo-se
dividir o intervalo da integração em dois.
Figura 4-22 - Extremos de integração para ponto fonte.
No caso da Figura 4-22, tem-se:
1
s
l
q
m
q
k
q
l
θ
k
θ
m
θ
ml
R
km
R
107
ln
~ ~
ln
~ ~
( ) ( )
( ) ( )
R R
kl
o n
ki ki ki
k l
R R
kn kn
k
R R
kl
l M
ki ki
k l
R R
kn kn
u K s d M d M d K s
N d N d
θ θ
θ θ
θ θ
θ θ
θ θ
ζ θ θ ζ θ θ
Ω
ª º
Ω = + +
« »
¬ ¼
+ +
³ ³ ³
³ ³
(4.93)
onde
1
s
é o ponto fonte. Desse modo, os valores extremos
, , ,
m M m M
R R
θ θ
ficam
devidamente definidos (Figura 4-23).
Figura 4-23 - Ponto fonte no vértice da célula.
Uma posição possível para o ponto fonte, não mostrada na figura anterior, é
aquela em que ele coincide com um dos cantos da célula (Figura 4-23). Na equação
integral (4.93), os extremos inferiores da variável R são iguais a zero, e os limites
abaixo devem ser considerados.
0
0
lim 0
lim 0
ki
r
ki
r
M
M
+
→
→
=
=
108
0
0
lim 0
lim 0
ki
r
ki
r
M
M
+
→
→
=
=
Chamando
1 1
e
ki ki
M N
as matrizes cujos extremos de integração são
e
kn kl
R R
e
2 2
e
ki ki
M N
quando os extremos são
ln
e
kn
R R
, pode-se escrever a equação (4.93)
como:
1 2
~ ~
1 2
( ) ( ) ( ) ( )
( ) ( ) ( ) ( )
l n
ki ki ki
k k l
l n
ki ki
k l
u K s d M d M d K s
N d N d
θ θ
θ θ
θ θ
θ θ
θ θ θ θ
θ ζ θ θ θ ζ θ θ
Ω
ª º
« »
Ω = + +
« »
¬ ¼
+ +
³ ³ ³
³ ³
(4.94)
A finalidade dessa forma compacta é facilitar a transformação das integrais em
θ
,
para integrais adimensionais nas variáveis
η η
1
e
2
, com:
[ ]
[ ]
1
l k
2
n l
-
2
- para o intervalo ,
- 2
2
- para o intervalo ,
- 2
l k
l k
n l
n l
θ θ
η θ θ θ
θ θ
θ θ
η θ θ θ
θ θ
§ ·
=
¨ ¸
© ¹
+
§ ·
=
¨ ¸
© ¹
(4.95)
As derivadas de
η η
1 2
e
são dadas por:
( )
( )
1
2
2
-
2
-
l k
n l
d d
d d
η θ
θ θ
η θ
θ θ
=
=
109
Assim, a equação (4.94) pode ser escrita:
1 1
1 1 2 2
~ ~
-1 -1
1 1
1 1 2 2
~ ~
-1 -1
- -
2 2
- -
2 2
l k n l
ki ki ki
k
l k n l
ki ki
u K d M d M d K
N d N d
θ θ θ θ
η η
θ θ θ θ
ζ η ζ η
Ω
ª º
Ω = + +
« »
« »
¬ ¼
+ +
³ ³ ³
³ ³
(4.96)
As integrais da equação (4.96) são calculadas pela quadratura de Gauss (Anexo
3).
Com os coeficientes obtidos, pode-se escrever a expressão como:
( )
*
~
~
1
2
11 12 11 12 11 12 1
21 22 21 22 21 22 2
1
2
( , ) ( ) ( )
k
ki
j
k
k
k l n l
n
u S q K q d q U
U
U
m m m m m m U
m m m m m m U
U
U
Ω
Ω =
½
½
° °
® ¾
¯ ¿
° °
° °
ª º
ª º ª º ª º ½
° °
=
« »
®® ¾ ¾
« » « » « »
« »
¬ ¼ ¬ ¼ ¬ ¼ ¯ ¿
° °
¬ ¼
° °
½
° °
® ¾
° °
¯ ¿
¯ ¿
³
(4.97)
onde as submatrizes
ij
m
são as submatrizes da matriz global M , da expressão (4.88)
4.2.4 OBTENÇÃO DA SOLUÇÃO NUMÉRICA
A equação (4.73) representa um sistema de equações diferenciais com
coeficientes constantes:
H U + M U = G P
110
Para se obter soluções para deslocamentos e forças, uma distinção entre os dois
tipos de condições de contorno deve ser feita:
1
Γ
será a posição do contorno em que os
deslocamentos são prescritos e
2
Γ
onde as forças são prescritas:
1
2
U
U =
U
½
® ¾
¯ ¿
1
2
P
P =
P
½
® ¾
¯ ¿
(4.98)
onde
1 2
U e P
são valores conhecidos.
Particionando as matrizes, tem-se:
11 1 12 2 11 1 12 2 11 1 12 2
H U H U M U M U G P G P
+ + + = +
(4.99)
21 1 22 2 21 1 22 2 21 1 22 2
H U H U M U M U G P G P
+ + + = +
Através da primeira das equações (4.99) escrevem-se as forças incógnitas
1
P
como função das outras variáveis:
11 1 11 1 12 2 11 1 12 2 12 2
1 1 1
1 11 11 1 11 12 2 11 11 1
1 1
11 12 2 11 12 2
G P = H U + H U + M U + M U - G P
P = G H U + G H U + G M U +
+ G M U - G G P
− − −
− −
(4.100)
Rearranjando a segunda equação de (4.99):
22 2 22 2 21 1 22 2 21 1 21 1
H U M U - G P = G P - H U - M
U
+
(4.101)
111
Substituindo
1
P
na expressão (4.101):
1 1
22 2 22 2 21 11 11 1 11 12 2
1 1 1
11 11 1 11 12 2 11 12 2
22 2 21 1 21 1
H U + M U - G ( G H U + G H U
+ G M U + G M U - G G P )
= G P - H U - M U
− −
− − −
=
(4.102)
Agrupando as submatrizes de maneira a se ter do lado direito os deslocamentos
e forças conhecidos, tem-se:
1 1
22 21 11 12 2 22 21 11 12 2
1 1
22 21 11 12 2 21 11 11 21 1
1
21 11 11 21 1
( - ) ( - )
( - ) ( )
( )
H G G H U M G G M U
G G G G P G G H H U
G G M M U
− −
− −
−
+ =
= + − +
+ −
(4.103)
De forma compacta:
2 2 2 1 1
ˆ ˆ
ˆ
ˆ ˆ ˆ ˆ
H U M U G P H U M U
+ = + +
(4.104)
Com:
1
22 21 11 12
1
22 21 11 12
ˆ
-
ˆ
-
H H G G H
M M G G M
−
−
=
=
112
1
22 21 11 12
1
21 11 11 21
1
21 11 11 21
ˆ
G = G - G G G
ˆ
ˆ
H = G G H - H
ˆ
ˆ
M = G G M - M
−
−
−
Para resolver o sistema de equações (4.104) será utilizado a integração no
tempo pelos Métodos de Newmark e Houbolt, descridos no Capítulo 8.
No problema de vibrações livres, onde
2
0
P
=
, a equação (4.104) torna-se:
ˆ ˆ
0
H U M U
+ =
(4.105)
Usando o mesmo procedimento descrito na formulação estática (item 6.1.5), os
deslocamentos e acelerações nodais são separados em três tipos:
( )
M M
U U
= deslocamentos (acelerações) em nós duplos considerados “master”;
( )
D D
U U
= deslocamentos (acelerações) em nós duplos considerados
“dependentes”;
( )
I I
U U
= deslocamentos (acelerações) em nós não duplos chamados
“independentes”.
A equação (4.105) pode ser reescrita como:
ˆ ˆ ˆ ˆ ˆ ˆ
0
ˆ ˆ ˆ ˆ ˆ ˆ
0
0
ˆ ˆ ˆ ˆ ˆ ˆ
DD DM DI DD DM DI
D
D
M
M
MD MM MI MD MM MI
I
I
ID IM II ID IM II
H H H M M M
U
U
U
U
H H H M M M
U
U
H H H M M M
ª º ª º
½
½
« » « »
° °
° °
°
« » « »
° °
° °
°
° ° ° ° °
« » « »
+ =
® ¾ ® ¾ ®
« » « »
° ° ° °
« » « »
° ° ° °
« » « »
° ° ° °
¯
¯ ¿
¯ ¿
« » « »
¬ ¼ ¬ ¼
½
°
°
°
¾
° °
° °
° °
¿
(4.106)
113
Através das relações:
D M
U U
=
e
D M
U U
=
O sistema de equações (4.106) pode ser condensado, eliminando-se os
deslocamentos e acelerações dependentes:
DD MD DM MM DI MI
M
I
ID IM II
DD MD DM MM DI MI
M
I
ID IM II
ˆ ˆ ˆ ˆ ˆ ˆ
H H H H H H
U
U
ˆ ˆ ˆ
H H H
ˆ ˆ ˆ ˆ ˆ ˆ
M M M M M M
0
U
0
U
ˆ ˆ ˆ
M M M
ª º
+ + + +
½« »
+
® ¾
« »
¯ ¿
« »
+
¬ ¼
ª º
+ + + +
½
« »
½
° °
+ =
® ¾ ® ¾
« »
° °
¯ ¿
¯ ¿
« »
+
¬ ¼
(4.107)
De forma compacta:
ˆ ˆ
0
C C C C
H U M U+ =
(4.108)
A solução harmônica para a equação (4.108) é da forma:
-
0
2 -
0
-
i t
i t
U U e
U U e
ω
ω
ω
=
=
(4.109)
onde
ω
é a freqüência.
Substituindo a equação (4.109) em (4.108):
- 2 -
0 0
ˆ ˆ
- 0
i t i t
H U e M U e
ω ω
ω
=
Simplificando:
114
2
0 0
ˆ ˆ
H U M U
ω
=
(4.110)
Colocando a equação (4.110) na forma padrão para cálculo de autovalores e
autovetores, obtém-se:
1
0 0
2
1
ˆ ˆ
H M U U
ω
−
=
(4.111)
Para o cálculo dos autovalores e autovetores foi utilizado o software Mathematica
(SOLIANI e BARRETO, 1997)
.
115
CAPÍTULO 5 MÉTODO DOS ELEMENTOS FINITOS
(MEF)
116
5.1 FORMULAÇÃO ESTÁTICA
Apresenta-se a composição básica do Método dos Elementos Finitos com
aplicação para meios contínuos da Mecânica dos Corpos Sólidos. Isto se faz a partir da
expressão de equilíbrio dos corpos elásticos, obtida pela aplicação do Princípio dos
Trabalhos Virtuais (PTV).
Para um corpo contínuo, elástico linear, pode-se escrever a energia potencial total,
funcional do problema, como:
1
- - -
2
T T
T T
i i
d U b d U f d U F
ε σ
Ω Ω Γ
Π = Ω Ω Γ
¦
³ ³ ³
(5.1)
onde:
ε
= vetor de deformações;
σ
= vetor de tensões;
U = vetor deslocamento de um ponto qualquer;
b = forças de volume;
f = forças distribuídas na superfície;
i
F
= força concentrada em um ponto qualquer.
O PTV estabelece que para qualquer campo de pequenos deslocamentos virtuais
compatíveis imposto sobre o corpo, o trabalho total das forças virtuais internas deve ser
igual ao trabalho total das forças externas, ou seja:
0
δ
Π =
(5.2)
117
Essa expressão representa a solução exata do problema contínuo, sendo
necessário encontrar as funções que definem o comportamento das variáveis
incógnitas, para as quais o problema é formulado. Ela pode ser expressa pela equação:
T T T T
i i
i
d U b d U f d U F
ε σ
Ω Ω Γ
Ω = Ω + Γ +
¦
³ ³ ³
(5.3)
onde a ´barra superior´ representa o estado devido a pequenos deslocamentos virtuais,
que satisfazem as condições de contorno do problema.
No Método dos Elementos Finitos, o funcional exato
Π
é substituído por um
funcional aproximado
a
Π
, onde as variáveis do problema são expressas em termos de
funções de interpolação. Para isso, subdivide-se o domínio de integração em n
elementos finitos, interligados pelos
N pontos nodais.
A expressão (5.3) pode, então ser escrita como:
( )
( ) ( )
( ) ( ) ( )
( ) ( ) ( ) ( ) ( ) ( )
T
n
T T
n n
n n n
n
T
n n n n n n
i
n i
i
d
U b d U f d U F
ε σ
Ω
Ω Γ
ª º
Ω =
¦
³
« »
¬ ¼
ª º
= Ω + Γ +
¦ ¦
³ ³
« »
¬ ¼
(5.4)
Os deslocamentos que ocorrem dentro de cada elemento são aproximados por
funções polinomiais e valores nodais, sendo que para o elemento n pode-se escrever:
( ) ( ) ( )
ˆ
( ) ( ) .
n n n
U q q U
φ
=
(5.5)
onde:
( )
( )
n
q
φ
=
matriz de interpolação dos deslocamentos;
118
( )
ˆ
n
U
= vetor que contém os deslocamentos nodais.
Do mesmo modo, as deformações também podem ser escritas em função dos
deslocamentos nodais, através da relação:
( ) ( ) ( )
ˆ
( ) ( )
n n n
q B q U
ε
=
(5.6)
onde:
( )
( )
n
B q
= matriz que relaciona os deslocamentos com as deformações.
A relação tensão-deformação é escrita através da Lei de Hooke:
( ) ( ) ( )
( )
n n n
q C
σ ε
=
(5.7)
Para o caso de estado plano de deformação específica, considerando material
elástico homogêneo, isótropo e com comportamento elástico e linear, a matriz
C
vale:
1 - 0
1 - 0
(1 ) (1- 2 )
1 - 2
0 0
2
E
C
ν ν
ν ν
ν ν
ν
ª º
« »
« »
=
« »
+
« »
« »
¬ ¼
(5.8)
com:
E = Módulo de Elasticidade;
ν
= Coeficiente de Poisson.
As equações (5.6) e (5.7) fornecem:
119
( ) ( ) ( ) ( )
ˆ
( ) ( )
n n n n
q C B q U
σ
=
(5.9)
Substituindo (5.5), (5.6) e (5.9) na equação (5.4), obtém-se:
( )
( ) ( )
( ) ( ) ( ) ( )
( ) ( ) ( ) ( ) ( ) ( )
ˆ
ˆ
T
n
T T
n n
T n n n n
n
T n n n n n n
C
n n
U B C B d U
U b d f d F U
φ φ
Ω
Ω Γ
ª º
Ω =
¦
³
« »
¬ ¼
ª º
½ ½
= Ω + Γ +
¦ ¦
³ ³
® ¾ ® ¾
« »
¯ ¿ ¯ ¿
¬ ¼
(5.10)
Para cada elemento finito n define-se um sistema de coordenadas locais (Figura
5-1), ao qual são referenciadas as variáveis nodais. Para se ter válidos os somatórios
da equação (5.10), é necessário rotacionar as matrizes do sistema de coordenadas
locais para o global, antes de efetuar a soma.
Figura 5-1 - Sistema local e global de coordenadas.
A matriz de transformação
T
, que relaciona o sistema global ao sistema local e
vice-versa, é dada por:
~
cos
cos
sen
T
sen
θ θ
θ θ
ª º
=
« »
« »
−
¬ ¼
(5.11)
com o ângulo
θ
definido na Figura 5-1.
( )
x u
( )
y v
θ
1
2
y
x
120
Assim:
( )
( )
= .
= .
n
T n
U T U
U T U
(5.12)
Fazendo os deslocamentos virtuais unitários, cada um a seu próprio tempo,
T
U
torna-se uma matriz identidade. Sendo
ˆ
U
, simplesmente
U
, a equação (5.10) pode ser
escrita como :
K U F
=
(5.13)
onde:
K
é a matriz de rigidez, dada por:
( )
( )
( ) ( ) ( ) ( )
T
n
n
n n n n
n n
K B C B d K
Ω
= Ω =
¦ ¦
³
(5.14)
F
é o vetor de forças nodais equivalentes, dado por:
C
F F f f
Ω Γ
= + +
(5.15)
( )
( )
( ) ( ) ( ) ( )
( ) ( ) ( ) ( )
T
n
T
n
n n n n
n n
n n n n
n n
f b d f
f f d f
φ
φ
Ω Ω
Ω
Γ Γ
Γ
= Ω =
¦ ¦
³
= Γ =
¦ ¦
³
Esse procedimento permite a montagem do sistema de equações considerando
elemento por elemento.
A equação (5.13) é resolvida pelas técnicas usuais do cálculo numérico.
121
5.2 FORMULAÇÃO DINÂMICA
Torna-se necessário considerar as forças de inércia na obtenção da equação de
equilíbrio de um corpo, quando as forças externas são aplicadas repentinamente.
Podem-se incluir as forças de inércia como parte das forças de volume, usando-se o
Princípio de D’ Albert.
Aproximando-se a aceleração pela mesma função utilizada para os
deslocamentos (equação(5.5)), escreve-se:
( )
( )
( ) ( ) ( ) ( )
- -
T
n
n
n n n n
n
f b U d f M U
φ ρ φ
Ω Ω
Ω
ª º
= Ω =
¦
³
¬ ¼
(5.16)
com:
( )
( ) ( ) ( ) ( ) ( )
T
n
n n n n n
n n
M d M
ρ φ φ
Ω
= Ω =
¦ ¦
³
onde:
( )
n
ρ
= densidade do material no elemento n ;
M
= matriz de massa;
A equação de equilíbrio fica:
K U M U F
+ =
(5.17)
Os problemas dinâmicos a serem analisados nesse trabalho são transientes.
Logo, a integração no tempo da equação (5.17) será feita numericamente. Foram
utilizados dois algoritmos Newmark e Houlbolt, que serão apresentados no Capítulo 8.
122
5.3 ELEMENTO FINITO UNIDIMENSIONAL
Nesse elemento, as incógnitas nodais são (Figura 5-2):
•
duas translações : os deslocamentos (u
e v) ;
•
uma rotação (
θ
).
Figura 5-2 - Elemento finito unidimensional.
Admitindo-se a deformação longitudinal
ε
constante ao longo do elemento finito,
os deslocamentos longitudinais podem ser representados por uma reta:
1 2
u a a x
= +
(5.18)
Então, os deslocamentos transversais têm que ser representados por um
polinômio o terceiro grau, do tipo:
2 3
3 4 5 6
v a a x a x a x
= + + +
(5.19)
Como os deslocamentos transversais v
e as rotações
dv
dx
θ
§ ·
=
¨ ¸
© ¹
não são
independentes, a componente v do deslocamento em qualquer ponto do elemento
1
u
1
1
2
u
1
v
2
v
1
θ
2
θ
L
123
finito, depende de quatro componentes de deslocamentos:
1 1
e
v
θ
no nó 1 e
2 2
e
v
θ
no nó 2 do elemento finito.
É mais conveniente trabalhar com funções aproximadoras explicitadas em relação
às incógnitas nodais. Assim, das igualdades (5.18) e (5.19) tem-se as funções
aproximadoras para o elemento finito em estudo:
( ) ( ) ( )
ˆ
n n n
U U
φ
=
(5.20)
com:
( )
n
u
U
v
½
° °
=
® ¾
° °
¯ ¿
( )
( )
2 3
2 3
( )
2 3
3 2
1 0
0 1 - 3 2
0 - 2
0
0 3 - 2
0 -
T
n
L
L
ξ
ξ ξ
ξ ξ ξ
ε
ξ
ξ ξ
ξ ξ
−
ª º
« »
« »
+
« »
« »
+
« »
=
« »
« »
« »
« »
« »
« »
¬ ¼
{
}
( )
1 1 1 2 2 2
ˆ
T
n
U u v u v
θ θ
=
(5.21)
Assumindo-se as funções de Bernoulli e desprezada a deformação por força
cortante, as matrizes de rigidez e massa passam a ser expressas em coordenadas
locais como:
124
~
1 0 0 1 0 0
0 2 3 0 2 3
0 3 4 0 3 5
1 0 0 1 0 0
0 2 3 0 2 3
0 3 5 0 3 4
A A
A A A A
A A A A
k
A A
A A A A
A A A A
−
ª º
« »
« »
−
« »
« »
−
« »
=
« »
−
« »
« »
− − −
« »
« »
−
¬ ¼
(5.22)
com:
3
2
1
12
2
6
3
4
4
2
5
E A
A
L
E I
A
L
E I
A
L
E I
A
L
E I
A
L
=
=
=
=
=
sendo:
A = área da seção transversal;
E = módulo de elasticidade longitudinal;
L = comprimento do elemento;
I = momento de inércia da seção transversal.
125
~
1 0 0 2 0 0
0 3 4 0 5 6
0 4 7 0 6 8
420
2 0 0 1 0 0
0 5 6 0 3 4
0 6 8 0 4 7
B B
B B B B
B B B B
A L
M
B B
B B B B
B B B B
ρ
ª º
« »
« »
−
« »
« »
−
« »
=
« »
« »
« »
−
« »
« »
− − −
¬ ¼
(5.23)
com:
B1 = 140
B2 = 70
B3 = 156
B4 = 22 L (5.24)
B5 = 54
B6 = 13 L
B7 = 4
2
L
B8 = 3
2
L
Para se obter um vetor de carga nodal equivalente, que seja global, supõe-se o
elemento finito inclinado de um ângulo
α
com a direção horizontal, solicitado por carga
vertical linearmente variável, conforme Figura 5-3.
126
Figura 5-3- Elemento inclinado com carga variável.
Decompondo o carregamento da Figura 5-3 em uma carga perpendicular ao
elemento, e outra tangente ao mesmo, obtém-se:
(
)
( )
1 2
2
1 2
- (1 - ) cos - cos
- (1 - ) cos - cos
u
v
f f sen f sen
f f f sen
ξ ξ α α ξ α α
ξ ξ α ξ α α
=
=
(5.25)
De forma matricial tem-se:
(
)
( )
1
2 2
2
-(1- ) cos - cos
-(1- )cos - cos
u
v
f
f
sen sen
f
f
ξ
ξ α α ξ α α
ξ
ξ α ξ α
½
ª º
½
° ° ° °
« »
=
® ¾ ® ¾
« »
° ° ° °
¯ ¿
¬ ¼
¯ ¿
(5.26)
O vetor de cargas nodais equivalentes é calculado pela expressão (5.15), na qual
a matriz
( )
T
n
ε
é dada pela equação (5.21).
Da mesma maneira, se houver forças volumétricas, como o peso próprio, pode-se
utilizar a igualdade (5.15).
( )
y v
1
2
α
( )
x u
≡
1
1
2
2
α
α
+
127
As cargas concentradas nos nós dos elementos finitos são diretamente
adicionadas ao vetor de cargas nodais equivalentes, considerando seus sinais corretos,
nas posições correspondentes aos graus de liberdade a que se relacionam.
128
CAPÍTULO 6 MÉTODOS DE INTEGRAÇÃO
NUMÉRICA
129
6.1 INTRODUÇÃO
A maneira mais usual para se obter uma resposta na análise dinâmica de um
sistema estrutural é utilizar um processo de integração direta no tempo. Esse processo
requer que o equilíbrio dinâmico do sistema seja obedecido em instantes discretos de
tempo, durante os quais são admitidas variações lineares ou constantes para os
deslocamentos, velocidades e acelerações.
Muitas possibilidades existem de métodos de integração direta no tempo. Esses
métodos podem ser implícitos e explícitos.
Métodos implícitos buscam soluções que satisfaçam as equações diferenciais no
tempo
t
, sendo conhecidas as soluções num tempo
t t
− ∆
. Esses métodos exigem que
o sistema de equações seja resolvido a cada passo de tempo. Porém intervalos
maiores de tempo podem ser utilizados.
Métodos explícitos não envolvem a solução de um conjunto de equações em
cada passo de tempo. Basicamente esses métodos usam a equação no tempo t para
prever uma solução num tempo
t t
+ ∆
. Para muitos estruturas, compostas de
elementos rígidos, é necessário que o intervalo de tempo
t
∆
seja muito pequeno para
se obter uma solução estável.
Foram utilizados nesse trabalho os métodos implícitos de Newmark e de
Houlbolt.
6.2 MÉTODO DE NEWMARK
A equação de equilíbrio dinâmico a ser resolvida é:
t t t t t t t t
M U C U K U F
+∆ +∆ +∆ +∆
+ + =
(6.1)
130
onde:
M
= matriz de massa;
C
= matriz de amortecimento;
K
= matriz de rigidez;
F
= vetor de cargas nodais;
U
= vetor deslocamento;
U
= vetor velocidade;
U
= vetor aceleração.
Sendo conhecidos o deslocamento
t
U
e a velocidade
t
U
de um ponto material no
instante
t
, pode-se determinar sua posição e velocidade no instante
t t
+ ∆
pelas
seguintes expressões:
( )
t t
t t t
t
U U U d
τ τ
+∆
+∆
= +
³
(6.2)
( ) ( )
t t
t t t t
t
U U t U t t U d
τ τ τ
+∆
+∆
= + ∆ + + ∆ −
³
(6.3)
Aproximando-se em séries de potência, a função aceleração
(
)
U
τ
pode ser
escrita:
131
( ) ( )
1
t t
t t t
t
U d U U t
τ τ γ γ
+∆
+∆
ª º
≅ − + ∆
¬ ¼
³
(6.4)
( ) ( )
2
1
2
t t
t t t
t
t t U d U U t
τ τ τ β β
+∆
+∆
ª º
§ ·
+ ∆ − ≅ − + ∆
¨ ¸
« »
© ¹
¬ ¼
³
(6.5)
Substituindo-se esses valores aproximados das integrais nas equações (6.2) e
(6.3) chega-se a:
(
)
1
t t t t t t
U U U U t
γ γ
+∆ + ∆
ª º
= + − + ∆
¬ ¼
(6.6)
2
1
2
t t t t t t t
U U t U U U t
β β
+∆ +∆
ª º
§ ·
= + ∆ + − + ∆
¨ ¸
« »
© ¹
¬ ¼
(6.7)
O integrador de Newmark é incondicionalmente estável se:
1
2
γ α
= +
e
2
1 1
4 2
β γ
§ ·
= +
¨ ¸
© ¹
para
0
≥
α
(6.8)
Se
0
α
=
na equação (8.8),
4
1
=
β
, o que nos mostra que a aceleração durante o
intervalo de tempo
t
∆
é constante e igual à média de
t t t
U e U
+∆
.
( )
1
0
2
t t t
U U U t
τ τ
+∆
ª º
= + ≤ ≤ ∆
¬ ¼
(6.9)
Com esses valores, obtém-se uma simplificação para as expressões (6.6) e (6.7):
1
2
t t t t t t
U U t U U
+∆ +∆
ª º
= + ∆ +
¬ ¼
(6.10)
132
( )
( )
2
1
4
t t t t t t t
U U t U t U U
+∆ +∆
= + ∆ + ∆ +
(6.11)
Isolando-se o valor de
t t
U
+∆
na equação (6.11), chega-se a:
( )
2
4 4
t t t t t t t
U U U U U
t t
+∆ +∆
= − − −
∆ ∆
(6.12)
Substituindo na equação (6.10), obtém-se:
( )
2
-
t t t t t t
U U U U
t
+∆ +∆
= −
∆
(6.13)
Na equação de equilíbrio dinâmico (6.1), substituindo-se os valores da aceleração
(6.12) e da velocidade (6.13) obtém-se:
( ) ( )
( ) ( )
2 2
2 2
1 1 1
2 4 4
1 1 1
2 4 4
t t t t
t t t
M t C t K U t F
M t C U t M t C U t MU
+∆ +∆
ª º
+ ∆ + ∆ = ∆ +
« »
¬ ¼
ª º ª º
+ + ∆ + ∆ + ∆ + ∆
« » « »
¬ ¼ ¬ ¼
(6.14)
O sistema (6.14) pode ser escrito de forma simplificada como:
t t
K U F
+∆
=
(6.15)
com:
( )
( ) ( ) ( )
2
2 2 2
1 1
2 4
1 1 1 1
4 2 4 4
t t t t t
K M t C t K
F t F M t C U t M t C U t MU
+∆
= + ∆ + ∆
°
°
°
®
°
ª º ª º
°
= ∆ + + ∆ + ∆ + ∆ + ∆
« » « »
°
¬ ¼ ¬ ¼
¯
(6.16)
133
Para a determinação de
t t
U
+∆
:
•
Monta-se a matriz de rigidez
K
, a matriz de massa
M
e a matriz de amortecimento
C
;
•
Determina-se a matriz
K
;
• Conhecidos os valores de
0 0 0
, e
U U F
, determina-se o valor das acelerações no
instante inicial
0
t
, pela equação:
0 0 0 0
- -
M U F C U K U
=
(6.17)
Para cada passo de tempo
t
∆
:
1. Calcula-se o vetor
F
;
2. Determina-se
t t
U
+∆
resolvendo-se o sistema de equações (6.15);
3. Determina-se o valor de
t t
U
+∆
com a expressão (6.13);
4. Determina-se o valor de
t t
U
+∆
com a expressão (6.12);
5. Retorna-se ao passo 1 para o próximo tempo
t t
+ ∆
.
6.3 MÉTODO DE HOUBOLT
O esquema de integração de Houbolt produz bons resultados para quase todos os
problemas. É um algoritmo implícito e incondicionalmente estável. É chamado de
134
Integrador de Passo Múltiplo, pois depende de mais um passo de tempo para a
determinação dos valores do passo seguinte.
A velocidade e a aceleração em um instante
t
são escritas como:
( )
2
1
11 18 9 2
6
t t t t t t t t t
U U U U U
t
+∆ +∆ −∆ − ∆
= − + −
∆
(6.18)
( )
( )
2
2
1
2 5 4
t t t t t t t t t
U U U U U
t
+∆ + ∆ −∆ − ∆
= − + −
∆
(6.19)
Substituindo-se as equações (6.18) e (6.19) na equação de equilíbrio (6.1), obtém-
se a expressão para o esquema de avanço no tempo:
( )
( )
2
2
2
2 11
3
5 3 4
2 3
t t
t t t t t t t
M t C t K U
t t
t F M t C U M C U M C U
+∆
+∆ −∆ − ∆
ª º
+ ∆ + ∆ =
« »
¬ ¼
∆ ∆
§ · § ·
= ∆ + + ∆ − + + +
¨ ¸ ¨ ¸
© ¹ © ¹
(6.20)
O intervalo de discretização no tempo
t
∆
deve obedecer a:
min
s
L
t
c
∆
∆ ≥ (6.21)
com
s
c
dado pela expressão (4.8) e
L
∆
sendo o menor comprimento do elemento de
contorno usado na elaboração da malha.
A obtenção dos deslocamentos nos dois primeiros passos, necessários para o
início da análise, normalmente é feita pela utilização de outro integrador. Daí em diante
o integrador de Houbolt é aplicado diretamente.
O sistema de equações (6.20) pode ser escrito de forma simplificada como :
135
t t
K U F
+∆
=
sendo:
( )
( )
2
2
2
2 11
3
5 3 4
2 3
t t t t t t t
K M t C t
t t
P t F M t C U M C U M C U
+∆ −∆ − ∆
°
= + ∆ + ∆
°
°
®
°
∆ ∆
§ · § ·
°
= ∆ + + ∆ − + + +
¨ ¸ ¨ ¸
°
© ¹ © ¹
¯
(6.22)
Para a determinação de
t t
U
+∆
:
•
Monta-se a matriz de rigidez
K
, a matriz de massa
M
e a matriz de amortecimento
C
;
• Determina-se a matriz
K
;
•
Conhecidos os valores de
0 0 0
, e
U U F
, determina-se o valor das acelerações no
instante inicial
0
t
, pela equação:
0 0 0 0
- -
M U F C U K U
=
;
• Com outro integrador, determinam-se os valores dos deslocamentos, velocidades e
acelerações nos dois primeiros passos.
Para cada passo de tempo
t
∆
:
1. Calcula-se o vetor
F
;
2. Determina-se
t t
U
+∆
resolvendo-se o sistema de equações (6.20);
136
3. Determina-se o valor de
t t
U
+∆
com a expressão (6.18);
4. Determina-se o valor de
t t
U
+∆
com a expressão (6.19);
5. Retorna-se ao passo 1 para o próximo tempo
t t
+ ∆
.
6.4 APLICAÇÃO DOS INTEGRADORES AO MÉTODO DOS
ELEMENTOS DE CONTORNO
No Método dos Elementos de Contorno há um crescimento do número de
incógnitas do problema, uma vez que a matriz de massa leva em conta não somente a
aceleração nos nós do contorno, mas também nos pontos internos. A equação para o
problema discretizado e sem amortecimento assume a forma:
t t t t
HU GP MU F
= − +
(6.23)
O sistema representado pela equação (6.23) pode ser escrito considerando-se
sub-matrizes referentes ao pontos do contorno (c) e aos pontos internos do domínio (i):
0 0
0
0
tc cc cic c
tc
tc
t
ti ic iii i
ti
U M MH G
P
U
F
U M M
H I G
U
§ ·
§ · § ·§ · § ·
§ ·
= − +
¨ ¸
¨ ¸ ¨ ¸¨ ¸ ¨ ¸
¨ ¸
© ¹
© ¹ © ¹© ¹ © ¹
© ¹
(6.24)
No instante inicial
0
0
t
=
são conhecidos os deslocamentos e velocidades nos
nós do contorno e pontos internos, bem como as forças de superfície e as forças de
corpo. As acelerações podem ser calculadas através da equação (6.23) , resolvendo-se
o sistema:
0 0 0 0
MU GP HU F
= − +
(6.25)
137
6.4.1 MÉTODO DE NEWMARK APLICADO AO MEC
Para o instante
t t
+ ∆
, a equação (6.23) pode ser escrita como:
t t t t t t t t
HU GP MU F
+∆ +∆ +∆ +∆
= − +
(6.26)
Substituindo-se o valores da aceleração (6.12) obtém-se:
( )
( )
( )
[ ]
( )
2
2 2 2
1
4
1 1
4 4
t t
t t t t t t t
M t H U
t GP t F M U t U t U
+∆
+∆ +∆
ª º
+ ∆ =
« »
¬ ¼
ª º
= ∆ + ∆ + + ∆ + ∆
« »
¬ ¼
(6.27)
O sistema (6.27) pode ser escrito de forma simplificada como:
(
)
(
)
2 2
4
t t t t t t t
HU t GP MU t F
+∆ +∆ +∆
= ∆ + + ∆
(6.28)
onde:
( )
( )
2
2
4
4 4
t t t
H t H M
U U t U t U
= ∆ +
= + ∆ + ∆
(6.29)
Para a determinação de
t t
U
+∆
:
•
Monta-se a matriz
H
e a matriz de massa
M
;
•
Determina-se a matriz
H
;
138
•
Conhecidos os valores dos deslocamentos, velocidades e forças nos nós do
contorno, determina-se o valor das acelerações no instante inicial
0
t
, pela equação
(6.25).
Para cada passo de tempo
t
∆
:
1. Calcula-se o vetor
U
;
2. Determina-se
t t
U
+∆
resolvendo-se o sistema de equações (6.28);
3. Determina-se o valor de
t t
U
+∆
com a expressão (6.13);
4. Determina-se o valor de
t t
U
+∆
com a expressão (6.12);
5. Retorna-se ao passo 1 para o próximo tempo
t t
+ ∆
.
6.4.2 MÉTODO DE HOUBOLT APLICADO AO MEC
Substituindo na equação (6.26) o valor de
t t
U
+∆
dado por (6.19) obtém-se:
( )
(
)
( ) ( ) ( )
2
2 2
2
2
5 4
t t
t t t t t t t t t
t H M U
t GP M U U U t F
+∆
+∆ −∆ − ∆ +∆
∆ + =
= ∆ − − + − + ∆
(6.30)
O sistema (6.30) pode ser escrito de forma simplificada como:
(
)
(
)
2 2
t t t t t t t
HU t GP MU t F
+∆ +∆ +∆
= ∆ − + ∆
(6.31)
onde:
139
( )
2
2
2
5 4
t t t t t t
H t H M
U U U U
−∆ − ∆
= ∆ +
= − + −
(6.32)
Para a determinação de
t t
U
+∆
:
•
Monta-se a matriz
H
e a matriz de massa
M
;
• Determina-se a matriz
H
;
•
Conhecidos os valores dos deslocamentos, velocidades e forças nos nós do
contorno, determina-se o valor das acelerações no instante inicial
0
t
, pela equação
(6.25).
• Com outro integrador, determinam-se os valores dos deslocamentos, velocidades e
acelerações nos dois primeiros passos.
Para cada passo de tempo
t
∆
:
1. Calcula-se o vetor
U
;
2. Determina-se
t t
U
+∆
resolvendo-se o sistema de equações (6.28);
3. Determina-se o valor de
t t
U
+∆
com a expressão (6.13);
4. Determina-se o valor de
t t
U
+∆
com a expressão (6.12);
5. Retorna-se ao passo 1 para o próximo tempo
t t
+ ∆
.
140
CAPÍTULO 7 COMBINAÇÃO DOS MÉTODOS
ELEMENTOS FINITOS - ELEMENTOS
DE CONTORNO
141
7.1 INTRODUÇÃO
Muitos problemas de engenharia apresentam interações ou acoplamento entre
diferentes partes. Vibrações em fundações de máquinas, abalos sísmicos sobre
edificações, conforto em edifícios nas proximidades de metrôs são exemplos de
sistemas que podem ser acoplados num mesmo problema. Cada parte é representada
por uma região física sobre a qual soluções numéricas particulares são aplicadas.
A idéia de combinar diferentes métodos numéricos é muito interessante e
apropriada em determinados tipos de problemas. Escolher qual o método mais
adequado a cada tipo de problema não é uma tarefa muito fácil.
O Método dos Elementos de Contorno (MEC) tem se mostrado muito
conveniente na análise de problemas com concentrações de tensões. No caso de
domínio infinito, o uso de soluções fundamentais que obedecem as condições no
infinito torna-se uma das maiores vantagens do MEC. Uma das principais
características desse método é a simplicidade dos dados de entrada que descrevem o
problema homogêneo.
Por outro lado, o MEF é mais indicado para problemas com domínio limitado, não
homogêneo ou anisotrópico e também para o estudo da não linearidade.
A técnica de acoplamento MEC – MEF é usada para usufruir efetivamente das
vantagens de cada um dos dois métodos.
As diferenças entre os diversos procedimentos para combinar os dois métodos
estão concentradas no tratamento dado à montagem do sistema de equações. Esses
sistemas de equações gerados pelo MEC e pelo MEF são expressos em termos de
diferentes variáveis, que necessitam de algumas considerações para serem acopladas.
142
Nesse trabalho, a estrutura reticulada, tratada por elementos finitos, vai ser
acoplada à estrutura bidimensional, tratada por elementos de contorno, através dos
pontos da interface (MEC equivalente).
O acoplamento é realizado utilizando-se a técnica das sub-regiões a partir da
imposição das condições de equilíbrio e compatibilidade geométrica aplicadas na
superfície de contato entre os meios, chamada de interface.
Por causa das vantagens e desvantagens que cada método apresenta, muitos
autores têm se interessado no desenvolvimento de algoritmos onde o acoplamento
MEC-MEF é considerado. Pode-se citar: Vallabhan (1987), Sousa (1992), Araújo
(1994), Coda e Venturini (2000), Ferro (1999), Almeida (2003), Ribeiro (2005) entre
outros.
7.2 EQUAÇÕES BÁSICAS DA FORMULAÇÃO MEC E MEF
Para uma região discretizada pelo Método dos Elementos de Contorno, o
sistema de equações a ser resolvido é dado por (8.28), com a aplicação do Integrador
de Newmark, ou (8.31) com a aplicação do Integrador de Houbolt. Uma forma
simplificada para esse sistema matricial é:
HU GP A
= +
(7.1)
onde o vetor
A
contém os valores independentes e os vetores
U
e
P
são as variáveis
atuais do problema.
O Método dos Elementos Finitos possui o sistema algébrico dado por:
*
KU MU G P B F
+ = + +
(7.2)
143
onde:
K
= matriz de rigidez;
M
= matriz de massa;
U
= deslocamentos;
U
= acelerações;
B
= força de volume equivalente;
F
= força concentrada;
*
G P
= forças distribuídas equivalentes.
O produto da matriz de transformação
*
G
pela matriz das forças de superfície
P
é escrito na sua forma expandida para poder considerar as forças de contato como
incógnitas do problema. A matriz de transformação
*
G
depende das mesmas funções
de interpolação definidas para os deslocamentos utilizadas para representar as forças
de superfície na interface das regiões.
As equações definidas para o Método dos Elementos de Contorno (equação
(7.1) e Método dos Elementos Finitos (equação (7.2)) são de naturezas diferentes e não
podem ser combinadas diretamente. Nesse trabalho, para combinar os dois métodos, a
montagem do sistema de equações é feita para a estrutura bidimensional, resolvida por
elementos de contorno. O sistema algébrico da parte do domínio discretizada por
elementos finitos é transformado em um sistema cuja forma é análoga à dos elementos
de contorno.
144
Considera-se a influência da estrutura reticulada, tratada por elementos finitos,
através das forças de superfície e deslocamentos dos pontos de contato entre as duas
estruturas.
Figura 7-1 – Coordenadas do elemento de contorno e do elemento finito.
7.3 TÉCNICA DE ACOPLAMENTO MEC – MEF
Seja um domínio formado por duas regiões
c
Ω
e
f
Ω
(Figura 7-2). A região
c
Ω
é
discretizada por elementos de contorno, e a região
f
Ω
por elementos finitos .
Figura 7-2 - Domínios
c
Ω
e
f
Ω
.
c
Ω
f
Ω
i
Γ
ff
Γ
cc
Γ
1
1
1
4
3
2
2
4
5
6
2
2
3
1
elemento
de
contorno
elemento
finito
145
As duas regiões tem em comum a interface
i
Γ
, sendo:
c cc i
f ff i
Γ = Γ + Γ
Γ = Γ + Γ
(7.3)
As equações (7.1) e (7.2) em um instante
0
t
, podem ser colocadas na forma:
H U G P
=
(7.4)
Escrevendo a equação (7.4) para o domínio
c
Ω
, tem-se:
c c
c c c c c
i i
c c
i i
U P
H H G G A
U P
½ ½
° ° ° °
ª º ª º
= +
® ¾ ® ¾
¬ ¼ ¬ ¼
° ° ° °
¯ ¿ ¯ ¿
(7.5)
Expressão semelhante obtém-se quando se escreve a equação (7.4) para o
domínio
f
Ω
:
f f
f f f f f
i i
f f
i I
U P
H H G G A
U P
½ ½
° ° ° °
ª º ª º
= +
® ¾ ® ¾
¬ ¼ ¬ ¼
° ° ° °
¯ ¿ ¯ ¿
(7.6)
O acoplamento da estrutura reticulada na estrutura bidimensional envolve
apenas deslocamentos e forças de superfície nos pontos da interface. Assim, o índice i
corresponde às coordenadas de translação dos pontos que estão na interface.
146
A compatibilidade de deslocamentos e as condições de equilíbrio em
i
Γ
são
dadas respectivamente por :
-
c f
i i i
c f
i i i
U U U
P P P
= =
= =
(7.7)
Com as equações (7.5), (7.6) e (7.7) chega-se a:
0 - 0
0 0
c
c c c c c c
f
i i
f f f f f f
i i
i
i
U
H H G G P A
U
H H G G P A
P
U
½
° °
° °
° °
ª º ª º ½ ½
° °
« » « » ° ° ° °
° °
= +
® ¾ ® ¾ ® ¾
« » « »
° ° ° ° ° °
« » « »
¬ ¼ ¬ ¼ ¯ ¿ ¯ ¿
° °
° °
° °
° °
¯ ¿
(7.8)
Se houver carregamento prescrito na interface i
(
)
ou
c f
i i
P P
, a equação (7.8)
fica:
c c c
i i i
f f f c
i i ii i
f
i
i
0 - 0 0
0 0 0
c
c
c c c
f
f
f f f
P
U
H H G G G A
P
U
H H G G G A
P P
U
P
½
½
° °
° °
° °
° °
° °
° °
ª º ª º ½
° °
° °
« » « » ° °
° ° ° °
= +
® ¾ ® ¾ ® ¾
« » « »
° ° ° ° ° °
« » « »
¬ ¼ ¬ ¼ ¯ ¿
° ° ° °
° ° ° °
° ° ° °
° °
° °
¯ ¿
¯ ¿
(7.9)
147
Resolve-se o sistema de equações para determinar as incógnitas,
deslocamentos e forças de superfície e, conseqüentemente os valores dos
deslocamentos dos pontos da interface MEC – MEF.
148
CAPÍTULO 8 APLICAÇÕES NUMÉRICAS
149
8.1 ANÁLISE ESTÁTICA
Nesta seção são apresentados alguns exemplos utilizando a formulação do
Método dos Elementos de Contorno para demonstrar a funcionalidade do programa
desenvolvido para o cálculo de deslocamentos e tensões.
Foi analisada a viga mostrada na
Figura 8-1, para a qual a teoria elementar da
flexão de vigas proporciona resultados com boas aproximações para os deslocamentos.
A viga tem seção retangular delgada e é engastada numa extremidade e livre na outra.
Os exemplos são adimensionais, tendo em vista apenas a finalidade de
comparação com soluções conhecidas. A largura t da viga foi assumida como unitária.
O Módulo de Elasticidade tem valor unitário e o Coeficiente de Poisson é nulo.
Figura 8-1 - Dimensões e vinculação da viga analisada.
8.1.1 EXEMPLO 1: VIGA SOLICITADA POR FORÇA NORMAL
UNIFORMEMENTE DISTRIBUÍDA CONSTANTE NA EXTREMIDADE
LIVRE
A discretização da viga da Figura 8-1 foi feita com vinte e oito nós, sendo nós
duplos nos cantos para melhor representar vinculações e cargas, conforme mostrado
na Figura 8-2. Na primeira análise os 28 nós foram agrupados em 24 elementos
lineares (Figura 8-3) e na seguinte em 12 elementos quadráticos (Figura 8-4).
1
t
=
3
h
=
12
L
=
150
Figura 8-2 – Viga com 28 nós no contorno.
Figura 8-3 - Viga com 28 nós no contorno em 24 elementos lineares.
Figura 8-4 - Viga com 28 nós no contorno em 12 elementos quadráticos.
• • •
10
9
• •
8
• •
7
• •
1
• • •
2
• •
3
• •
4
• •
•
½
°
•
¾
°
¿
•
½
°
•
¾
°
¿
•
5
6
•
°
•
®
°
¯
•
°
•
®
°
¯
•
11
12
x
y
9
12
21
22
x
y
20
• •
19
•
18
•
17
•
16
•
15
•
14
•
13
•
1
• •
2
•
3
•
4
•
5
•
6
•
7
•
8
•
}
•
•
}
•
}
•
}
•
10
11
{
•
•
{
•
{
•
{
•
24
23
• • • • • • • • •
10
15
x
y
16
17
18
19
20
21
22
23
• • • • • • • • •
1
2
3
4
5
6
7
8
9
•
•
•
•
•
11
12
13
14
•
•
•
•
•
28
27
26
25
24
151
A viga é solicitada por uma força de tração distribuída na extremidade direita
(P=30). A vinculação é feita por apoios móveis nos nós 24, 25 27 e 28 (restringem
deslocamento na direção x) e apoio fixo no nó 26 (restringe deslocamentos nas
direções x e y), como pode ser visto na Figura 8-5.
Figura 8-5 - Vinculação e carregamento da viga.
Os deslocamentos horizontais dos nós indicados na Tabela 8-2, obtidos com a
utilização do MEC, são idênticos aos valores calculados com a equação analítica que
define o deslocamento axial de uma viga engastada, dado por:
=
N L
L
E A
∆
(8.1)
onde:
L
∆
= deslocamento na direção x do ponto analisado;
N
= força normal resultante na seção transversal;
L
= comprimento da viga;
E
= Módulo de Elasticidade Longitudinal do material;
10
x
y
•
•
•
•
•
11
12
13
14
•
•
•
•
•
28
27
26
25
24
p
152
A
= área da seção transversal.
Os valores analíticos da tensão normal, para o caso da tração, são:
N
A
σ
=
(8.2)
Tabela 8-1- Deslocamentos e tensões nos nós da extremidade livre da viga em balanço com força
normal uniformemente distribuída ( elementos linear e quadrático)
Nó Coordenada x Coordenada y Deslocamento (x)
Tensão Normal
1 0,00 0,00 0,00 10,00
2 1,50 0,00 15,00 10,00
3 3,00 0,00 30,00 10,00
4 4,50 0,00 45,00 10,00
5 6,00 0,00 60,00 10,00
6 7,50 0,00 75,00 10,00
7 9,00 0,00 90,00 10,00
8 10,50 0,00 105,00 10,00
9 12,00 0,00 120,00 10,00
12 12,00 1,50 120,00 10,00
Embora o exemplo processado seja simples, os resultados obtidos mostram que
a formulação apresentada é adequada para sólidos elásticos solicitados à tração.
A aproximação descontínua adotada elimina os problemas surgidos na análise
de nós de canto. As integrações utilizadas, tanto a numérica quanto a analítica
mostram-se eficientes junto à formulação.
153
8.1.2 EXEMPLO 2: VIGA SOLICITADA POR FORÇA CORTANTE NA
EXTREMIDADE LIVRE
Foi feita uma análise da viga mostrada na Figura 8-1, com a mesma vinculação
de exemplo anterior, sujeita a um esforço aplicado na extremidade livre, na forma de
uma força cortante.
Para estudar a influência da discretização na resposta, duas
malhas foram
estudadas. A primeira delas com 28 nós em 24 elementos lineares (Figura 8-3) e 28 nós
em 12 elementos quadráticos (Figura 8-4). A segunda com 40 nós em 36 elementos
lineares Figura 8-7) e 40 nós em 18 elementos quadráticos (Figura 8-8).
Figura 8-6 - Viga com 40 nós no contorno.
Figura 8-7 - Viga com 40 nós no contorno em 36 elementos lineares.
x
y
N
1
• •
N
2
•
N
3
•
N
4
•
N
5
•
N
6
•
N
7
•
N
8
•
N
9
•
N
10
•
N
11
•
N
12
•
}
•
•
}
•
}
•
}
•
}
•
}
•
13
14
15
16
17
18
P
19
• •
P
20
•
P
21
•
P
22
•
P
23
•
P
24
•
P
25
•
P
26
•
P
27
•
P
28
•
P
29
•
P
30
•
{
•
•
{
•
•
{
•
•
{
•
•
{
•
•
{
•
•
31
32
33
34
35
36
10
15
x
y
16
17
18
19
20
21
22
23
1
2
3
4
5
6
7
8
9
11
12
13
14
28
27
26
24
• • • • • • • • • • • • •
•
•
•
•
•
•
•
•
•
•
•
•
•
•
29
30
31
32
33
34
35
36
37
38
39
40
25
• • • • • • • • • • • • •
154
Figura 8-8 - Viga com 40 nós no contorno em 18 elementos quadráticos.
Para melhor representar o problema real, o carregamento foi aplicado com
variação parabólica. Os valores da carga em cada nó foram calculados pela equação
da parábola com área de valor 30 (equação 8.3), que representa a carga total aplicada.
Figura 8-9 - Representação do carregamento parabólico.
A equação da parábola é
2
( )
p y ay by c
= + +
. As condições do problema em
estudo são:
·para
0 ( ) 0 0
y p y c
=
= ∴ =
·para
3 ( ) 0 3
y p y b a
=
= ∴ = −
x
y
1
• • •
2
• •
3
• •
4
• •
5
• •
6
• •
10
• • •
11
• • •
12
• • •
13
• • •
14
• • •
15
• • •
{
•
•
•
{
•
•
•
{
•
•
•
16
17
18
}
•
•
•
}
•
•
•
}
•
•
•
7
8
9
3
( )
p y
y
155
·área =30
( )
3
2
0
6,67
3 30
20
a
ay ay dy
b
= −
°
− = ∴
®
=
°
¯
³
Chega-se então à função que representa o carregamento:
2
( ) 6,67 20
p y y y
= − +
(8.3)
O resultado analítico é o da viga de Timoshenko, que considera também o efeito
da força cortante na equação da elástica. Com a origem do eixo
x
mostrada na Figura
8-10 e sendo c o fator de forma, a equação da elástica
( )
v x
é dada por:
( )
3
3
2 - 3 +
6
PL x x cP
v L x
E I L L GA
ª º
§ ·
= + −
« »
¨ ¸
© ¹
« »
¬ ¼
(8.4)
A flecha máxima ocorre na extremidade livre, como mostrado na Figura 8-10, e
é dada por
3
max
1,2
3
PL PL
v
E I GA
= +
(8.5)
Figura 8-10 - Elástica da viga
Aplicando-se a equação da linha elástica a partir da teoria de vigas, pode-se
calcular o deslocamento transversal da mesma ao longo do seu comprimento e
comparar com os valores encontrados para as discretizações estudadas neste trabalho.
p
max
v
x
L
156
Na Tabela 8-2 apresentam-se os valores dos deslocamentos verticais (direção
y
)
dos pontos ao longo do eixo
x
situados na parte inferior da viga, para a malha com 28
nós (Figura 8-2). A ordenada x se refere ao sistema de eixos da Figura 8-2. A análise
dos resultados pode ser feita através da Figura 8-11
.
O erro da diferença entre o valor do deslocamento analítico e o deslocamento
obtido pelo MEC, é calculado como:
Teórico-MEC
Erro = 100
Teórico
§ ·
¨ ¸
© ¹
(8.6)
Tabela 8-2 - Viga em balanço com força cortante na extremidade livre com 28 nós no contorno.
Deslocamentos
Ordenada
x
Analítico
Elemento linear (L) Erro (%)
Elemento quadrático (Q) Erro (%)
1,5 -208,5 -178,4235 14,4252
-216,917 4,036
3 -732 -660,8231 9,7236 -745,904 1,8995
4,5 -1525,5 -1385,054 9,2066 -1534,224 0,5719
6 -2544 -2314,8529 9,0074 -2531,633 0,4861
7,5 -3742,5 -3407,7036 8,9458 -3729,722 0,3414
9 -5076 -4622,9836 8,9247 -5054,755 0,4185
10,5 -6499,5 -5921,571 8,8919 -6460,615 0,5983
12 -7968 -7257,5196 8,9167 -7924,218 0,5495
157
-9000
-8000
-7000
-6000
-5000
-4000
-3000
-2000
-1000
0
0 1,5 3 4,5 6 7,5 9 10,5 12
Ordenada x
Deslocamentos
Figura 8-11 - Elástica da viga em balanço com força cortante na extremidade livre: comparação
dos resultados com elementos lineares e quadráticos (28 nós no contorno).
Na Tabela 8-3 apresentam-se os valores dos deslocamentos verticais (direção
y
)
dos pontos ao longo do eixo
x
situados na parte inferior da viga, para a malha com 40
nós no contorno (Figura 8-6). A ordenada x da Tabela 8-3 se refere ao sistema de eixos
da Figura 8-6. A análise dos resultados pode ser feita através da Figura 8-12.
158
Tabela 8-3 - Viga em balanço com força cortante na extremidade livre com 40 nós no contorno.
Deslocamentos
Ordenada
x
Analítico
Elemento linear (L) Erro (%)
Elemento quadrático (Q) Erro (%)
1 -101,7775
-109,1801 7,2733 -107,458 5,5813
2 -350,222 -354,7085 1,2810 -339,528 3,0535
3 -732 -716,33001 2,1407 -750,954 2,5893
4 -1233,778
-1198,6004 2,8512 -1215,425 1,4875
5 -1842,222
-1782,6293 3,2348 -1821,87 1,1048
6 -2544 -2456,2606 3,4489 -2522,354 0,8509
7 -3325,778
-3206,4857 3,5869 -3301,709 0,7237
8 -4174,222
-4020,6068 3,6801 -4148,786 0,6094
9 -5076 -4886,011 3,7429 -5047,582 0,5599
10 -6017,778
-5789,5485 3,7926 -5988,38 0,4885
11 -6986,222
-6719,3068 3,8206 -6953,025 0,4752
12 -7968 -7660,26001 3,8622 -7933,103 0,4380
159
Figura 8-12 - Elástica da viga em balanço com força cortante na extremidade livre: comparação
dos resultados com elementos lineares e quadráticos (40 nós no contorno).
Os valores das tensões normais nas duas discretizações são apresentados na
Tabela 8-4 junto aos valores analíticos, obtidos por:
M
ı y
I
=
(8.7)
-9000
-8000
-7000
-6000
-5000
-4000
-3000
-2000
-1000
0
0 1 2 3 4 5 6 7 8 9 10 11 12
Ordenada x
Deslocamentos
160
onde:
M
= momento fletor na seção em estudo;
I
= momento de inércia da seção transversal;
y
= abcissa do ponto analisado.
Tabela 8-4 - Viga em balanço com força cortante na extremidade livre.
Tensão Normal
x y
Analítico
Elemento Linear (L)
Erro (%)
Elemento Quadrático(Q) Erro (%)
3 0,75
90 80,4325 10,6306
89,27002 0,8111
6 0,75
60 53,6489 10,5852
59,4723 0,8795
9 0,75
30 26,8872 10,3760
29,7332 0,8893
Os valores das tensões de cisalhamento nas duas discretizações são
apresentados na Tabela 8-5 junto aos valores analíticos, obtidos por:
V S
b I
τ
=
(8.8)
onde:
V
= força cortante na seção em estudo;
I
= momento de inércia da seção transversal;
S
= momento estático da seção transversal;
b
= largura da seção transversal.
161
Tabela 8-5 - Viga em balanço com força cortante na extremidade livre.
Tensão de Cisalhamento
x y
Analítico
Elemento Linear (L)
Erro (%)
Elemento Quadrático(Q) Erro (%)
3 0,75
11,25 10,2804 8,6187 10,5592 6,1404
6 0,75
11,25 10,2482 8,9049 10,4732 6,9049
9 0,75
11,25 10,2687 8,7227 10,6329 5,4853
Os resultados obtidos permitem dizer que a formulação do Método dos
Elementos de Contorno apresentada é coerente com o problema analisado. Observa-se
ainda que, à medida que a discretização envolve mais pontos nodais, os valores
obtidos na análise se aproximam daqueles utilizando-se a teoria clássica de vigas.
8.1.3 EXEMPLO 3: VIGA SOLICITADA POR MOMENTO CONCENTRADO NA
EXTREMIDADE LIVRE
Foi feita uma análise da viga mostrada na Figura 8-1, com a mesma vinculação de
exemplo anterior, sujeita a um esforço aplicado na extremidade livre, na forma de um
momento concentrado de valor 15.
As discretizações utilizadas são as mesmas do exemplo anterior. A primeira
análise com 28 nós no contorno, elementos lineares (Figura 8-3) e elementos
quadráticos (Figura 8-4). A segunda análise com 40 nós no contorno, elementos
lineares (Figura 8-7) e elementos quadráticos (Figura 8-8).
O carregamento aplicado nos nós de modo a se obter o efeito do momento
concentrado na extremidade livre da estrutura está mostrado na Figura 8-13 (a) para a
viga discretizada com 28 nós; na Figura 8-13 (b) com 40 nós.
162
Figura 8-13 - Carregamento aplicado nos nós da extremidade livre.
Os deslocamentos analíticosforam calculados por:
2
2
1 2
2
ML x x
v
E I L L
ª º
§ ·
= − +
« »
¨ ¸
© ¹
« »
¬ ¼
(8.9)
onde:
M
= momento fletor na seção em estudo;
I
= momento de inércia da seção transversal;
L
= comprimento da viga.
10
10
10
10
5
5
6,67
6,67
3,37
3,37
(
)
a
(
)
b
3
163
Figura 8-14 - Elástica da viga
Na Tabela 8-6 apresentam-se os valores dos deslocamentos verticais (direção
y) dos pontos situados ao longo do eixo x, na parte inferior da viga a partir da
extremidade vinculada, para a discretização com 28 nós. A análise dos resultados pode
ser feita através da Figura 8-15.
Tabela 8-6 - Viga em balanço com momento concentrado na extremidade livre com 28 nós no
contorno.
Deslocamentos
Ordenada x
Analítico
Elemento linear (L) Erro (%)
Elemento quadrático (Q) Erro (%)
1,5 -7,5 -6,8748 8,3360 -7,5 0
3 -30 -27,6854 7,7153 -30 0
4,5 -67,5 -64,0668 5,0862 -67,5 0
6 -120 -115,0871 4,0941 -120 0
7,5 -187,5 -180,7772 3,5855 -187,5 0
9 -270 -261,1957 3,2609 -270 0
10,5 -367,5 -356,387 3,0239 -367,5 0
12 -480 -467,8207 2,5374 -480 0
x
M
2
max
2
ML
v
EI
=
L
164
-600
-500
-400
-300
-200
-100
0
0 1,5 3 4,5 6 7,5 9 10,5 12
Ordenada x
Deslocamentos
,,
Figura 8-15 - Elástica da viga em balanço com momento na extremidade livre: comparação dos
resultados com elementos lineares e quadráticos ( 28 nós no contorno).
Na Tabela 8-7 apresentam-se os valores dos deslocamentos verticais (direção y)
dos pontos situados ao longo do eixo x na parte inferior da viga a partir da extremidade
vinculada, para a discretização com 40 nós. A análise dos resultados pode ser feita
através da Figura 8-16.
165
Tabela 8-7 - Viga em balanço com momento concentrado na extremidade livre com 40 nós no
contorno.
Deslocamentos
Ordenada x
Analítico Elemento linear (L)
Erro (%)
Elemento quadrático (Q)
Erro (%)
2 3,3333 3,1786 4,6410 3,3333 0
3 13,3333 12,6252 5,3108 13,3333 0
4 30 29,0444 3,1853 30 0
5 53,3333 52,0846 2,3413 53,3333 0
6 83,3333 81,6391 2,0330 83,3333 0
7 120 117,8264 1,8113 120 0
8 163,33333
160,6156 1,6639 163,33333 0
9 213,33333
210,0161 1,5550 213,33333 0
10 270 266,0321 1,4696 270 0
11 333,33333
328,6769 1,3969 333,33333 0
12 403,333333
397,9541 1,3337 403,333333 0
166
-600
-500
-400
-300
-200
-100
0
0 1 2 3 4 5 6 7 8 9 10 11 12
Ordenada x
Deslocamentos
L
Q
Figura 8-16 - Elástica da viga em balanço com momento na extremidade livre: comparação dos
resultados com elementos lineares e quadráticos (40 nós no contorno).
As tensões normais que compõem a Tabela 8-8 foram comparadas à tensão
normal analítica (equação (8.7)), na qual o momento M é constante em toda a viga.
Pode-se notar que com a utilização de uma malha mais refinada, as respostas tornam-
se mais precisas, como era de se esperar.
167
Tabela 8-8- Viga em balanço com momento concentrado na extremidade livre.
Coordenada
Tensão Normal
x Analítico 28 nós 40 nós
12,00 10,00 9,77 9,91
11,00 10,00 9,92
10,00 10,00 9,93
9,00 10,00 9,46 9,93
8,00 10,00 9,92
7,00 10,00 9,91
6,00 10,00 9,14 9,90
5,00 10,00 9,89
4,00 10,00 9,88
3,00 10,00 9,04 9,87
2,00 10,00 9,86
1,00 10,00 9,82
0,00 10,00 9,25 9,81
8.2 CONDENSAÇÃO ESTÁTICA
A solução proposta neste trabalho para se utilizar os nós duplos em problemas
dinâmicos e poder resolver o sistema de equações é um rearranjo das matrizes
168
geradas. Esse rearranjo, chamado de condensação, elimina as linhas e colunas
dependentes, como mostrado na seção 4.1.5.
O objetivo dessa análise estática é comprovar a validade do algoritmo criado
para executar essa técnica, antes de aplicá-la à análise dinâmica.
Foi estudada a viga engastada numa extremidade e livre na outra, sujeita a um
esforço aplicado na extremidade livre, na forma de uma força cortante de variação
parabólica na seção 10.1.2. A discretização é a mostrada na Figura 8-17 com 18
elementos.
O objetivo é calcular os deslocamentos segundo o eixo
y
, dos pontos localizados
na extremidade inferior ao longo do eixo
x
e compará-los aos valores obtidos no
exemplo 2, o qual foi feito sem a manipulação das matrizes.
Figura 8-17 - Discretização com 18 elementos .
Os resultados estão apresentados na Tabela 8-9, na qual foram colocados
também os valores analíticos, obtidos com a equação (8.4).
x
y
1
• • •
2
• •
3
• •
4
• •
5
• •
6
• •
10
• • •
11
• • •
12
• • •
13
• • •
14
• • •
15
• • •
{
•
•
•
{
•
•
•
{
•
•
•
16
17
18
}
•
•
•
}
•
•
•
}
•
•
•
7
8
9
169
Tabela 8-9 - Viga em balanço com força cortante na extremidade livre com 18 elementos.
Deslocamentos
Ordenada x
Analítico Com Condensação Erro(%) Sem Condensação Erro(%)
1 -101,7775 -129,22 26,96 -127,458 25,23
2 -350,222 -376,45 7,49 -375,528 7,23
3 -732 -768,324 4,96 -770,954 5,32
4 -1233,778 -1210,275 1,90 -1215,425 1,49
5 -1842,222 -1813,93 1,54 -1821,87 1,10
6 -2544 -2512,72 1,23 -2522,354 0,85
7 -3325,778 -3298,41 0,82 -3301,709 0,72
8 -4174,222 -4120,6068 1,28 -4148,786 0,61
9 -5076 -5026,361 0,98 -5047,582 0,56
10 -6017,778 -5973,23 0,74 -5988,38 0,49
11 -6986,222 -6942,37 0,63 -6953,025 0,48
12 -7968 -7915,78 0,66 -7933,103 0,44
170
Figura 8-18 – Elástica da viga em balanço com força cortante na extremidade livre: comparação
dos resultados com condensação e sem condensação das matrizes.
Após a comparação dos resultados (Figura 8-18), verifica-se que a diferença dos
valores dos deslocamentos obtidos com a condensação de matrizes proposta é muito
pequena em relação aos deslocamentos obtidos da maneira usual do Método dos
Elementos de Contorno. Isso demonstra a eficiência da formulação e da análise
realizada.
Essa técnica será usada nos exemplos seguintes, para a análise dinâmica de
estruturas.
-8000
-7500
-7000
-6500
-6000
-5500
-5000
-4500
-4000
-3500
-3000
-2500
-2000
-1500
-1000
-500
0
0 1 2 3 4 5 6 7 8 9 10 11 12
Ordenada x
Deslocamentos
171
8.3 MEC – VIBRAÇÃO LIVRE DE UMA VIGA EM BALANÇO
A análise de vibrações em estruturas ocupa um grande espaço, pois abrange o
estudo do desempenho estrutural de pontes, edifícios, pisos industriais entre outros. As
características dinâmicas de um sistema material podem fornecer importantes critérios
de avaliação da sua integridade.
Para a viga da Figura 8-1, engastada em uma extremidade e livre na outra,
foram analisados os cinco primeiros modos de vibração. Com a finalidade de se
comparar os resultados obtidos com valores já determinados por outros autores, as
dimensões da viga são altura = 6 e comprimento = 24. O material que a constitui possui
Coeficiente de Poisson = 0,2 e a relação Módulo de Elasticidade Longitudinal por
Massa Específica = 10.
A viga foi modelada para a análise pelo MEC de quatro maneiras diferentes e,
para cada modelagem é mostrada a divisão do domínio em células.
1. MEC-dual com 6 elementos, MEC-integração direta com 10 células:
x
y
5
4
1
2
6
°
°
°
°
®
°
°
°
°
¯
3
½
°
°
°
°
¾
°
°
°
°
¿
1
2
3
4
5
6
7
8
9
10
11
172
Figura 8-19 - Discretização do item 1.
2. MEC-dual com 10 elementos, MEC-integração direta com 18 células:
Figura 8-20 - Discretização do item 2.
3. MEC-dual com 20 elementos, MEC-integração direta com 32 células:
2
7
12
3
4
5
6
1
8
9
10
11
13
x
y
9
1
10
°
°
°
°
®
°
°
°
°
¯
5
½
°
°
°
°
¾
°
°
°
°
¿
1
3
5
6
7
8
12
13
14
15
8
7
6
2
3
4
2
4
9
11
10
16
17
1
3
5
6
7
8
12
13
14
15
2
4
9
11
10
16
17
173
Figura 8-21 - Discretização do item 3.
4. MEC-dual com 30 elementos, MEC-integração direta com 72 células:
x
y
17
1
19
°
®
°
¯
10
½
°
¾
°
¿
1
5
11
10
12
13
17
21
22
24
28
18
16
15
14
13
12
2
3
4
6
5
7
8
11
9
½
°
¾
°
¿
20
°
®
°
¯
2
3
4
6
7
8
9
14
15
16
18
19
20
23
25
26
27
29
30
31
1
5
11
10
12
13
17
21
22
24
28
2
3
4
6
7
8
9
14
15
16
18
19
20
23
25
26
27
29
30
31
x
y
N
1
28
°
®
°
¯
13
½
¾
¿
11
32
12
13
17
28
1
5
2
3
4
6
7
8
9
14
15
16
29
30
31
10
33
34
N
10
N
11
N
12
N
4
N
5
N
7
N
8
N
9
N
2
N
3
N
6
14
½
¾
¿
15
½
¾
¿
21
P
16
18
19
20
P
17
P
18
22
24
23
25
26
27
P
26
P
19
P
20
P
21
P
22
P
23
P
24
P
25
P
27
29
°
®
°
¯
30
°
®
°
¯
35 36 37 38 39 40 41 42 43 4
4 45
46 47 48 49 50 51 52 53 54 5
5 56
174
Figura 8-22 - Discretização do item 4.
A Tabela 8-10 apresenta uma comparação entre os resultados obtidos para os
períodos dos modos de vibração pela dupla reciprocidade com a condensação das
matrizes e pela integração direta com a divisão do domínio em células, para as quatro
discretizações apresentadas acima.
Tabela 8-10 – Períodos dos Modos de Vibração: Teóricos, MEC-Integração Direta, MEC-
Reciprocidade Dual
Discretização
MEC Teórico 6,21 1,24 0,96 0,54 0,34
Células 10 células 5,887790
1,082900 0,943738 0,529268 0,283176 1
Dual 6 elementos 5,58741 1,07943 0,93625 0,51999 0,28093
Células 18 células 6,156750
1,183820 0,955947 0,535613 0,308355 2
Dual 10 elementos
6,06576 1,17594 0,949302 0,53189 0,306212
Células 32 células 6,197080
1,201870 0,957827 0,545227 0,327453 3
Dual 20 elementos
6,142649
1,195891 0,953169 0,544139 0,54939
Células 72 células 6,208037
1,220195 0,959711 0,555014 0,337734 4
Dual 30 elementos
6,18349 1,21534 0,95568 0,525177 0,33605
11
32
12
13
17
28
1
5
2
3
4
6
7
8
9
14
15
16
29
30
31
10
33
34
21
18
19
20
22
24
23
25
26
27
175
Figura 8-23 - Vibração livre.
0,2 0
0,4 0
0,6 0
0,8 0
1,0 0
1,2 0
Ite m 1 Iite m 2 Ite m 3 Item 4
D is c re tiz a çã o
Período
5,40
5,60
5,80
6,00
6,20
6,40
cé lula
te oria
dua l
1
5
4
3
2
1,30
176
O processo da integração direta é mais eficiente, uma vez que há uma melhor
aproximação das densidades do domínio pelo uso de células internas. Porém esse
método acarreta um volume de dados e um esforço computacional maior.
Nota-se que o refinamento da malha de elementos de contorno melhora os
resultados, mas uma precisão maior é obtida com a utilização de pontos no interior do
domínio, que representa melhor as propriedades dinâmicas do sistema.
Conclui-se que as freqüências naturais dadas pelo MEC utilizando a
reciprocidade dual apresentam características de convergência próximas às obtidas
com a utilização de células.
8.4 MEC – VIBRAÇÃO LIVRE: ARCO
O arco analisado foi dividido em 12 elementos no contorno e 7 pontos internos
(Figura 8-24). As dimensões geométricas e as constantes físicas são: raio interno
0
2,5
r =
, raio externo
1
5
r
=
,
4
10
E
ρ
= e
0,25
ν
=
(PARTRIDGE, BREBBIA e
WROBEL,1992).
Figura 8-24 - Discretização do Arco - MEC
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
0
r
1
r
30
31
177
Figura 8-25 - Discretização do Arco - MEF
Figura 8-26 - Freqüências Naturais para o Arco.
Nota-se uma boa aproximação entre os resultados obtidos com o Método dos
Elementos de Contorno utilizando a reciprocidade dual proposta nesse trabalho, com
aqueles calculados pelo Método dos Elementos Finitos (Figura 8-26).
7,181 7,292
MEC MEF
10,864 10,72
MEC MEF
16,89 16,49
MEC MEF
18,43 18,35
MEC MEF
178
Figura 8-27 – Comparação das freqüências naturais para o arco: MEC e MEF.
As diferenças máximas entre os valores das freqüências calculados através do
MEF (108 nós) e do MEC (20 elementos) são respectivamente 1,5%, 1,3%, 2,4% e
0,4% (Figura 8-27). Conclui-se que o método utilizado é eficiente, pois mesmo com
poucos nós, os valores foram muito próximos.
6
8
10
12
14
16
18
20
1 2 3 4
Frequências Naturais
MEF
MEC
179
8.5 – ANÁLISE TRANSIENTE DE PROPAGAÇÃO DE ONDAS
ELÁSTICAS: VIGA ENGASTADA SUBMETIDA A CARREGAMENTO
TRANSVERSAL SÚBITO EM SUA EXTREMIDADE
A propagação de ondas determina um movimento no sistema, chamado resposta
dinâmica. Essa resposta depende de muitos fatores, que se fundamentam nas
características físicas e geométricas do sistema e varia de acordo com a distribuição
espacial e temporal da excitação nele aplicada.
Quando as cargas aplicadas são localizadas e de curta duração, o trem de
ondas formado pode exigir um refinamento adequado do modelo para a sua precisa
representação.
Também a escolha do esquema de integração numérica para se obter a solução
é parte fundamental do problema.
Analisa-se, nesse exemplo, o movimento da viga engastada apresentada na
Figura 8-28, submetida a um carregamento transversal súbito em sua extremidade, cuja
intensidade e comportamento temporal são mostrados na Figura 8-29.
Figura 8-28 – Viga engastada submetida a carregamento súbito.
1
t m
=
2
h m
=
4
L m
=
180
Figura 8-29 – Intensidade e comportamento temporal do carregamento.
Foram empregados 24 elementos de contorno e 3 pontos internos, conforme
Figura 8-23.
Figura 8-30 – Viga com 24 elementos no contorno.
As propriedades do material são: Coeficiente de Poisson
0,25
ν
=
, Módulo de
Elasticidade
40.000
E MPa
=
, Velocidades
346,41 /
d
c m s
=
e
200 /
s
c m s
=
.
São apresentados na Tabela 8-11 os valores dos deslocamentos do ponto
central da extremidade livre ao longo do tempo, desde o instante inicial até t = 0,1
segundo, com intervalos de t = 0,004 s. São comparados os valores analíticos e os
calculados pelo MEF, com aqueles obtidos pelo MEC com a formulação apresentada.
9
12
21
22
x
y
20
• •
19
•
18
•
17
•
16
•
15
•
14
•
13
•
1
• •
2
•
3
•
4
•
5
•
6
•
8
•
8
•
}
•
•
}
•
}
•
}
•
10
11
{
•
•
{
•
{
•
{
•
24
23
•
•
•
1
0,02
( )
Força kN
( )
Tempo s
181
Para uma análise mais completa, na Figura 8-31 são mostrados os valores
dos deslocamentos do ponto central até o tempo de t = 0,16 segundos em intervalos
0,001
t s
∆ =
.
Tabela 8-11 – Deslocamento do ponto central da extremidade livre ao longo do tempo.
Deslocamentos
Tempo
Teórico M E C Erro(%) MEF Erro(%)
0,004 0,04860224 0,05226077
7,53
0,0514581
5,88
0,012 0,5041906 0,477965
5,20
0,5181267
2,76
0,02 1,089705 1,030833
5,40
1,079698
0,92
0,028 1,54238 1,451879
5,87
1,534321
0,52
0,036 1,785773 1,650723
7,56
1,780359
0,30
0,044 1,879435 1,8118
3,60
1,910198
1,64
0,052 2,153219 2,075624
3,60
2,133426
0,92
0,06 2,478661 2,369699
4,40
2,403785
3,02
0,068 2,366467 2,325224
1,74
2,350587
0,67
0,076 1,881481 1,910756
1,56
1,879705
0,09
0,084 1,106628 1,126939
1,84
1,153044
4,19
0,1 2,338198 2,240765
4,17
2,295071
1,84
182
Figura 8-31 - Deslocamento do ponto central da extremidade livre ao longo do tempo.
Nota-se a boa concordância entre os resultados no intervalo de tempo analisado.
O MEC leva vantagem pelo menor tempo computacional necessário para a obtenção da
solução.
8.6 MEC - ANÁLISE TRANSIENTE DE PROPAGAÇÃO DE ONDAS
ELÁSTICAS: VIGA ENGASTADA SUBMETIDA A CARREGAMENTO
LONGITUDINAL SÚBITO EM SUA EXTREMIDADE
Para melhor avaliar a consistência numérica da formulação do MEC foi estudado
o problema de propagação de ondas, onde um grande número de modos de vibração
está presente. Neste exemplo são estudados os resultados obtidos para a propagação
-3,00
-2,00
-1,00
0,00
1,00
2,00
3,00
0,000
0,010
0,020
0,030
0,040
0,050
0,060
0,070
0,080
0,090
0,100
0,110
0,120
0,130
0,140
0,150
0,160
Tempo (s)
Deslocamento (0,0001 m)
M E F
M E C
183
de ondas longitudinais elásticas em uma barra fixa em uma das extremidades e
submetida a um carregamento longitudinal súbito na outra extremidade. Na Figura 8-32
mostra-se a intensidade e comportamento do carregamento e as características
geométricas do problema.
Figura 8-32 - Barra com carregamento longitudinal e variação da carga com o tempo.
Considera-se uma malha composta por 18 elementos quadráticos com nós
duplos nos cantos. Para simular o comportamento da barra sob tração, utilizou-se o
coeficiente de Poisson nulo. O Módulo de Elasticidade do material é E = 100 e o
amortecimento viscoso é nulo (CARRER, 1991).
O tempo total do teste foi de 19,2 segundos, ou seja, foram utilizados 64
iterações com intervalo de tempo t = 0,3 s. Foi estudado o deslocamento de um ponto
central da extremidade livre.
0,03
( )
Força kN
( )
Tempo s
12
m
3
m
x
y
1
• • •
2
• •
3
• •
4
• •
5
• •
6
• •
10
• • •
11
• • •
12
• • •
13
• • •
14
• • •
15
• • •
{
•
•
•
{
•
•
•
{
•
•
•
16
17
18
}
•
•
•
}
•
•
•
}
•
•
•
7
8
9
184
Na Figura 8.33 são apresentados os resultados analíticos, pelo Método dos
Elementos Finitos e pelo Método dos Elementos de Contorno com a utilização do
esquema Newmark de integração no tempo.
Na Figura 8-34 os deslocamentos calculados pelo MEC correspondem a
utilização do esquema de Houbolt e para melhor visualização apresenta-se a Tabela
8-12 com resultados até o tempo de 4,8 segundos.
0
0,0005
0,001
0,0015
0,002
0,0025
0,003
0,0
1,2
2,4
3,6
4,8
6,0
7,2
8,4
9,6
10,8
12,0
13,2
14,4
15,6
16,8
18,0
19,2
Tempo(s)
Deslocamentos(m)
MEF
Analítico
MEC
Figura 8-33 – Deslocamento do ponto central da extremidade: Newmark.
185
0,0000
0,0005
0,0010
0,0015
0,0020
0,0025
0,0030
0,0
0,9
1,8
2,7
3,6
4,5
5,4
6,3
7,2
8,1
9,0
9,9
10,8
11,7
12,6
13,5
14,4
15,3
16,2
17,1
18,0
18,9
Tempo(s)
Deslocamento(m)
Analítico
MEF
MEC
Figura 8-34 - Deslocamento do ponto central da extremidade livre: Houbolt.
Observa-se que os resultados obtidos com o emprego do método de Newmark
nessa aplicação são piores do que os obtidos com Houbolt. Tal característica já havia
sido verificada anteriormente (LOEFLER,1988) e parece independer da formulação do
MEC adotada. No esquema Houbolt, a vantagem está relacionada com o
amortecimento fictício, que é equivalente ao truncamento da participação dos modos
superiores na resposta dinâmica.
A avaliação dos resultados mostrados na Tabela 8-12 permite concluir que o
emprego da integração no tempo de Houbolt é a técnica mais adequada para a solução
desse problema.
186
Tabela 8-12 – Deslocamento no ponto central da extremidade livre: Houbolt.
Tempo Analítico M E C Erro(%) M E F Erro(%)
0,3 0,0003 0,000247261 17,58 0,000237722 20,76
0,6 0,0006 0,000662447 10,41 0,000652353 8,73
0,9 0,0009 0,000930048 3,34 0,00091064 1,18
1,2 0,0012 0,001212471 1,04 0,00115472 3,77
1,5 0,0015 0,001552455 3,50 0,00151465 0,98
1,8 0,0018 0,00186892 3,83 0,001833172 1,84
2,1 0,0021 0,002190602 4,31 0,00208068 0,92
2,4 0,0024 0,002322511 3,23 0,002307017 3,87
2,7 0,0021 0,002153462 2,55 0,002270819 8,13
3,3 0,0018 0,001929816 7,21 0,00183661 2,03
3 0,0015 0,001611722 7,45 0,00138472 7,69
3,6 0,0012 0,001283456
6,95
0,00119336
0,55
3,9 0,0009 0,000997894
10,88
0,00098643
9,60
4,2 0,0006 0,000654431
9,07
0,00057701
3,83
4,5 0,0003 0,000289663
3,45
0,00023816
20,61
187
8.7 ACOPLAMENTO MEF - MEC – ANÁLISE TRANSIENTE:
CARREGAMENTO TRANSVERSAL
Esse exemplo tem como objetivo verificar a precisão dos procedimentos
implementados. Analisa-se uma viga engastada com seção retangular submetida a um
carregamento transversal súbito, conforme mostrado na Figura 8-35.
Figura 8-35 - Viga com carregamento transversal e variação da carga com o tempo.
A malha utilizada contém 34 nós no contorno e 3 elementos finitos, dispostos na
forma apresentada na Figura 8-36.
Figura 8-36 - Discretização da viga.
F
( )
força kN
( )
tempo s
9
12
m
3
m
11
32
12
13
17
28
1
5
2
3
4
6
7
8
9
14
15
16
29
30
31
10
33
34
21
18
19
20
22
24
23
25
26
27
1
2
3
188
Figura 8-37 – Deslocamento do ponto central da extremidade livre.
Como a finalidade dessa análise é a comprovação dos resultados, considera-se o
Módulo de Elasticidade e a Densidade com valores unitários. O tempo total analisado
foi de 600 s, com intervalo de integração temporal de
t
∆
= 18,18 s.
Na Tabela 8-13 mostram-se os deslocamentos do ponto central da extremidade
livre obtidos utilizando-se o acoplamento Método dos Elementos de Contorno – Método
dos Elementos Finitos, o Método dos Elementos de Contorno com células e o Método
dos Elementos Finitos. Observando e comparando os resultados apresentados na
Figura 8-37 pode-se dizer que a resposta obtida pela utilização da formulação acoplada
MEC-MEF está em ótima concordância com os outros métodos analisados.
-5000
-4500
-4000
-3500
-3000
-2500
-2000
-1500
-1000
-500
0
0,0
36,4
72,7
109,1
145,5
181,8
218,2
254,5
290,9
327,3
363,6
400,0
436,4
472,7
509,1
545,5
581,8
tempo (s)
deslocamento (0,0001 m)
MEF
MEC-
MEF
MEC
célula
189
Tabela 8-13 - Deslocamento no ponto central da extremidade livre
Tempo M E F MEC - MEF MEC Células
0,0 0,000 0,000 0,000
36,4 -715,995 -734,334 -723,030
54,5 -1320,743 -1334,984 -1315,246
90,9 -3003,270 -2911,820 -2956,201
127,3 -4217,362 -4332,633 -4264,991
163,6 -4533,255 -4646,802 -4573,651
200,0 -3487,464 -3861,886 -3600,579
236,4 -1977,583 -2254,659 -2018,226
272,7 -462,059 -810,671 -696,813
309,1 -128,443 -50,335 -88,244
345,5 -784,481 -604,898 -674,199
381,8 -2412,473 -1956,074 -2224,376
418,2 -3845,955 -3570,350 -3712,265
454,5 -4560,360 -4454,513 -4482,963
490,9 -3995,020 -4426,018 -4296,127
527,3 -2545,618 -3195,630 -2846,860
545,5 -1654,822 -2456,008 -2116,981
190
CAPÍTULO 9 CONCLUSÕES
191
O programa desenvolvido permitiu a análise de diversas categorias de
problemas, não só análise de problemas estáticos, como também vibrações livres de
meios finitos e casos transientes, através dos métodos dos Elementos de Contorno,
Elementos Finitos e do acoplamento entre eles.
Foram estudados dois processos para análise de problemas dinâmicos com o
Método dos Elementos de Contorno, utilizando o conceito de matriz de massa.
A formulação com Reciprocidade Dual se ajusta perfeitamente à filosofia do
MEC, pois reduz o problema dinâmico de um corpo à análise de seu contorno,
eliminando a necessidade da resolução das integrais de domínio. Isso resulta numa
implementação computacional mais simples, pois o volume de dados é muito menor.
Por utilizar a solução fundamental do problema estático, diminui o esforço
computacional na análise dos problemas transientes, além de permitir a determinação
direta das freqüências naturais de uma estrutura. Esse método propõe o uso de uma
série de funções particulares para a função deslocamento, funções estas dependentes
da geometria e do tempo.
A escolha de uma classe de funções a ser utilizada na transformação da integral
de domínio em integral de contorno é um dos pontos chave do processo da
reciprocidade dual. Desde a introdução dessa técnica, as funções radiais, escritas em
termo da distância Euclidiana entre dois pontos, têm sido as mais utilizadas, uma vez
que os resultados obtidos com essas funções levam à convergência das soluções.
Recentes pesquisas sugerem o estudo da utilização de funções chamadas TPS (thin
plate splines) e MQs (multiquadráticas). Estudos sobre a precisão e a eficiência das
respostas obtidas devem motivar pesquisadores a explorar algumas dessas
possibilidades.
Nesse trabalho, em estruturas com cantos ou angulosidades, foram empregados
nós duplos, o que exigiu o desenvolvimento de uma técnica para manipulação das
192
matrizes. A singularidade da matriz de transformação dos deslocamentos impedia sua
inversão, necessária para a resolução do sistema de equações. Além se superar essa
dificuldade a técnica de condensação das matrizes diminuiu o número de equações do
sistema gerado e portanto o tempo despendido na sua solução.
A resposta dinâmica para casos com carregamentos harmônicos ou periódicos
de baixa freqüência é mais precisa, pois esses casos excitam os modos de vibração
mais baixos. Para casos onde a distribuição do carregamento é localizada e sua
aplicação é súbita, a solução numérica é difícil, pois há a participação de todos os
modos de vibração. Apesar disso, os resultados obtidos nos exemplos apresentados
podem ser considerados bastante satisfatórios.
O esquema Houbolt de avanço no tempo foi o mais apropriado para a utilização
com elementos de contorno, embora boa parte das imprecisões possa ser atribuída à
questão da discretização temporal. A continuidade do trabalho com a reciprocidade dual
deve abordar pesquisa de técnicas alternativas de avanço no tempo, que possam vir a
ser mais adequadas à formulação do MEC.
Na análise com elementos de contorno foi também usado o processo da
integração direta, que aproxima os termos inerciais incógnitos da equação integral por
polinômios interpoladores parametrizados em células. Essa técnica leva a integração de
todo o domínio em forma discretizada. Embora fuja do propósito do MEC, que é
discretizar apenas o contorno, os resultados obtidos são bastante precisos, motivo pelo
qual foi aqui apresentado. Com a utilização dos nós duplos surgiu uma descontinuidade
no interior do corpo. Apesar de não impedir a convergência dos resultados para o valor
exato, procurou-se resolver esse problema colocando-se pontos internos próximos aos
nós e dividindo esse trecho do domínio num número maior de células.
O acoplamento dos dois métodos possibilitou o aproveitamento dos benefícios
de cada técnica. Foi usada a divisão do domínio em sub-regiões, sendo que na
interface comum são admitidas condições de compatibilidade de deslocamentos e de
193
equilíbrio de forças. Para a determinação das incógnitas, utiliza-se o sistema de
equações gerado pelo MEC.
O desenvolvimento e a aplicação da técnica de condensação de matrizes na
formulação com o processo da reciprocidade dual para representar o fenômeno da
propagação de ondas em meios contínuos foi bem sucedido, como demonstram as
análises apresentadas. O bom desempenho do MEC em problemas transientes
motivaram o acoplamento com o MEF, para se tirar vantagens das características dos
dois métodos.
194
REFERÊNCIAS BIBLIOGRÁFICAS
195
AGNANTIARIS, J. P., POLYZOS, D., BESKOS, D. E. Some studies on dual
reciprocity BEM for elastodynamics analysis.
Computational Mechanics
, Springer-
Verlag,v.17, 1996, p 270-277.
AHMAD, S. ASCE, A.M., BANERJEE, K.P
. Free vibration analysis by BEM using
particular integrals.
Journal of Engineering Mechanics, ASCE, 1986 ,112 , p.682-695.
ALMEIDA, F.P.A. Aplicação do acoplamento entre o MEC e o MEF para o estudo
da interação dinâmica elastoplástica entre o solo e estruturas. Tese (Doutorado) –
Escola de Engenharia de São Carlos, Universidade de São Paulo, 2003.
ALMEIDA, F.P. A; CODA, H.B. Análise dinâmica tridimensional de meios
contínuos via MEC aplicando matriz de massa: influência da densidade de células.. In:
IBERIAN LATIN-AMERICAN CONGRESS ON COMPUTATIONAL METHODS IN
ENGINEERING
, 22, 2001, Campinas, 2001.
ALMEIDA, F.P.A., CODA, H. B. & MESQUITA, E. Soil- structure interaction by
TDBEM-FEM coupling: a stable procedure. Iberian Latin-American Congress on
Computational Methods in Engineering, 2004.
ANSYS. User’s manual: procedure. s.L.p.: SAS IP, 1995. v.1. (Revision
________. User’s manual:
comands. s.L.p.: SAS IP, 1995. v.2. (Revision
________. User’s manual: elements. s.L.p.: SAS IP, 1995. v.3. (Revision 5.2).
________. User’s manual: theory. s.L.p.: SAS IP, 1995. v.4. (Revision 5.2).
ARAÚJO, F. C. Formulação Acoplada FEM/BEM no domínio tempo para a
resolução de problemas elastodinâmicos 3D.
Congresso Latino Americano sobre
Métodos Computacionais para Engenharia
, Minas Gerais, Brasil, 1994, p 621-629.
196
ASCHER, U.M., MATTHEIJ, R.M.M., RUSSELL, R.D. Numerical solution of
boundary value problems for ordinary differential equations
. Philadelphia: SIAM,
1995. 595p.
ASSAN, A.E.
Método dos elementos finitos: primeiros passos. Campinas:
Editora Unicamp, 1996. 298p.
BALAS, J., SLADEK, J., SLADEK, V. Stress analysis by boundary element
methods.
Amsterdam: Elsevier, 1989. 683p.
BANERJEE, P.K. The boundary elemet methods in engineering. 2ed.
London: McGraw-Hill, 1994. 496p.
________. KOBAYASHI, S.
Advanced dynamic analysis by boundary
element methods.
London: Elsevier, 1992. 410p.
________. AHMAD, S., MANOLIS, G.D. Advanced elastodynamic analysis in
boundary element methods. In: BESKOS, D.E., ed.
Boundary element method in
Mechanics. Amsterdam: North Holland, 1987. V. 3, p 257 – 284
________. BUTTERFIELD, R.
Boundary element methods in engineering
science.
London: McGraw-Hill, 1981
BARBIRATO, J. C. C., VENTURINI, W. S. Método dos elementos de contorno
com a reciprocidade dual para a análise transiente tridimensional da mecânica do
fraturamento.
Cadernos de Engenharia de Estruturas, São Carlos, v.7, n. 24, p.
63-88, 2005.
BATHE, K.J. Finite Element procedures in engineering analysis, Englewood
Cliffs Prentice–Hall, 1982 735p.
197
________. WILSON, E.L. Numerical methods in finite element analysis.
Englewood Cliffs: Prentice-Hall, 1976. 528p.
BECKER, A. A.,
The Boundary Element Method in Engineering: a Complete
Course, McGraw Hill, London, UK,1992.
BECKER, E.B., COREY, G.F., ODEN, T.J.
Finite elements: an introduction.
The University of Texas at Austin, Prentice-Hall Inc., Englewood Cliffs, New Jersey,
1981. 258p.
BEER, G. Implementation of combined boundary element – finite element
analysis with applications in geomechanics. In: BANEJEE, P.K., WATSON, I.O., ed.
Development in boundary elements. London: Applied. Science Publication, 1986.
v.4, p.11 – 225
BESKOS, D., ed.
Boundary element methods in structural analysis.
New
York: American Society of Civil Engineers, 1989. 341p.
________. Boundary element methods in dynamic analysis.
Applied.
Mechanics. Review
,.: v.40, p.1 - 23 ,1987
________. Boundary element methods in dynamic analysis: part II (1986-1996).
Applied Mechanics Reviews,
v.50, n.3, p.149 - 164, Mar.1997.
________. Boundary element methods in mechanics. Amsterdam: North-
Holland, 1987.v. 3, 3p.
BETTI, E. Teoria dell Elasticita.
Il Nuovo Cimento. Ser-3 ,1872. p 6-10.
BOYCE, W.E., DIPRIMA, R.C.
Equações diferenciais e problemas de valores
de contorno.
5.ed. Rio de Janeiro: Guanabara Dois, 1998. 387p.
198
BREBBIA, C: A The boundary element method for engineers. Petech Press ,
London, 1984, 189 p.
________ DOMINGUEZ, J.
Boundary elements na introductory course.
New York: McGraw-Hill, 1989. 293p.
________. FERRANTE, A.J. The finite element technique. Editora da
UFRGS, Universidade Federal do Rio Grande do Sul, 1975. 409p.
________. TELLES, J.C.F., WROBEL, L.C.
Boundary element techniques
:
theory and applications in engineering. Berlin: Springer-Verlag, 1984. 464p.
________. WALKER, S.
Boundary Element Techniques in Engineering,
London, Butterworth & Co , 1980, 210 p.
________. CHUANG, P. Applications of the boundary element method for
solving elastodynamics. In: CAMAK, A.S., ABDEL-GHAFFAR, A.M., BREBBIA, C.A.,
ed.
Problems in soil dynamics and earthquake engineering. Rotterdam: A.A.
Balbema1982, p 381 – 408
________. TELLES, J.C.F., WROBEL, L.C. Boundary element techniques:
theory and applications in engineering.
Berlin: Springer, 1984. 462p.
CARRER, J.A.M.
Técnicas implícitas para análise elastoplástica estática e
dinâmica com o método dos elementos de contorno.
Tese (Doutorado em
Engenharia Civil) – COPPE, Universidade Federal do Rio de Janeiro, 1991. 95p.
CAUCHY, A. L.
Exercices de Mathématique, ´Sur l’équilibre et le mouvement
d’un systemé de points matériels par des forces d’attraction ou de répulsion
mutuelle.1828. apud Love, A. E. H.
A treatise on the mathematical theory of
elasticity
. 4
th
ed, New York, Dover, 1944.
199
CHEN, G., ZHOU, J. Boundary element methods. London: Academic, 1992.
646p.
CHEN, C. S. , BREBBIA, C. A. , POWER, H. Dual Reciprocity Method using
compactly supported radial basis functions.
Communication in Numerical Methods in
Engineering
, v. 15, 1999, p 137-150.
CHENG, A.H.-D., LAFE, O. , GRILLI, S. Dual-reciprocity BEM based on global
interpolation functions.
Engineering Analysis with Boundary Elements,
Great
Britain
., v.13, n.4, p.303-311, April , 1994.
CHOI, C.K., AHN, J.S. A coupled formulation of finite and boundary element
methods in two-dimensional elasticity.
Computers and Structures, Great Britain., v.32,
n.6, p.1429-1437, Mês,
1989.
CODA, H.B.
Análise tridimensional transiente de estruturas pela
combinação entre o método dos elementos de contorno e o método dos
elementos finitos.
Tese (Doutorado em Engenharia Civil) - Escola de Engenharia de
São Paulo, 1993. 201p.
________.
Análise da vibração livre de meios elásticos bidimensionais pelo
método dos elementos de contorno.
Dissertação (Mestrado em Engenharia de
Estruturas) – Escola de Engenharia de São Carlos, Universidade de São Paulo, 1990.
130p.
________. VENTURINI, W.S. Alternative boundary element formulation for
elastodynamics. In:
Internacional Conference on Boundary Elements in
Engineering
,Sapporo, Japan, Sept. 1990. Proceedings. Southampton, CMP, Berlin,
v.1, n.12, p.517-534.
________. VENTURINI, W.S. Avaliação numérica de integrais singulares em
elementos de contorno para problemas elastostáticos e elastodinâmicos, In:
200
Congresso Ibero Americano sobre Métodos Computacionais para Engenharia, 15.,
1994, Belo Horizonte, p 641 – 650
________. VENTURINI, W.S. Formulação não singular do MEC para análise de
problemas dinâmicos no domínio tempo, In:
Congresso Ibero Americano sobre
Métodos Computacionais para Engenharia
, 15., 1994, Belo Horizonte, p 631 – 640
________. VENTURINI, W.S. Vibração livre de meios elásticos bidimensionais
pelo método dos elementos de contorno. In.
Congresso Ibero
Americano sobre
Métodos Computacionais para Engenharia ,
11., 1990, Rio de Janeiro, v.2, p.781-
790
__________ VENTURINI, W. S. Interação estática solo-estrutura através do
acoplamento entre o MEC e o MEF. In:
Simpósio Interação Estrutura-Solo em
Edifícios
, 1, 2000, São Carlos.
___________COOK, R.D.
Concepts and applications of finite element
analysis.
2.ed. New York: John Wiley, 1981. 403p.
________. MALKUS, D.S., PLESMA, M.E.
Concepts and applications of
finite element analysis.
3.ed. New York: John Wiley, 1989. 402p.
CRUSE, T.A. A direct formulation and numerical solution of the general transient
elastodynamic problem II ,
Journal of Mathematical Anallysis and Apllications, v.22,
1968. p 341-355
________.., RIZZO, F.J A direct formulation and numerical solution of the
general transient elastodynamic problem I ,
Journal of Mathematical Anallysis and
Apllications
, v.22, 1968.
p 244 - 259
DESAI, C.S. Elementary finite element method. Englewood Cliffs: Prentice-
Hall, 1979.
434p.
201
________. ABEL, J.F. Introduction to the finite element method: a numerical
method for engineering analysis. New York: Van Nostrand Reinhold, 1972. 477p.
DOMINGUEZ, J., ALARCON, E. Elastodynamics. In: BREBBIA, C.A., ed.
Progress in boundary element methods. New York Halsted Press , 1981. p 213 –
257
ERINGEN, A C., SUHUBI, E. Elastodynamics. New York : Academic, 1975
FENNER, R.T.
Finite element methods for engineers.
London: Macmillan,
1975.
171p.
FERRO, N. C. P. (1999).
Uma combinação MEC/MEF para análise de
interação solo-estrutura
.São Carlos. Tese (Doutorado) – Escola de Engenharia de
São Carlos, Universidade de São Paulo.
FONTES, J.A.
Implementação do método dos elementos de contorno para
elasticidade bidimensional com o uso de microcomputadores.
Dissertação
(Mestrado em Ciências) – COPPE, Universidade Federal do Rio de Janeiro, 1987.
169p.
FREDHOLM, I. Sur une classe d’es fonctionelles.
Acta Math. n 27,1903. p.365-
390
FRISWELL, M.I., MOTTERSHEAD, J.E.
Finite element model updating in
structural dynamics.
Dordrecht: Kluwer, 1995. 286p.
GIPSON, G.S.
Boundary element fundamentals: basic concepts and recent
developments in the Poisson equation. Southampton: Computational Mechanics
Publications, 1987. 287p. (Topics in Engineering, v.2)
202
GOLBERG, M.A. et al. Some comments on the use of radial basis functions in
the dual reciprocity method.
Computational Mechanics,
Springer Verlag., v.22, n.1,
July, 1998, p.61-69.
________. et al. Some comments on the use of radial basis functions in the dual
reciprocity method.
Computational Mechanics, Springer Verlag , v.21, n.2, Mar. 1998,
p.141-148.
________. CHEN,C. S., BOWMAN, H. Some recent results and proposals for the
use of radial basis functions in the BEM.
Engineering Analysis with Boundary
Elements
, Great Britain, v.23, 1999, p 285-296.
GOULD, P.L.
Introduction to linear elasticity. New York: Springer-Verlag,
1983. 157p.
GREEN,G.
An essay on the application of mathematical analysis to the
theories of eletricity and magnetism
. Nottingham, 1828.
HALL, W.S.
The boundary element method. Dordrecht: Kluwer, 1994. 227p.
HARTMANN, F.
Introduction to boundary elements: theory and applications.
Berlin: Springer, 1989. 416p.
________. OWEN, D.R.J.
Finite element programming.
Orlando: Academic,
1977.
305p.
HINTON, E.H., OWEN, D.R.J.
An introduction to finite element
computations.
Swansea: Pineridge, 1979. 385p.
HRENIKOFF, A. Solution of problems in elasticity by the framework method.
Journal Applied Mechanics
, v. A8, 1941, p 169-175.
203
HUEBNER, K.H., THORTON, E. A., BYROM, T.G. The finite element method
for engineers.
3.ed. New York: J.Wiley, 1995. 627p.
HUGHES, T.J.R. The finite element method: linear static and dynamic finite
element analysis Englewood Cliffs: Prentice Hall, 1987, 803 p.
KAMIYA, N., ANDOH, E., NOGAE, K. Eigenvalue analysis by the boundary
element method: new developments.
Engineering Analysis with Boundary
Elements
., v.12, n.3, p.151-162, 1993.
KIKUCHI, N. Finite element methods in mechanics. Cambridge: Cambridge
University, 1986. 418p.
KITAHARA, M.
Boundary integral equation methods in eigenvalue problems
of elastodynamics and thin plates.
Amsterdam: Elsevier, 1985. 281p.
LACERDA, L.A., VELOSO, M.J.G., MANSUR, W.J. O método dos elementos de
contorno para elastodinâmica bidimensional no domínio da freqüência com emprego de
elementos quadráticos. In:
Congresso Ibero Latino Americano sobre Métodos
Computacionais para Engenharia,
15., Belo Horizonte, 1994. Anais... Belo
Horizonte: CILANCE, 1994. P.601-610.
LEI, X.,QINGHUA, D. Recent engineering applications in China on coupling
boundary elements with finite elements.
Applied Mech. Review ,v. 50, n. 12, 1997, p
731-740.
LOEFFLER, C.F. Comparações entre métodos elastodinâmicos estabelecidos
com o método dos elementos de contorno. In:
Congresso Ibero Latino Americano
sobre Métodos Computacionais para Engenharia,
15., 1994, Belo Horizonte p 611 –
620
204
________. Uma formulação alternativa do método dos elementos de
contorno aplicada a problemas de campo escalar.
Tese (Doutorado em Ciências) –
COPPE, Universidade Federal do Rio de Janeiro, 1988. 156 p
.
________. MANSUR, W.J. Analysis of time integration schemes for
boundary element applications to transient wave propagation problems.
Boundary element techniques:
applications in stress analysis and heat transfer.
Computational Mechanics Publication, 1987.
LOVE, A.E.H. A treatise on the mathematical theory of elasticity. 4.ed. New
York: Dover, 1944. 309p.
MACKERLE, J., BREBBIA, C.A. The boundary element reference book.
Berlin: Springer, 1988. 382p.
MANOLIS, G.D., BESKOS, D.E.
Boundary element methods in
elastodynamics.
London: Unwin Hyman, 1988. 282p.
________. , BESKOS, D.E. Dynamic response of lined tunnel by an
isoparametric boundary element method.
Comp. Method Appl. Mech. Eng. v.36,1983.
p 291-307
MANSUR, W.J., BREBBIA, C.A. Formulation of the boundary element method
for transient problems governed by the scalar wave equations.
Applied Mathematics
Modeling
v. 6, 1982.
________. BREBBIA, C.A. Further development on the solution of the transient
scalar wave equation
.
________. ., BREBBIA, C.A. Boundary integral formulation of mass matrices for
dynamic analysis. In: BREBBIA, C.A., ed.
Topics in boundary element research.
Berlin: Springer Verlag, 1985. v.2, p.191 – 208.
205
________. ., BREBBIA, C.A Solution of parabolic and hyperbolic time
dependent problems using boundary elements.
Computational and Mathematics with
Applications,
1986 , 12B, p.1061-1072
________. BREBBIA, C.A. Transient boundary element elastodynamics using
the dual reciprocity method and modal superposition. In:
Boundary Elements
Conference
, 8., Tokyo, 1986, Berlin: Springer v.1., p.435-443.
________. BREBBIA, C.A., A new approach to free vibration analysis using
boundary elements.
Boundary Element Methods in Engineering, Berlin, September,
1982 p 313 – 326
MARTIN, H.C., GRAHAM F.C. Introduction to finite element analysis: theory
and apllication. New York: McGraw-Hill, 1973. 385 p.
McHENRY, D. A lattice analogy for the solution of plane stress problem .
Journal
Inst. Civ. Eng.
,v.21, 1943, p 59-82.
NARDINI, D. ,BREBBIA, C. A. A new approach to free vibration analysis using
boundary elements .In
. Boundary Elements methods in Engineering: Proceedings
of theourth International Seminar
, Southampton, England, September, 649 p,1982.
________.,BREBBIA, C. A. Boundary integral formulation of mass matrices for
dynamic analysis. In.BREBBIA, C. A. ed
. Topics in boundary element research-2 ,
Berlin, Springer- Verlag, 1985. p 191-208.
________.,BREBBIA, C. A. Transient boundary element elastodynamic using the
Dual Reciprocity Method and Modal Superposition. In
Boundary Element VIII
Conference,
Southampton, 1986, p435-443.
NAVIER.
Mém. Acad. Sciences, Paris, 7, 1827. apud Love, A. E. H. A treatise
on the mathematical theory of elasticity
. 4
th
ed,new York, Dover, 1944.
206
NEWMARK, N. M. Numerical methods of analysis in bars plates and elastic
bodies. In :
Grinter, L. E. ed . Numerical Methods in Analysis in Engineering,
Macmillan, 1949.
NIKU, S.M., ADEY, R.A. Computational aspect of the dual reciprocity method for
dynamics.
Engineering Analysis with Boundary Elements, Great Britain., v.18, n.1,
p.43-61, July, 1996.
________. KOBAYASHI, S., KITAHARA, M. Applications of the boundary
integral equation method to ergevalue problems of elasto dynamics. In:BREBBIA C. A
,ed..
Boundary element methods in engineering , Springer Verlag, New York, 1982,
p 291-311
________. KOBAYASHI, S., KITAHARA, M. Determination of ergevalues by
boundary element methods, developments. In: BANERJEE, P.K., SHAW, R.P., ed.
Developments in boundary element methods. London: Applied Science Publications,
1982. v.2., cap.6, p.143 – 176
________. KOBAYASHI, S., KITAHARA, M. Applications of the boundary
integral equation method to eigenvalue problems of elastodynamics
Boundary Element
Methods in Engineering
, Berlin, September, 1982 p 297 – 311
NOWAK, A.J., NEVES, A.C., ed.
The multiple reciprocity boundary element
method.
Southampton: Boston: C M P, 1994, 240p.
________. PARTRIDGE, P.W. Comparison of the dual reciprocity and the
multiple reciprocity methods.
Engineering Analysis with Boundary Elements, v.10,
n.4
, p.155-160, 1992.
ODEN, J. T. Finite elements: special problems in solid mechanics. Englewood
Cliffs: Prentice-Hall, 1984. 273p.
207
PALERMO Jr., L., Análise de peças de seção delgada como associação de
placas pelo método dos elementos de contorno
, Tese (Doutorado na área de
Engenharia de Estruturas), Escola de Engenharia de São Carlos, São Carlos, 1989,148
p.
PARTRIDGE, P.W.,BREBBIA,C.A .,WROBEL,L.C.
The dual reciprocity
boundary element method. London: Elsevier, 1992. 276p.
______. SENSALE, B. Hybrid approximation functions in the dual reciprocity
boundary element method.
Communications in Numerical Methods in Engineering.,
v.13, n.2, p.83-94, Feb. 1997
PAULA, F.A., TELLES, J.C.F., MANSUR, W.J. Combination of boundary
elements and finite elements to solve two-dimensional elasticity problems. In:
Boundary element techiniques:
application in stress analysis and heat transfer.
Computational. Mechanics. Publications, 1987. P 163 - 176.
RAMALHO, M.A. Sistema para análise de estruturas considerando interação
com o meio elástico.
Tese (Doutorado na área de Engenharia de Estruturas) - Escola
de Engenharia de São Carlos, Universidade de São Paulo, 1990. 389p.
RAO, S.S.
The finite element method in engineering.
2.ed. Oxford:
Pergamon, 1989.
643 p.
REDDY, J.N.
An introduction to the finite element method. New York:
McGraw-Hill, 1984. 495p.
RIBEIRO, D.B. (2005).
Análise da interação solo-estrutura via acoplamento
MEC- MEF
. São Carlos. 121p. Dissertação (Mestrado) – Escola de Engenharia de São
Carlos, Universidade de São Paulo.
208
RIZZO, F. J. An integral equation approach to boundary value problems of
classical elastostatics.
Quart. Appl. Math.,
n.25 1967. p 83-95
SAADA, A.S.
Elasticity theory and applications. Krieger Publishing Company,
2. Ed. , Malabar , Flórida , USA, 1993, 321 p.
SANTIAGO, J.A.F. Implementação do método dos elementos de contorno
para elasticidade bidimensional com o uso de microcomputadores.
Dissertação
(Mestrado em Ciências) - Rio de Janeiro, 169p.
SEGERLIND, L.J.
Applied finite element analysis. 2.ed. New York: J. Wiley,
1984, 427 p.
SIQUEIRA, E. F. N., CARRER, J.A.M., MANSUR, W. J. Propagação de ondas
elásticas em 2D: Acoplamento MEC/MEF,
Computacional Methods in Engineering,
São Paulo, Brasi, 1999, p 148.1-148.16.
SMITH, I.M., GRIFFITHS, D.V.
Programming the finite element method. 2.ed.
New York: John Wiley, 1988. 469p.
SOLIANI, R., BARRETO, F.M.S. Guia de Consulta do Mathematica, Faculdade
de Engenharia Civil, Universidade Estadual de Campinas, 1997, 63p.
SOUSA, E.A.C.
Acoplamento do método dos elementos finitos e elementos
de contorno para tratamento de problemas estacionários da elastodinâmica.
Dissertação (Mestrado em Engenharia Mecânica) - Faculdade de Engenharia Mecânica,
Universidade Estadual de Campinas1992. 91p.
STEIN, E., WENDLAND, E., WOLFANG, L., ed.
Finite element and boundary
element techniques from mathematical and engineering point of view.
Wien:
Springer, 1988. 333p.
209
STOKES, G. G. On the dynamical theory of diffraction. Transactions of the
Cambridge Philosophical Society,
v.9, 1849, p.1
SZABO, B.A., BABUSKA, I.
Finite element analysis. New York: John Wiley,
1991.
368p.
TANG, W.
Transforming domain into boundary integrals in BEM: a
generalized approach. Berlin: Springer, 1988. 209p.
TIMOSHENKO, S. ,YOUNG, D. H. ,WEAVER, W. JR.
Vibrations Problems in
Engineering,
5 ed., 1990, New York: John Wiley & Sons, 521p
________. GOODIER, J.N.
Teoria da elasticidade. 3.ed. Rio de janeiro:
Guanabara Dois, 1980. 545p.
TOTTENHAM, H., BREBBIA, C., ed. Finite element techniques in structural
mechanics.
Southampton: Stress Analysis Pub., 301p.
VALLABHAN, C.V.G. Coupling of BEM/FEM techonology in overview. In:
BREBBIA. C.A., VENTURINI, W.S., ed. Boundary element techniques: application in
stress analysis and heat transfer.
Computational. Mechanics. Publications, 1987
VETURINI, W.S. Boundary element method in geomechanics. Berlin:
Springer,Verlag 1983. 246p.
Von
ESTORFF, O., KAUSEL, E. Coupling of boundary and finite elements for
soil structure interaction problems.
Earthquake Engineering and Structural
Dynamics,
v.18, n. 7, p. 1065 - 1075, October, 1989.
________. PRABUCHI, M.J. Dynamic response in the time domain by coupled
boundary and finite elements
. Computational Mechanics.
Springer Verlag, 1990 v. 6,
p 35 – 46.
210
WARBURTON, G. B. The dynamical behaviour of structures, 2. ed. Oxford,
Pergamon Press, 1976.
WEN, P.H., ALIABADI, M.H., ROOKE, D.P. A new method for transformation of
domain integrals to boundary integrals in boundary element method.
Communications
in Numerical Methods in Engineering,
v.14, n.11, p.1055-1065, Nov. 1998.
ZHU, S., ZHANG, Y. Solving general field equations in infinite domains with dual
reciprocity boundary element method.
Engineering Analysis with Boundary
Elements,
Great Britain., v.12, n.4, p.241-250, October, 1993.
ZIENKIEWCZ, O.C. The finite element method in engineering science. 2.ed.
London: McGraw-Hill, 1971. 521p.
________. CHEUNG, Y. K.
The finite element method in structural and
continuum mechanics:
numerical solution of problems in structural and continuum
mechanics. London: McGraw-Hill, 1967. 274p.
________. KELLY, D.W., BESS, P. The coupling of the finite element method
and boundary solution.
International Journal for Numerical. Methods. in
Engineering
de acordo com a NBR 6023 da Associação Brasileira de Normas Técnicas
(ABNT) de 1989.
211
ANEXO 1: NOTAÇÃO INDICIAL CARTESIANA
212
Resumem-se aqui algumas regras básicas:
• Convenção de Soma de Einstein
: quando um índice aparece repetido em uma
grandeza ou uma expressão, isto significa que o índice deve assumir todos os
valores no seu intervalo de variação e as grandezas definidas somadas entre si. Por
exemplo, a expressão:
3
1
i ij j ij j
j
a x a x
ξ
=
= =
¦
(A.1)
contém, no lado direito, o índice j repetido.
Tomando-se os valores de i = 1, 2, 3 por vez , obtém-se três equações:
1 1 11 1 12 2 13 3
2 2 21 1 22 2 23 3
3 3 31 1 32 2 33 3
j j
j j
j j
a x a x a x a x
a x a x a x a x
a x a x a x a x
ξ
ξ
ξ
= = + +
= = + +
= = + +
(A.2)
De maneira análoga, o módulo do vetor
x
, cujas componentes são
1 2 3
, ,
x x x
é
dado por:
2 2 2
1 2 3
i i
x x x x x x
= + + =
(A.3)
•
Delta de Kronecker
: é um operador cuja definição é dada por:
1
0
ij
ij
se i j
se i j
δ
δ
= =
°
®
°
= ≠
¯
(A.4)
213
Tomando-se os valores i = 1, 2, 3 obtém-se nove valores, dados por:
11 22 33
12 21 13 31 23 32
1
0
δ δ δ
δ δ δ δ δ δ
= = =
= = = = = =
• Convenção para Derivadas Parciais
: a vírgula, acompanhada por índices,
representa a derivada parcial com relação à coordenada dada pelo índice:
,
2
,
i
i j
j
i
i j k
j k
f
f
x
x x
ϕ
ϕ
∂
=
∂
∂
=
∂ ∂
(A.5)
214
ANEXO 2: CÁLCULO DO JACOBIANO DA
TRANSFORMAÇÃO PARA O ELEMENTO LINEAR
215
A expressão (4.9) nos permite calcular o valor do Jacobiano:
2 2
1 2
x x
d
J
d
η η η
∂ ∂§ · § ·
Γ
= + =
¨ ¸ ¨ ¸
∂ ∂
© ¹ © ¹
(A.6)
Para o elemento linear, tem-se:
x
=
1 2
1 2
0 0
0 0
ϕ ϕ
ϕ ϕ
ª º
« »
¬ ¼
1
1
2
1
1
2
1
2
x
x
x
x
½
° °
° °
° °
° °
® ¾
° °
° °
° °
° °
¯ ¿
(A.7)
Logo:
1 1 2 2
1 1 1
x x x
ϕ ϕ
= +
(A.8)
com
1 2
1 1
(1- ) (1 )
2 2
e
ϕ η ϕ η
= = +
216
Substituindo os valores de
1 2
e
ϕ ϕ
, obtém-se as ordenadas
1 2
x e x
:
1 2
1 1 1
1 1 2 2
1 1 1 1 1
1 1
(1- ) (1 )
2 2
1 1 1 1
2 2 2 2
x x x
x x x x x
η η
η η
= + +
= − + +
(A.9)
A derivada
1
x
η
∂
∂
será dada por:
( )
1 2 2 1
1
1 1 1 1
1 1 1
2 2 2
x
x x x x
η
∂
= − + = −
∂
(A.10)
Do mesmo modo,
( )
1 2 2 1
2
2 2 2 2
1 1 1
2 2 2
x
x x x x
η
∂
= − + = −
∂
(A.11)
No triângulo retângulo da figura acima:
2 1
1 1
2 1
2 2
cos
x x l
x x l sen
θ
θ
− =
− =
(A.12)
Substituindo (A.12) em (A.10) e (A.11), chega-se.
1
cos
2
x
l
θ
η
∂
=
∂
e
2
2
x
l
sen
θ
η
∂
=
∂
(A.13)
217
O Jacobiano é dado por (A.6):
( ) ( )
( )
2 2
2 2
cos cos
4 4
l l
J l lsen sen
θ θ θ θ
= + = + =
(A.14)
218
ANEXO 3: INTEGRAÇÃO NUMÉRICA DE GAUSS
219
A expressão utilizada para integração numérica unidimensional de Gauss é:
1
1
1
( ) ( )
n
k k
k
I f x dx f x w
=
−
=
¦
³
(A.15)
onde:
k
x
é a coordenada do k-ésimo ponto de integração;
k
w
é o fator peso associado a cada ponto;
n o número total de pontos de integração.
Estão listados na tabela a seguir os valores para x
k
e w
k
:
n
±
±±
± x
k
w
k
4 0.
0.
861
339
136
981
311
043
594
584
053
856
0.
0.
347
652
854
145
845
154
137
862
454
546
6 0.
0.
0.
932
661
238
469
209
619
514
386
186
203
466
083
152
265
197
0.
0.
0.
171
360
467
324
761
913
492
573
934
379
048
572
170
139
691
8 0.
0.
0.
0.
960
796
525
183
289
666
532
434
856
477
409
642
497
413
916
495
536
627
329
650
0.
0.
0.
0.
101
222
313
362
228
381
706
683
536
034
645
783
290
453
877
378
376
374
887
362
10 0. 973 906 528 517 172 0. 066 671 344 308 688
220
0.
0.
0.
0.
865
679
433
148
063
409
395
874
366
568
394
338
688
299
129
981
985
024
247
631
0.
0.
0.
0.
149
219
269
295
451
086
266
524
349
362
719
224
150
515
309
714
581
982
996
753
12 0.
0.
0.
0.
0.
0.
981
904
769
587
367
125
560
117
902
317
831
233
634
256
674
954
498
408
246
370
194
286
998
511
719
475
305
617
180
469
0.
0.
0.
0.
0.
0.
047
106
160
203
233
249
175
939
078
167
492
147
336
325
328
426
536
045
386
995
543
723
538
813
512
318
346
066
355
403
16 0.
0.
0.
0.
0.
0.
0.
0.
989
944
865
755
617
458
281
095
400
575
631
404
876
016
603
012
934
023
202
408
244
777
550
509
991
073
387
355
402
657
779
837
649
232
831
003
643
227
258
637
0.
0.
0.
0.
0.
0.
0.
0.
027
062
095
124
149
169
182
189
142
253
158
628
595
156
603
450
459
523
511
971
988
519
415
610
411
938
682
255
816
395
044
455
754
647
492
533
576
002
923
068
20 0.
0.
0.
993
963
912
128
971
234
599
927
428
185
277
251
094
913
325
0.
0.
0.
017
040
062
614
601
672
007
429
048
139
800
334
152
387
109
221
0.
0.
0.
0.
0.
0.
0.
839
746
636
510
373
227
076
116
331
053
867
706
785
526
971
906
680
001
088
851
521
822
460
726
950
715
141
133
218
150
515
827
419
645
497
0.
0.
0.
0.
0.
0.
0.
083
101
118
131
142
149
152
276
930
194
688
096
172
753
741
119
531
638
109
986
387
576
817
961
449
318
472
130
705
240
518
176
382
603
726
222
ANEXO 4: OBTENÇÃO DA SOLUÇÃO DESLOCAMENTO
223
O problema estático hipotético em um domínio infinito a ser analisado é:
lim,
0
j
m li
f
σ δ
+ =
(A.16)
A solução em termos de deslocamentos deste problema será chamada
j
li
ϕ
.
Pode-se escrever:
, ,
( ) 0
j
li jj li ji li
G G f
ϕ λ ϕ δ
+ + + =
(A.17)
ou:
2
~ ~
~ ~
( ) ( . ) 0
j
G G f
ϕ λ ϕ
∇ + + ∇ ∇ + =
(A.18)
O operador diferencial
∇
, chamado divergente, é definido como:
1 2 3
1 2 3
i i i
x x x
∂ ∂ ∂
∇ = + +
∂ ∂ ∂
G G G
(A.19)
O operador diferencial escalar
2
∇
, chamado Laplace, é definido como:
2 2 2
2
2 2 2
1 2 3
x x x
∂ ∂ ∂
∇ = + +
∂ ∂ ∂
(A.20)
Pode-se escrever
ϕ
como função de um vetor auxiliar
A
, chamado vetor de
Galerkin:
2
1
= 2 (1 - ) A - ( . A)
2 G
ϕ ν
ª º
∇ ∇ ∇
¬ ¼
(A.21)
224
Substituindo-se
ϕ
da equação (A.21) em (A.17), resulta:
2 2
(1 - ) ( ) 0
j
A f
ν
∇ ∇ + =
(A.22)
Com a função escolhida
1
j
f r
= +
, a expressão (A.22) torna-se:
2 2
(1 - ) ( ) (1 ) 0
A r
ν
∇ ∇ + + =
(A.23)
A resolução de (A.23) no problema de estado plano de deformação é feita de
modo simples, com o operador Laplaciano
2
∇
escrito em coordenadas cilíndricas à
uma variável, ou seja:
( )
2
1
( )
d d
r
r dr dr
ª º
∇ =
« »
¬ ¼
(A.24)
Resolvendo (A.23):
2 2
1
( ) -
1- 1 -
r
A
ν ν
∇ ∇ = −
(A.25)
Aplicando o operador (A.24):
2
1 1
- -
1- 1 -
d d r
r A
r dr dr
ν ν
ª º
ª º
∇ =
« »
« »
¬ ¼
¬ ¼
(A.26)
~
1 1 1
- -
1- 1 -
d d d d r
r r A
r dr dr r dr dr
ν ν
½
ª º
§ ·
=
® ¾
¨ ¸
« »
© ¹
¬ ¼
¯ ¿
225
2
1
- -
1- 1 -
d d d d r r
r r A
dr dr r dr dr
ν ν
½
ª º
§ ·
=
® ¾
¨ ¸
« »
© ¹
¬ ¼
¯ ¿
Integrando:
2 3
1
1
- -
2(1- ) 3(1 - )
d d d r r
r r A C
dr r dr dr
ν ν
½
ª º
§ ·
= +
® ¾
¨ ¸
« »
© ¹
¬ ¼
¯ ¿
2
1
1
- -
2(1- ) 3(1 - )
C
d d d r r
r A
dr r dr dr r
ν ν
ª º
§ ·
= +
¨ ¸
« »
© ¹
¬ ¼
Integrando:
2 3
1 2
~
1
- - ln
4(1- ) 9(1 - )
d d r r
r A C r C
r dr dr
ν ν
§ ·
= + +
¨ ¸
© ¹
3 4
1 2
- - ln
4(1- ) 9(1 - )
d d r r
r A C r r C r
dr dr
ν ν
§ ·
= + +
¨ ¸
© ¹
Integrando:
4 5 2 2
1
2
2 3
- - ln -
16(1- ) 45(1 - ) 2 4
2
d r r r r
r A C r
dr
r
C C
ν ν
§ ·
= + +
¨ ¸
© ¹
+ +
226
3 4
1
3
2
- - ln -
16(1- ) 45(1 - ) 2 4
2
d r r r r
A C r
dr
C
r
C
r
ν ν
§ ·
= + +
¨ ¸
© ¹
+ +
Integrando uma última vez, chega-se ao vetor
A
:
4 5 2 2 2
1
2
2 3 4
1
- - ln - -
64(1- ) 225(1 - ) 2 2 4 8
ln
4
r r r r r
A C r
r
C C r C
ν ν
ª º
§ ·
= + +
« »
¨ ¸
© ¹
¬ ¼
+ + +
(A.27)
Com a matriz
A
conhecida, acha-se a solução da equação (A.27), utilizando a
definição dos operadores (A.19) e (A.20). Assim, obtém-se:
3
2
, , , ,
(1- 2 ) 10
3 -
(5 - 4 ) 30 (1 - ) 3
li l i li l i
r
r r r r r
G G
ν ν
ψ δ
ν ν
ª º
§ ·
= + −
¨ ¸
« »
© ¹
¬ ¼
(A.28)
227
ANEXO 5: ASPECTOS DA IMPLEMENTAÇÃO DO
SOFTWARE
228
Para a implementação e validação dos desenvolvimentos teóricos apresentados
no presente trabalho, foram elaborados programas computacionais, com base nas
formulações apresentadas. O sistema foi desenvolvido na linguagem FORTRAN, não
havendo a preocupação de explorar a capacidade dos equipamentos.
O programa desenvolvido possui a opção de se trabalhar com MEC, MEF ou
acoplamento dos dois métodos, como se fossem três programas independentes.
Na análise pelo Método dos Elementos de Contorno (MEC), utiliza-se a solução
fundamental de Kelvin e o elemento estudado é o linear, isto é,
u e p possuem variação
linear sobre os elementos. Podem ser feitas análises estáticas e dinâmicas (problemas
de vibrações livres e problemas transientes). No estudo dos problemas dinâmicos foram
implementados os métodos de integração numérica de Newmark e Houbolt.
Na formulação do Método dos Elementos Finitos (MEF) foram utilizados
elementos de barra, e também é possível realizar análises estáticas e dinâmicas.
Por fim, é feito o acoplamento dos dois métodos através da técnica das sub
regiões. A ligação é feita grau de liberdade a grau de liberdade. As matrizes geradas
nos programas do MEC e MEF são alocadas para o acoplamento entre os métodos.
No programa principal é feita a criação dos arquivos para armazenar as matrizes
das sub-regiões. Em seguida são lidos os dados para cada sub-região e também das
interfaces, se houver. Processa-se a sub-rotina de elementos de contorno e/ou a sub-
rotina de elementos finitos. Finalmente, montam-se as matrizes para cada sub-região,
já com a contribuição dos nós das interfaces calculadas pelo MEF acopladas às
matrizes do MEC. Resolve-se então o sistema de equações, que fornece a reposta no
domínio de elementos de contorno.
A seguir apresenta-se o fluxograma do programa principal.
229
Leitura de Dados das
Sub-Regiões e Interface
Não
Sim
Leitura dos Parâmetros
Temporais
Tipo de
Sub-Região
Montagem das Matrizes
das Sub-Regiões
Resolução do Sistema
Saída dos Resultados
Finito
Contorno
Sub-Rotina
MEF
Sub-Rotina
MEC
Início
Estático
Fim
Livros Grátis
( http://www.livrosgratis.com.br )
Milhares de Livros para Download:
Baixar livros de Administração
Baixar livros de Agronomia
Baixar livros de Arquitetura
Baixar livros de Artes
Baixar livros de Astronomia
Baixar livros de Biologia Geral
Baixar livros de Ciência da Computação
Baixar livros de Ciência da Informação
Baixar livros de Ciência Política
Baixar livros de Ciências da Saúde
Baixar livros de Comunicação
Baixar livros do Conselho Nacional de Educação - CNE
Baixar livros de Defesa civil
Baixar livros de Direito
Baixar livros de Direitos humanos
Baixar livros de Economia
Baixar livros de Economia Doméstica
Baixar livros de Educação
Baixar livros de Educação - Trânsito
Baixar livros de Educação Física
Baixar livros de Engenharia Aeroespacial
Baixar livros de Farmácia
Baixar livros de Filosofia
Baixar livros de Física
Baixar livros de Geociências
Baixar livros de Geografia
Baixar livros de História
Baixar livros de Línguas
Baixar livros de Literatura
Baixar livros de Literatura de Cordel
Baixar livros de Literatura Infantil
Baixar livros de Matemática
Baixar livros de Medicina
Baixar livros de Medicina Veterinária
Baixar livros de Meio Ambiente
Baixar livros de Meteorologia
Baixar Monografias e TCC
Baixar livros Multidisciplinar
Baixar livros de Música
Baixar livros de Psicologia
Baixar livros de Química
Baixar livros de Saúde Coletiva
Baixar livros de Serviço Social
Baixar livros de Sociologia
Baixar livros de Teologia
Baixar livros de Trabalho
Baixar livros de Turismo
