PMR3401-Apostila10

Prévia do material em texto

PMR3401 - Mecânica Computacional para Mecatrônica 
MONTAGEM DO SISTEMA GLOBAL 
 
Antes do sistema global ser montado, deve-se estabelecer um esquema de numeração global para 
especificar a topologia do sistema. A tabela abaixo define a conectividade dos elementos da malha da 
figura da página 55. Como trata-se de um caso unidimensional o esquema de numeração pode parecer 
trivial, mas para problemas de duas ou três dimensões é o único meio de se especificar quais nós 
pertence a quais elementos. 
 
TOPOLOGIA DO SISTEMA PARA O ESQUEMA DE SEGMENTAÇÃO POR ELEMENTOS FINITOS DO PROBLEMA DA PÁGINA 55 
 
 Numeração dos nós 
Elementos Local Global 
1 1 1 
 2 2 
2 1 2 
 2 3 
3 1 3 
 2 4 
4 1 4 
 2 5 
 
 Uma vez que a topologia é fixada, a equação do elemento pode ser escrita para cada elemento 
usando coordenadas globais. Pode-se então adicioná-las uma a uma para montar o sistema global. O 
processo é ilustrado abaixo: 
 
 
PMR3401 - Mecânica Computacional para Mecatrônica 
CONDIÇÕES DE CONTORNO 
 
Note que quando as equações são montadas as condições de contorno se cancelam. Dessa maneira, o 
resultado final de [F] nas equações acima tem condições de contorno presentes apenas no primeiro e no 
último nó. Como T1 e T5 são dados, essas condições naturais nas extremidades das barras, dT(x1)/dx e 
dT(x5)/dx, são incógnitas do problema. Portanto as equações podem ser reescritas como: 
 
 
 
SOLUÇÃO 
 
A solução para a equação acima é a seguinte: 
 
( ) ( ) 3475.25324575.17366 54321 −===== xdx
dTTTTx
dx
dT
 
 
PÓS PROCESSAMENTO 
 
Os resultados podem ser exibidos graficamente. A figura abaixo mostra a comparação dos resultados da 
solução do problema pelo MEF e da solução exata. Note que o MEF acompanha a tendência geral da 
solução exata e, nesse caso fornece uma solução exata nos nós. Entretanto há uma discrepância no 
interior de cada elemento devido à natureza linear das funções de forma. 
 
 
PMR3401 - Mecânica Computacional para Mecatrônica 
 
PROBLEMAS BIDIMENSIONAIS 
 
Embora o desenvolvimento matemático aumente consideravelmente, a extensão dos conceitos do MEF 
de uma para duas dimensões é imediata. Deve-se seguir os mesmos passos apresentados até agora. 
 
DISCRETIZAÇÃO 
 
Em duas dimensões utilizam-se basicamente dois tipos de elementos simples: triângulos e quadriláteros. 
Nossa análise será limitada a elementos triangulares do tipo apresentados na figura abaixo. 
 
 
 y
 x
 2
 3
 1
 
 
MATRIZ DO ELEMENTO 
 
Analogamente ao caso unidimensional, o próximo passo é desenvolver as equações para aproximar a 
solução dentro do elemento. Para um elemento triangular, a solução mais simples é uma função 
polinomial de 1a ordem: 
 
 u(x,y) = a0 + a1,1 x + a 1,2 y (32) 
 
sendo u(x,y) a variável de estado, ai,j, os coeficientes e x e y, as variáveis independentes. Os valores 
dessa função u nos nós do triângulo (x1,y1), (x2,y2) e (x3,y3) são: 
 
u1 = a0 + a1,1 x1 + a1,2 y1 
u2 = a0 + a1,1 x2 + a1,2 y2 (33) 
u3 = a0 + a1,1 x3 + a1,2 y3, 
 
ou em forma matricial: 
 










=




















3
2
1
2,1
1,1
0
33
22
11
1
1
1
u
u
u
a
a
a
yx
yx
yx
, 
 
cuja solução é 
 
PMR3401 - Mecânica Computacional para Mecatrônica 
( ) ( ) ( )[ ]
( ) ( ) ( )[ ]
( ) ( ) ( )[ ]1233122312,1
2131323211,1
1221331132233210
2
1
2
1
2
1
xxuxxuxxu
A
a
yyuyyuyyu
A
a
yxyxuyxyxuyxyxu
A
a
e
e
e
−+−+−=
−+−+−=
−+−+−=
 (34) 
 
sendo Ae a área do elemento triangular, dada por: 
 
( ) ( ) ( )[ ]1221311323322
1 yxyxyxyxyxyxAe −+−+−= (35) 
As equações (34) podem ser substituídas em (32). Após agrupamento de termos, chega-se a: 
 
 u = N1u1 + N2u2 + N3u3 , (36) 
sendo 
 
( ) ( ) ( )[ ]
( ) ( ) ( )[ ]
( ) ( ) ( )[ ]yxxxyyyxyx
A
N
yxxxyyyxyx
A
N
yxxxyyyxyx
A
N
e
e
e
122112213
311331132
233223321
2
1
2
1
2
1
−+−+−=
−+−+−=
−+−+−=
 
 
 Como no caso unidimensional, a equação (36) fornece uma maneira de se estimar valores de u 
internos ao elemento com base nos valores nodais. A figura abaixo mostra as funções de forma e as 
funções de interpolação. Note que a soma das funções de interpolação é sempre igual a 1. 
 
 Analogamente ao caso unidimensional, há vários métodos 
disponíveis para se desenvolver as equações nos elementos baseado na 
EDP e nas funções de aproximação. As equações resultantes são 
consideravelmente mais complicadas que a (26). Entretanto, como as 
funções de aproximação são normalmente polinômios de ordem 1, os 
termos da matriz elementar final consistirão de polinômios de ordem 
baixa e constantes. 
 
CONDIÇÕES DE CONTORNO E MONTAGEM DO SISTEMA GLOBAL 
 
A imposição das condições de contorno e a montagem dos sistema 
global também ficam mais complicadas quando o MEF é aplicado em 
problemas bi e tridimensionais. Entretanto, como com a construção da 
matriz do elemento, a dificuldade se concentra mais na mecânica do 
processo que na complexidade conceitual. Por exemplo, o 
estabelecimento da topologia dos sistema que era trivial para uma 
dimensão torna-se um passo importantíssimo em duas e três dimensões. 
Em particular, a escolha de um dado esquema de numeração irá definir a 
estrutura de banda da matriz global e portanto a eficiência com a qual 
ele poderá ser resolvido. A figura abaixo mostra um esquema de 
numeração que foi desenvolvido para a solução, pelo MEF, de um 
problema de uma placa aquecida. 
 
 
 
PMR3401 - Mecânica Computacional para Mecatrônica 
 
TOPOLOGIA DO SISTEMA PARA O ESQUEMA DE SEGMENTAÇÃO 
POR ELEMENTOS FINITOS DO PROBLEMA AO LADO 
 
Numeração dos nós
Elementos Local Global
1 1
2 21
3 7
1 1
2 72
3 6
1 2
2 33
3 8
etc.
 
 
APLICAÇÃO DO MÉTODO DE RESÍDUOS PONDERADOS E GALERKIN EM 
DUAS DIMENSÕES 
 
 Seja a equação de Poisson em duas dimensões, sendo T a função temperatura: 
 
0),(2 =+∇ yxfT 
 
Solução aproximada: T~ 
 
Substituindo na equação de Poisson resulta o resíduo: 
 
RyxfT =+∇ ),(~2 
 
Aplicação do M.R.P. e Galerkin (Wi = Ni): 
 
∫∫ =
eA
i dARN 0 , i =1, 2, ..., m, sendo m, o número total de nós. 
 
Então 
 ( ) 0),(~2 =+∇∫∫ dAyxfTN
eA
i , i =1, 2, m 
 ( ) ∫∫∫∫ −=∇
ee A
iA i
dAyxfNdATN ),(~2 
 
Para se obter a forma fraca, aplica-se o Teorema de Green (integração por partes vetorial): 
 
( ) ( ) ( )∫∫∫∫∫ ∇⋅∇−⋅∇=∇ S TCS dSvudnvudSvu l
r2
 
 
Então 
 
PMR3401 - Mecânica Computacional para Mecatrônica 
( ) ( ) ( )∫∫∫∫∫ ∇⋅∇−⋅∇=∇
eee A
T
iC iA i
dATNdnTNdATN ~~~2 lr 
 
O primeiro termo do lado direito corresponde a uma integral de linha sobre o contorno de cada 
elemento, sendo que: 
 
( )
n
T
nT
∂
∂
=⋅∇
~
~ r
 
 
( ) ( ) n
c
b
T
A
TnNnTTNT
j j
j
j
ej
jj
j
jj
rrr
⋅





=⋅∇=⋅∇⇒= ∑∑∑
===
3
1
3
1
3
1 2
1~~
 
 
 
Esse termo fornece a condição de contorno natural do problema (condição de Neumann, ou derivada 
normal). 
 
Do segundo termo resultará a matriz de rigidez (“stiffness matrix”) do elemento, que para elementos 
triangulares de primeira ordem possui solução analítica, como segue: 
 
[ ]ycxba
A
N iii
e
i ++= 2
1
, 
 
12321312213
31213231132
23132123321
,,
,,
,,
xxcyybyxyxa
xxcyybyxyxa
xxcyybyxyxa
−=−=−=
−=−=−=
−=−=−=
 
 






=∇
i
i
e
i
c
b
A
N
2
1
, ∑∑∑
===






=∇=∇⇒=
3
1
3
1
3
1 2
1~~
j j
j
j
ej
jj
j
jj cb
T
A
TNTTNT 
 
( ) ( ) [ ] [ ]∫∫∫∫ ∑∫∫






















+





+





⋅=






















⋅=∇⋅∇
=
eee A
ii
e
A j j
j
jii
e
A
T
i dA
c
b
T
c
b
T
c
b
Tcb
A
dA
c
b
Tcb
A
dATN
3
3
3
2
2
2
1
1
12
3
1
2 4
1
4
1~
 
O termo dentro do integrando não depende de x e y, podendo ser colocado para fora da integral, e 
eA
AdA
e
=∫∫ . Portanto a expressão acima fica: 
 
[ ] ( ) ( ) ( )[ ]333222111
3
3
3
2
2
2
1
1
1 4
1
4
1
ccbbTccbbTccbbT
Ac
b
T
c
b
T
c
b
Tcb
A iiiiiie
ii
e
+++++=













+





+





⋅ . 
 
 Para cada elemento triangular, i = 1, 2, 3, e tem-se três equações, como segue: 
 
( ) ( ) ( )[ ]
( ) ( ) ( )[ ]
( ) ( ) ( )[ ]







+++++
+++++
+++++
333332323213131
323232222212121
313132121211111
4
1
4
1
4
1
ccbbTccbbTccbbT
A
ccbbTccbbTccbbT
A
ccbbTccbbTccbbT
A
e
e
e
 
 
PMR3401 - Mecânica Computacional para Mecatrônica 
que escritas em forma matricial ficam: 
 










⋅










+
++
+++
3
2
1
3333
32322222
313121211111
.
4
1
T
T
T
ccbbsim
ccbbccbb
ccbbccbbccbb
Ae
, sendo 1, 2 e 3 os vértices de cada triângulo da malha 
de elementos finitos. 
 
O vetor de carregamento dado por: 
 
∫∫
eA
idANyxf ),( , 
 
admitindo que f(x,y) é constante dentro de cada elemento e vale f, fica: 
[ ] [ ]
eeee
e eeeee
GiGiieGieGiei
e
A A iA ii
e
A iii
e
A iA i
YcXbafAYcAXbAa
A
f
ydAcxdAbdAa
A
fdAycxba
A
fdANfdANyxf
++=++=
=



 ++=++== ∫∫ ∫∫∫∫∫∫∫∫∫∫
22
2
)(
2
),(
 
onde: 
3
;
3
321321 yyyYxxxX
ee GG
++
=
++
= 
Substituindo na equação acima e após alguma álgebra, temos que: 










=∫∫ =
1
1
1
3
),( 3,2,1 eA ii
fAdANyxf
e
. 
ou seja, o carregamento uniforme é distribuído eqüitativamente pelos nós do elemento no caso do 
elemento triangular. 
 
Somando-se o termo de contorno, que dá as condições de contorno naturais do problema, chega-se a: 
 
∫ ∑ ⋅






⋅





+










=
eC j j
j
j
e
i
e dn
c
b
T
A
NfA lr
3
12
1
1
1
1
3
 (37) 
 
 Como a normal nr aponta sempre para fora do elemento, para elementos internos ao domínio, 
com suas arestas não pertencendo ao contorno externo, ou seja, na interface entre dois elementos, o 
segundo termo acima se anula. Ele permanecerá apenas nas arestas que estão sobre o contorno externo 
do domínio. Se o contorno externo C do domínio é isolado termicamente do meio externo, ou seja, não 
há troca de calor com o exterior, então a integral acima é nula, e sobra apenas o primeiro termo. 
 
 
 
 
 
 
 
 
 
 
 
 
PMR3401 - Mecânica Computacional para Mecatrônica 
 
 
MONTAGEM DA MATRIZ GLOBAL 
 
A montagem da matriz global é realizada de forma expedita a partir das matrizes de cada 
elemento finito, de maneira análoga ao caso unidimensional. 
 Para cada um dos elementos associamos duas numerações para seus vértices: uma global, 
resultante da numeração seqüencial de todos os vértices do domínio, e outra local, de 1 a 3, seguindo o 
sentido horário, como mostrado na figura abaixo. 
 
 y
 x
 2
 3
 1
 r
 s
 t
 
 
 Duas matrizes serão construídas: a matriz de rigidez e o vetor de carregamento locais do 
elemento, como descrito no item anterior. A matriz do elemento é transportada para a matriz global de 
dimensão m×m (sendo m o número total de nós do domínio) seguindo a relação de correspondência 
entre as numerações local e global, isto é, o elemento (1,1) da matriz local é somado ao elemento (r,r) 
da matriz global; o elemento (1,2) é somado ao elemento (r,s), e assim sucessivamente. A matriz 
resultante, assim como a matriz local, resulta simétrica. O exemplo a seguir ilustra esse processo. 
 
 1
 2
 3
 4
 1
 2
 3
 4
 5
 6
 
 
Numeração local 1 2 3 
Elemento 1 1 2 4 
Elemento 2 1 4 3 
Elemento 3 3 4 6 
Elemento 4 3 6 5 
 
 
 
 
 
PMR3401 - Mecânica Computacional para Mecatrônica 
 
 
 A figura a seguir mostra as configurações da matriz global e do vetor de carregamento para o 
problema acima. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
PMR3401 - Mecânica Computacional para Mecatrônica 
 
 
INTRODUÇÃO DAS CONDIÇÕES DE CONTORNO DE DIRICHLET 
 
Após a montagem do sistema global de equações algébricas torna-se necessária a imposição das 
condições de contorno no problema do tipo Dirichlet, ou seja, valores conhecidos de T em alguma 
porção do contorno do domínio, o que não foi realizado até então. 
 
 Para impormos tal condição devemos efetuar algumas alterações no sistema original de modo 
que, após a solução desse sistema, resulte nos vértices de triângulos situados nas fronteiras o valor de T 
imposto. 
 
 Seja p um nó cujo valor de T é conhecido e vale T . Para que após a solução do sistema 
tenhamos Tp = T basta fazermos na matriz global e no vetor de ações as seguintes alterações. Sendo Kij 
um termo genérico da matriz de rigidez e Fj um termo genérico do vetor de carregamento, faz-se: 
 
Kpp = 1, Kpj = 0 e Fp = T , para j ≠ p 
 
 
 O procedimento anterior elimina a simetria da matriz global, o que não é interessante do ponto 
de vista da solução do sistema. Para se recuperar a simetria basta fazer: 
 
Fjnovo = Fjantigo − Kjp .T , 
Kjp = 0, para j ≠ p 
 
 Como exemplo, suponhamos que o nó 4 do exemplo anterior tenha um valor de T conhecido 
igual a T . O sistema de equações após a introdução das condições de contorno será ( T=Φ no 
esquema): 
T1
T2
T3
T4
T5
T6
 
 
 
T