Introducao ao Metodo de Elementos Finitos - Prof Samuel da Silva (UNIOESTE-PR)

Prévia do material em texto

Universidade Estadual do Oeste do Paraná
UNIOESTE/Campus de Foz do Iguaçu
Centro de Engenharias e Ciências Exatas - CECE
Introdução ao Método
dos Elementos Finitos
Notas de Aulas
Prof. Dr. Samuel da Silva
Foz do Iguaçu, 2009.
Prefácio
Este texto apresenta as notas de aulas da disciplina optativa Introdução ao
Método dos Elementos Finitos (FEM) do curso de graduação em Engenharia
Mecânica do Centro de Engenharias e Ciências Exatas da Universidade Es-
tadual do Oeste do Paraná, Campus de Foz do Iguaçu. Inúmeros excelentes
livros textos são disponíveis sobre FEM (como por exemplo [1], [2] ou [6]),
porém a biblioteca da UNIOESTE é muito pobre no tema, não contendo pra-
ticamente nenhuma referência sobre o assunto. Sendo assim, eu espero que
esta apostila diminua esta deficiência, porém, sem a pretensão de substituir
um livro texto. Tudo que os alunos vão encontrar nesta apostila se refere a
temas e assuntos clássicos imensamente divulgados e disponíveis. Tentei ao
máximo transcrever o tema de maneira natural com rigor matemático, mas
não de forma exaustiva e focando sempre aplicações práticas. Neste sentido,
ao longo do texto inúmeros exemplos são resolvidos de forma computacional
para ilustração e fixação dos conceitos básicos de FEM. Por fim, o objetivo
é servir como uma referência básica para guiar o estudo de FEM em um
nível introdutório e acessível aos alunos da graduação e com boa qualidade
gráfica1. Espero contar com o apoio dos alunos e demais colaboradores para
melhorar este texto constantemente, visto que ele foi feito de forma rápida
e nesta primeira versão pode estar sujeito a falhas, sendo assim, sugestões,
correções e comentários são muito bem vindos. Boa leitura e estudo!
Prof. Dr. Samuel da Silva
fevereiro de 2009.
1O texto foi redigido com o LATEX2ε.
2
Sumário
Lista de Figuras 6
Lista de Tabelas 9
1 Introdução 10
1.1 Exemplos de aplicação . . . . . . . . . . . . . . . . . . . . . . 11
1.2 Etapas na solução de um problema via FEM . . . . . . . . . . 13
1.3 Discretização por elementos finitos . . . . . . . . . . . . . . . 13
2 Fundamentos Matemáticos Básicos de FEM 17
2.1 Análise vetorial . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.1.1 Produto escalar . . . . . . . . . . . . . . . . . . . . . . 18
2.1.2 Produto vetorial . . . . . . . . . . . . . . . . . . . . . 18
2.1.3 Gradiente . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.1.4 Divergente . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.1.5 Rotacional . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.1.6 Teorema da divergência . . . . . . . . . . . . . . . . . 20
2.1.7 Teorema de Green-Gauss . . . . . . . . . . . . . . . . . 20
2.2 Análise matricial . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.3 Equações diferenciais . . . . . . . . . . . . . . . . . . . . . . . 23
2.4 Tensores cartesianos . . . . . . . . . . . . . . . . . . . . . . . 24
2.5 Exercícios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3 Elementos Finitos Unidimensionais 26
3.1 Solução exata de problemas de deformação axial em barra uni-
forme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.2 Aproximação via método de Galerkin . . . . . . . . . . . . . . 29
3.2.1 Aplicação do método de Galerkin na solução do pro-
blema da barra axial . . . . . . . . . . . . . . . . . . . 31
3.2.2 Solução linear . . . . . . . . . . . . . . . . . . . . . . . 32
3.2.3 Solução quadrática . . . . . . . . . . . . . . . . . . . . 33
3
3.3 Forma de elementos finitos de soluções assumidas . . . . . . . 34
3.3.1 Funções de interpolação linear para problemas de se-
gunda ordem . . . . . . . . . . . . . . . . . . . . . . . 35
3.3.2 Funções de ponderação de Galerkin na forma de ele-
mentos finitos . . . . . . . . . . . . . . . . . . . . . . . 37
3.3.3 Funções de interpolação Hermitianas para um ele-
mento com dois nós . . . . . . . . . . . . . . . . . . . . 37
3.4 Solução de elementos finitos de problemas de deformação axial 39
3.4.1 Solução linear assumida . . . . . . . . . . . . . . . . . 39
3.4.2 Equações de elementos usando o método de Galerkin . 40
3.5 Exemplos de problemas unidimensionais em engenharia . . . . 44
3.5.1 Transferência de calor . . . . . . . . . . . . . . . . . . 45
3.5.2 Fluxo de potência . . . . . . . . . . . . . . . . . . . . . 46
3.5.3 Transferência de massa . . . . . . . . . . . . . . . . . . 46
3.5.4 Eletricidade . . . . . . . . . . . . . . . . . . . . . . . . 47
3.6 Função no Matlab para resolver um PVC unidimensional . . . 48
3.7 Exercícios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4 Elementos de Treliça, Viga e Frame 55
4.1 Treliças planas . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.1.1 Exemplo de implementação prática . . . . . . . . . . . 60
4.2 Deformação transversal em vigas . . . . . . . . . . . . . . . . 71
4.2.1 Equação diferencial para deflexão em vigas . . . . . . . 72
4.2.2 Condições de contorno em vigas . . . . . . . . . . . . . 76
4.2.3 Exemplo de solução exata em vigas . . . . . . . . . . . 78
4.2.4 Diagramas de cortante e momento usando o método
das seções . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.2.5 Elemento de viga com dois nós . . . . . . . . . . . . . 82
4.2.6 Solução cúbica assumida . . . . . . . . . . . . . . . . . 83
4.2.7 Matriz de rigidez local em vigas . . . . . . . . . . . . . 84
4.2.8 Exemplo de aplicação com carregamento concentrado . 87
4.2.9 Exemplo de aplicação com carregamento distribuído . . 98
4.2.10 Viga com carregamento trapezoidal . . . . . . . . . . . 105
4.3 Frames . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
4.3.1 Exemplo de aplicação com carregamento concentrado . 113
4.3.2 Exemplo de aplicação com carregamento distribuído . . 122
4.4 Considerações finais . . . . . . . . . . . . . . . . . . . . . . . . 129
4.5 Exercícios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
4
5 Análise Dinâmica de Estruturas via FEM 142
5.1 Matriz de massa para elementos estruturais . . . . . . . . . . 142
5.2 Matriz de massa para barra . . . . . . . . . . . . . . . . . . . 144
5.3 Matriz de massa para elemento de treliça plana . . . . . . . . 145
5.4 Matriz de massa para elemento de viga . . . . . . . . . . . . . 146
5.5 Matriz de massa para elemento de frame . . . . . . . . . . . . 147
5.6 Análise modal analítica - resposta livre . . . . . . . . . . . . . 147
5.7 Algoritmo de Newmark . . . . . . . . . . . . . . . . . . . . . . 151
5.8 Exemplo de aplicação . . . . . . . . . . . . . . . . . . . . . . . 156
6 Elementos Finitos Bidimensionais 163
6.1 Integração por partes em duas dimensões - Teorema de Green-
Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 165
6.1.1 Teorema da divergência de Gauss . . . . . . . . . . . . 165
6.1.2 Teorema de Green-Gauss . . . . . . . . . . . . . . . . . 166
6.2 Equações de elementos finitos usando o método de Galerkin . 167
6.3 Elementos finitos retangulares . . . . . . . . . . . . . . . . . . 172
6.4 Exemplo: equação de Laplace em um domínio quadrado . . . 178
6.5 Elementos finitos triangulares . . . . . . . . . . . . . . . . . . 179
7 Introdução à Mecânica dos Sólidos Computacional 184
7.1 Revisão básica de elasticidade . . . . . . . . . . . . . . . . . . 185
7.1.1 O conceito de tensão . . . . . . . . . . . . . . . . . . . 185
7.1.2 Tensão principal e direção principal . . . . . . . . . . . 187
7.1.3 Critérios de falhas . . . . . . . . . . . . . . . . . . . . 188
7.1.4 Relações constitutivas . . . . . . . . . . . . . . . . . . 190
7.2 Equação de elementos finitos em sólidos elásticos . . . . . . . 192
7.2.1 Problema no estado plano de tensões . . . . . . . . . . 195
7.2.2 Elemento triangular com três nós . . . . . . . . . . . . 195
Referências Bibliográficas 196
5
Lista de Figuras
1.1 Modelo FEM de um blocode motor [4]. . . . . . . . . . . . . . 14
1.2 Tipos e formas de elementos [2]. . . . . . . . . . . . . . . . . . 15
1.3 Diferença entre o contorno físico e a geometria do contorno do
modelo de elementos finitos [2]. . . . . . . . . . . . . . . . . . 16
1.4 Malhas válida e inválida para elementos com 4 nós [2]. . . . . 16
3.1 Barra uniforme carregada axialmente [2]. . . . . . . . . . . . . 28
3.2 Elemento simples com dois nós para um problema de segunda
ordem [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36
3.3 Elemento simples com dois nós para um problema de quarta
ordem [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.4 Elemento linear de barra para o problema de deformação axial
[2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.5 Barra biengastada discretizada com três elementos [3]. . . . . 44
3.6 Barra não-uniforme carregada axialmente [2]. . . . . . . . . . . 49
3.7 Barra uniforme carregada axialmente [2]. . . . . . . . . . . . . 50
3.8 Coluna de um prédio modelada com quatro elementos lineares
[2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.9 Placa com fonte de calor [2]. . . . . . . . . . . . . . . . . . . . 52
3.10 Transferência de calor através de uma aleta [2]. . . . . . . . . 53
3.11 Perfil de velocidades de um fluído viscoso escoando entre duas
placas [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.1 Coordenadas global e local para uma barra axial orientada no
plano [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.2 Treliça plana a ser analisada [4]. . . . . . . . . . . . . . . . . . 60
4.3 Visualização da malha implementada com o Matlab R©. . . . . 63
4.4 Flexão em vigas [2]. . . . . . . . . . . . . . . . . . . . . . . . . 73
4.5 Forças e momentos agindo em um elemento diferencial da viga
[2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
4.6 Condições de contorno comuns em vigas [2]. . . . . . . . . . . 76
4.7 Viga pinada e apoiada em uma mola. . . . . . . . . . . . . . . 78
6
4.8 Soluções exatas para a viga com várias constantes. . . . . . . . 80
4.9 Exemplo de viga com carregamento. . . . . . . . . . . . . . . . 81
4.10 Elemento de viga com flexão. . . . . . . . . . . . . . . . . . . 82
4.11 Carregamento distribuído concentrado nos nós. . . . . . . . . . 85
4.12 Viga com carregamento concentrado aplicado. . . . . . . . . . 88
4.13 Deslocamento transversal da viga solucionado via FEM com
funções de interpolação Hermitiana. . . . . . . . . . . . . . . . 96
4.14 Momento fletor M(x) calculado usando FEM. . . . . . . . . . 97
4.15 Cortante V (x) calculada usando FEM. . . . . . . . . . . . . . 97
4.16 Tensão normal σ(x) devido a flexão calculada usando FEM. . 98
4.17 Viga não-uniforme com carregamento distribuido. . . . . . . . 99
4.18 Deslocamento transversal da viga com carregamento distri-
buído solucionado via FEM com funções de interpolação Her-
mitiana. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
4.19 Momento fletor M(x) na viga com carregamento distribuído
calculado usando FEM. . . . . . . . . . . . . . . . . . . . . . . 106
4.20 Cortante V (x) na viga com carregamento distribuído calcu-
lada usando FEM. . . . . . . . . . . . . . . . . . . . . . . . . 106
4.21 Carregamento equivalente devido a uma carga distribuída de
forma trapezoidal em um elemento. . . . . . . . . . . . . . . . 107
4.22 Elemento de frame no plano. . . . . . . . . . . . . . . . . . . . 109
4.23 Frame no plano com carregamento concentrado. . . . . . . . . 113
4.24 Malha do frame. . . . . . . . . . . . . . . . . . . . . . . . . . . 114
4.25 Frame no plano com carregamento distribuído. . . . . . . . . . 123
4.26 Treliça 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
4.27 Treliça 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
4.28 Treliça 3. Considere d = 2.0 m e L = 1000 N. . . . . . . . . . 132
4.29 Treliça 4. Considere a = 3.0 m . . . . . . . . . . . . . . . . . . 132
4.30 Treliça 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133
4.31 Treliça 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
4.32 Treliça 7. α = 30o e P = 500 N. . . . . . . . . . . . . . . . . . 134
4.33 Treliça 8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
4.34 Viga 1 com w0 = 10 kN/m, L = 5 m, a = 2 m, b = 4 m e
M = 20 kNm. . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
4.35 Viga 2 com L = 4 m e w0 = 10 kN/m. . . . . . . . . . . . . . 135
4.36 Viga 3 com a = 1 m, L = 4 m e w0 = 10 kN/m. . . . . . . . . 136
4.37 Viga 4 com a = 2 m, b = 3 m, L = 4 m, P = 10 kN e w0 = 10
kN/m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
4.38 Viga 5 com a = 1 m e w0 = 10 kN/m. . . . . . . . . . . . . . . 137
4.39 Viga 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
4.40 Viga 7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
7
4.41 Viga 8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
4.42 Frame 1, q = 10 kN/m, L = 2 m, E = 210 GPa, A = 4× 10−2
m2 e I = 4× 10−4 m4. . . . . . . . . . . . . . . . . . . . . . . 139
4.43 Frame 2, q = 10 kN/m, L = 2 m, E = 210 GPa, A = 4× 10−2
m2 e I = 4× 10−4 m4. . . . . . . . . . . . . . . . . . . . . . . 139
4.44 Frame 3, P = 12 kN, q = 3 kN/m, E = 200 GPa, A =
4.95× 10−3 m2 e I = 125.3× 10−6 m4. . . . . . . . . . . . . . 139
4.45 Frame 4 com barras quadradas de 200 mm × 200 mm, E = 10
GPa e P = 20 kN. . . . . . . . . . . . . . . . . . . . . . . . . 140
4.46 Frame 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
4.47 Frame 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
4.48 Frame 7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141
4.49 Barragem de concreto. . . . . . . . . . . . . . . . . . . . . . . 141
5.1 Esquema de aceleração média constante de Newmark. . . . . . 153
5.2 Resposta para o sistema. . . . . . . . . . . . . . . . . . . . . . 160
5.3 Resposta para o sistema. . . . . . . . . . . . . . . . . . . . . . 161
5.4 1.o e 2.o Modo de vibrar da viga calculado analiticamente. . . 162
5.5 3.o e 4.o Modo de vibrar da viga calculado analiticamente. . . 162
6.1 Domínio de solução de um problema 2D. . . . . . . . . . . . . 165
6.2 Elemento finito retangular com quatro nós. . . . . . . . . . . . 173
6.3 Elemento triangular com três nós. . . . . . . . . . . . . . . . . 179
7.1 Simulação em realidade virtual de ingestão de água nos moto-
res de um avião. . . . . . . . . . . . . . . . . . . . . . . . . . . 185
7.2 Teste em vôo em condição real. . . . . . . . . . . . . . . . . . 186
8
Lista de Tabelas
3.1 Graus de liberdade (DOF) necessário nos nós. . . . . . . . . . 35
3.2 Viscosidade µ do fluído. . . . . . . . . . . . . . . . . . . . . . 53
4.1 Deslocamentos nodais (mm). . . . . . . . . . . . . . . . . . . . 70
4.2 Forças externas aplicadas nos nós (kgf). . . . . . . . . . . . . . 70
4.3 Deformações �, tensões σ e forças internas axiais F nas barras. 71
77
4.5 Deslocamentos transversais e rotações na viga. . . . . . . . . . 96
4.6 Esforços atuantes nos nós do frame. . . . . . . . . . . . . . . . 122
4.7 Esforços atuantes nos nós do frame com carregamento distri-
buído. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129
5.1 Algumas frequencias naturais e fatores de amortecimento na
viga. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
9
Capítulo 1
Introdução
O método dos elementos finitos (FEM)1 é uma ferramenta numérica po-
derosa para resolver equações diferencias parciais. Muitos problemas físicos e
de engenharia em meios contínuos são descritos por equações diferenciais par-
ciais. A solução destes problemas na sua forma analítica (fechada) de forma
exata só é possível para sistemasmuito simples. Assim, para sistemas mais
complexos envolvendo geometrias e condições de contorno mais sofisticadas
não é possível se obter uma solução exata. Nestes casos deve-se optar por
procedimentos de aproximação com precisão aceitável para a aplicação de
engenharia em questão.
Inúmeros métodos de precisão para solução destes problemas são usa-
dos em engenharia entre eles pode-se destacar: método dos elementos de
contorno, método das diferenças finitas, método dos volumes finitos, mé-
todo de Galerkin, método de Rayleigh-Ritz e o método dos elementos finitos.
Deve ficar claro ao estudante que nenhum destes métodos pode ser conside-
rado superior ao outro. Isto depende do tipo de aplicação, solução desejada,
capacidade computacional, etc. que um engenheiro tem em mãos no mo-
mento de resolver um problema de engenharia. O FEM acabou se tornando
o mais popular de todos, sobretudo pelo aparecimento de diferentes pacotes
de software comercias sobre o assunto, como por exemplo o ANSYS, NAS-
TRAN/PATRAN, ADAMS, ABAQUS, etc.
A ideia básica do FEM é realizar uma divisão do domínio de integração de
uma estrutura ou sistema de interesse em um conjunto de pequenas regiões,
chamadas de elementos finitos transformando o domínio de contínuo para
discreto. Esta divisão do domínio é conhecida como malha ou grid, que nada
mais é do que o conjunto de elementos finitos resultante da discretização. A
malha é formada de elementos compostos de faces e nós, que são pontos de
1Finite Element Methods que vem do inglês.
10
intersecção e ligação entre os elementos. A grande "sacada"do FEM é não
buscar uma função admissível que satisfaça as condições de contorno para
todo o domínio, o que pode ser praticamente impossível em um problema
complexo, e sim buscar estas soluções em cada elemento separadamente.
Suponha que o funcional para um elemento seja Ψi, sua soma sobre a malha
com n elementos corresponde ao funcional de todo o domínio:
Ψ =
n∑
i=1
Ψi (1.1)
Para cada um dos elementos i existe uma função de interpolação (aproxi-
madora) u de ordemm descrita em função dos nós dos elementos (parâmetros
nodais αj) e por funções de forma (φ). A função de interpolação é descrita
como:
u =
m∑
j=1
αjφj (1.2)
O funcional da eq. (1.1) fica sendo descrito por:
Ψ(αj) =
n∑
i=1
Ψ(αj)i (1.3)
Aplicando as condições de estacionariedade geral leva um sistema de equa-
ções algébricas lineares. A solução do sistema de equações fornece os valores
dos parâmetros nodais αj. Os parâmetros nodais podem estar associados a
deslocamentos, forças internas, tensões, temperaturas, pressão, etc. e dep-
dende da formulação do elemento usado. Todos estes pontos serão melhor
discutidos no decorrer deste curso. Porém de antemão o que se espera é que o
aluno compreenda que o FEM é uma busca por uma solução local que possa
ser generalizada para todo o domínio. Várias abordagens do método FEM
são usadas.
1.1 Exemplos de aplicação
O FEM têm inúmeras aplicações nos diferentes ramos da ciência, em
especial em aplicações estruturais. Historicamente, as primeiras utilizações
de FEM em engenharia foram em aplicações aeronáuticas e de estruturas
civis, daí o grande avanço tecnológico de FEM nas empresas deste setor.
Seria impossível o Brasil atingir um alto nível de competência em projetos
de aeronaves sem o uso consistente de ferramentas envolvendo elementos
finitos.
11
Entre as áreas que usam FEM em projeto e análise se destacam:
• Estruturas oceânicas e navios.
• Veículos rodoviários e ferroviários.
• Hidrogeradores.
• Estruturas aeroespaciais e aviões.
• Mecânica estrutural.
• Mecânica dos fluídos computacional.
• Condução de calor.
• Eletromagnetismo.
A lista acima é imensa e serve apenas para mostrar as aplicações básicas.
Uma vez que FEM envolve ferramentas matemáticas das mais simples (en-
volvendo algebra vetorial) até as mais avançadas (como teoremas integrais)
o uso de pacotes comercias, como o NASTRAN, para análise é muito corri-
queiro. Em virtude do conhecimento que estes programas contém por trás de
seu código fonte, o seu preço é alto, o que faz com que apenas empresas de
grande porte tenham condições de ter as licenças comerciais destes softwa-
res. Contudo, deve ficar claro que um engenheiro que não sabe modelar um
problema via FEM sem o computador não saberá como proceder tendo uma
máquina e os mais avançados dos programas. As facilidades gráficas de fer-
ramentas CAD, CAE, CAM traz a sensação que para ser um engenheiro de
projetos basta "decorar"meia dúzia de comandos para se dizer especialista
em FEM. Porém, isto é um conceito errado. O autor do livro [4] cita um
exemplo interessante: Imagine que você está muito doente e procura um mé-
dico que não é um grande especialista na sua enfermidade. O médico diz
para não se preocupar, pois ele tem um programa onde "basta"digitar na
entrada os sintomas que ele fornece na saída os diagnósticos com a profilaxia
adequada. Se você não tem algum problema psiquiátrico grave, provavel-
mente você não irá confiar neste médico. Agora já imaginou entrar em uma
aeronave projetada por um engenheiro com está visão! Sendo assim, o ideal
é o estudante ter uma base sólida em FEM conhecendo os princípios básicos
do método. Isto permite que ele use pacotes comerciais com maior rigor de
análise e que saiba interpretar as soluções e gráficos e, por que não, ser ca-
paz de programar seus elementos em rotinas próprias. Quem usa softwares
e nunca estudou FEM de forma convencional não se pode dizer que saiba o
que é o método.
12
Apenas para ilustrar uma aplicação prática, a fig. (1.1) mostra a análise
da região de um bloco de motor [4].O bloco é modelado com elementos sóli-
dos usando elementos tetraédricos parabólicos em virtude de sua geometria
complexa. A meta foi calcular o panorama de tensões, que é mostrado na
mesma fig. na estrutura visando analisar sua resistência mecânica e tole-
rância a falha. Assim, o FEM é uma ferramenta útil e imprescindível em
projetos modernos de engenharia.
1.2 Etapas na solução de um problema via
FEM
O FEM é um procedimento bem metódico dividido em várias etapas:
1. Desenvolvimento das equações do elemento.
2. Discretização do domínio de solução dentro de uma malha de elementos
finitos.
3. Montagem das equações do elemento.
4. Introdução das condições de contorno (restrições físicas e geométricas).
5. Solução para os nós desconhecidos.
6. Cálculo da solução e das quantidades (grandezas) em cada elemento.
Muitas vezes estas etapas são misturadas ou trabalhadas de forma si-
multânea. Durante este curso cada uma destas etapas serão estudadas com
calma. O foco de problemas abordados será direcionado a aplicações de en-
genharia mecânica envolvendo problemas estruturais e, eventualmente, tér-
micos.
1.3 Discretização por elementos finitos
O primeiro passo de um método FEM é escolher qual elemento utilizar.
Estes elementos podem ser do tipo unidimensional (1D), como os elementos
de barra e viga, bidimensional (2D), como os elementos de placa, e tridimen-
sionais (3D) como os elementos sólidos. A fig. (1.2) mostra alguns exemplos
destes elementos. Os programas comerciais de FEM possuem bibliotecas com
centenas de elementos finitos que podem ser empregados em simulações. Ape-
sar do senso comum acreditar que elementos 3D são sempre superiores aos
13
Fig. 1.1: Modelo FEM de um bloco de motor [4].
elementos unidimensionais, isto não é verdade. A escolha de um elemento
deve ser condicionada ao tipo de geometria e de aproximação de solução
14
que se deseja obter. Formulações de alguns elementos podem ter resultado
superiores a de outros elementos no processo de aproximação. Neste curso
pretende-se estudar alguns elementos básicos e clássicos como os elementos
de barra, viga de Euler-Bernoulli e de placa de Kirchoff.
Fig. 1.2: Tipos e formas de elementos [2].
Como járessaltado anteriormente, no FEM uma solução aproximada é
assumida em cada nó através de uma função de interpolação, que envolve
funções de forma e parâmetros nodais. Cada nó tem seus graus de liberdade
que podem ser deslocamento, temperatura, pressão, voltagem, etc. que nor-
malmente são incógnitas. O processo resultante da montagem dos elementos
finitos no domínio global conduz em um sistema de equações algébricas de
grande dimensão. Do ponto de vista matemático, FEM é uma forma espe-
cial dos métodos de aproximação de Galerkin e Rayleigh-Ritz utilizados para
encontrar soluções de equações diferenciais.
A qualidade (acurácia) da aproximação é diretamente proporcional a
quantidade de elementos usados. O custo computacional também é ligado
ao número de elementos, uma vez que o sistema de equações se torna maior.
Em um problema FEM uma estrutura pode ter uma malha com múltiplos
tamanhos de discretização em regiões onde se necessita de maior acurácia
(malha mais refinada). Já regiões onde não se tem muito interesse podem
usar malhas mais grosseiras. Os contornos curvilíneos são exemplos onde
malhas finas devem ser usadas, conforme a fig. 1.3.
A escolha de um elemento, numero de elementos, etc. deve ser pautada no
tipo de solução e capacidade computacional disponível. Uma boa alternativa
é começar com malhas grosseiras para se ter noção do tipo de solução obtido
e depois refinar a malhar conforme desejo, obtendo assim uma economia
computacional e de tempo que pode resultar em maior produtividade.
15
Fig. 1.3: Diferença entre o contorno físico e a geometria do contorno do
modelo de elementos finitos [2].
Um erro comum que usuário inexperientes em FEM cometem se refere a
conectividade dos nós. Uma malha deve ser formada por elementos que se
conectam através de nós. A interface dos nós portanto deve ser tal que per-
mita que a malha seja devidamente "fechada"entre elementos adjacentes. A
fig. (1.4) apresenta um exemplo de malha válida e inválida para um elemento
com quatro nós. A malha fica inválida quando usando uma configuração com
três elementos, uma vez que o nó 4 dos elementos (2) e (3) não está conectado
com nenhum nó do elemento (1).
Fig. 1.4: Malhas válida e inválida para elementos com 4 nós [2].
Outra forma de economizar tempo seria analisando as simetrias entre con-
dições de contorno e domínio de solução. Se a malha for construída de forma
simétrica o resultado também o seria. Um exemplo de estrutura simétrica é
formado por esquadria metálicas de galpões formados por conjuntos idênticos
de treliças. Porém deve-se ter cuidado, pois se a estrutura for simétrica a
alguma linha de centro, a escolha de um elemento não adequado pode fazer
com que este procedimento não seja correto.
16
Capítulo 2
Fundamentos Matemáticos
Básicos de FEM
Para se aplicar FEM é necessária uma base sólida em procedimentos
matemáticos que vão dos mais simples, como manipulação de matrizes, até os
mais avançados, envolvendo por exemplo teoremas de cálculo vetorial. Este
capítulo tem como meta revisar alguns destes conceitos que serão utilizados
no decorrer do curso. Não será dado nenhum rigor matemático, muito menos
provas dos teoremas mostrados. O aluno interessado pode consultar livros
básicos de cálculo para reforçar alguns conceitos.
2.1 Análise vetorial
Para descrever uma grandeza vetorial, como uma força, deslocamento,
fluxo, etc. é necessário se definir três componentes: módulo, direção e sen-
tido. Um vetor a pode ser descrito em coordenadas cartesianas em função
de vetores unitários (i, j, k):
a = axi + ayj + azk (2.1)
Os conceitos de cálculo vetorial são muito usados em FEM. A ideia bá-
sica do cálculo vetorial é considerar cada ponto no espaço como uma função
vetorial, o que forma um campo vetorial. Um campo vetorial pode ser um
deslocamento, fluxo de um fluído, força gravitacional ou eletromagnética, etc.
Já um campo escalar significa associar cada ponto no espaço com um funci-
onal escalar. Um exemplo de campo escalar é um campo de temperatura em
um ponto no espaço, campo de pressão, etc.
O operador diferencial (del)∇ é muito usado para definir operações mate-
máticas fundamentais em campos escalares e vetoriais. O operador diferencial
17
∇ é dado por:
∇ ≡ ∂
∂x
i +
∂
∂y
j +
∂
∂z
k (2.2)
e representa um operador diferencial de 1.o ordem. Um operador de 2.o
ordem é conhecido como Laplaciano e dado por:
∇2 ≡ ∂
2
∂2x
i +
∂2
∂2y
j +
∂2
∂2z
k (2.3)
O operador ∇ é usado para definir três operações básicas envolvendo
campos escalares e vetoriais:
• Gradiente
• Divergente
• Rotacional
Estas operações são usadas na definição de teoremas fundamentais de
integrais de vetores tais como o Teorema da Divergência e o Teorema de
Green-Gauss. Estes dois teoremas são a base matemática para compreender
o método de Galerkin, que por sua vez é uma das bases fundamentais de
FEM. Antes de aprofundar nestas questões é interessante revisar operações
básicas envolvendo vetores, que são mostradas a seguir.
2.1.1 Produto escalar
O produto escalar entre dois vetores a e b é definido por:
a · b = |a| (|b|cos(θ)) = b · a = |b| (|a|cos(θ)) (2.4)
sendo | · | o módulo do vetor e θ o ângulo entre eles. O resultado da
operação de produto escalar é um escalar. Note que i · i = 1 e que i · j = 0 e
assim por diante, uma vez que os vetores unitários definem uma base e são
ortonormais (θ = 90o).
2.1.2 Produto vetorial
O produto vetorial entre dois vetores a e b é definido por:
a× b =
∣∣∣∣∣∣
i j k
ax ay az
bx by bz
∣∣∣∣∣∣ (2.5)
18
O resultado da operação de produtor vetorial é um vetor perpendicular
ao plano onde estão contidos os vetores a e b. Note que i × i = 0 e que
i× j = k. Importante observar que i× j 6= j× i.
2.1.3 Gradiente
O gradiente de uma função escalar Φ(x, y, z)1 é dado por:
∇Φ =
(
∂
∂x
i +
∂
∂y
j +
∂
∂z
k
)
Φ =
(
∂Φ
∂x
i +
∂Φ
∂y
j +
∂Φ
∂z
k
)
(2.6)
Note que o resultado do gradiente é um vetor. Esta operação representa
uma diferença entre níveis de um campo escalar, representando a variação de
uma grandeza escalar por unidade de espaço. O significado físico pode ser
interpretado como a diferença de temperatura nas faces de um bloco, para
este tipo de aplicação.
2.1.4 Divergente
O divergente já é uma operação envolvendo um campo vetorial dado por
uma função vetorial do tipo a(x, y, z)2e calculado por:
∇ · a =
(
∂
∂x
i +
∂
∂y
j +
∂
∂z
k
)
· (axi + ayj + azk) (2.7)
o que leva a seguinte expressão:
∇ · a = ∂ax
∂x
+
∂ay
∂y
+
∂az
∂z
(2.8)
note que ∇ · a 6= a · ∇ uma vez que o operador ∇ deve agir sobre a.
O divergente pode ser interpretado como um escalar que mostra se um
campo vetorial está se expandindo ("fonte") ou comprimindo ("ralo"). É
uma medida de magnitude da dispersão de um campo vetorial.
2.1.5 Rotacional
O rotacional representa um vetor resultante entre o produto vetorial en-
volvendo o operador diferencial ∇ e um campo vetorial a(x, y, z). Seu re-
1Esta função pode representar fisicamente um campo de pressão, temperatura, etc. no
espaço.
2Que pode representar fisicamente um campo de deslocamento, fluxo de um fluído,
força, etc.
19
sultado pode ser escrito na forma de um tensor cartesiano. Esta operação é
calculada como:
∇× a =
∣∣∣∣∣∣
i j k
∂
∂x
∂
∂y
∂
∂z
ax ay az
∣∣∣∣∣∣ (2.9)
O rotacional tem este nome pois esta operação representa uma trans-
formação linear de coordenadas (rotação) do campo vetorial a(x, y, z) que
visa observar suas características nestas novas coordenadas. A representação
acima só vale para representações em coordenadas retangulares.
2.1.6 Teorema da divergência
O Teorema da divergência é definido como:∫
V
(∇ · a) dV =
∫
S
(a · n) dS (2.10)
sendo V um volume, uma superfície de área S e n um vetorortonormal
à esta superfície S3. O teorema da divergência relaciona o divergente total
de um campo vetorial a em um volume V com o fluxo total deste campo
vetorial atravessando uma superfície S.
2.1.7 Teorema de Green-Gauss
Muitos problemas de engenharia podem ser escritos em uma forma unidi-
mensional4 e considerando as derivadas de funções escalares Ψ e Φ com um
valor k constante. Assim:
d
dx
(
k
dΦ
dx
Ψ
)
= k
d2Φ
dx2
Ψ + k
dΦ
dx
dΨ
dx
(2.11)
Aplicando integral de ambos os lados de a até b temos:
k
∫ b
a
d
dx
(
dΦ
dx
Ψ
)
dx = k
∫ b
a
d2Φ
dx2
Ψdx+ k
∫ b
a
dΦ
dx
dΨ
dx
dx (2.12)
Notanto que o lado esquerdo da eq. (2.12) forma uma integral perfeita
tem-se que:
3a · n pode representar um fluxo.
4O próximo capítulo irá revisar alguns destes problemas.
20
k
∫ b
a
d
(
dΦ
dx
Ψ
)
= k
∣∣∣∣dΦdxΨ
∣∣∣∣b
a
(2.13)
Substituindo a eq. (2.13) em (2.12) e rearranjando tem-se:
k
∫ b
a
d2Φ
dx2
Ψdx = k
∣∣∣∣dΦdxΨ
∣∣∣∣b
a
− k
∫ b
a
dΦ
dx
dΨ
dx
dx (2.14)
Considerando que a = βb o teorema da divergência, eq. (2.10), pode ser
reescrito como: ∫
V
(∇ · βb) dV =
∫
S
(βb · n) dS (2.15)
Uma vez que ∇ · βb = β∇ · b +∇β · b tem-se que:∫
V
(β∇ · b) dV =
∫
S
(βb · n) dS −
∫
V
(∇β · b) dV (2.16)
A eq. (2.16) é um resultado clássico do teorema de Green-Gauss. Em
FEM a eq. (2.14) é uma extensão da eq. (2.16) sendo que Φ e Ψ são
matrizes representando funções de interpolação (funções aproximadores) dos
elementos empregados em uma discretização. Em um momento oportuno
nos próximos capítulos este ponto será revisto com mais detalhes.
2.2 Análise matricial
O FEM emprega bastante em sua formulação o uso e a manipulação
de matrizes. Muitos destas manipulações são triviais, como por exemplo,
cálculo de determinante, mínimo de uma matriz, cofatores, adjuntos, etc.
Outros são mais avançados, como por exemplo, técnicas para inversão de
matrizes visando solucionar sistemas lineares de grande dimensão.
Uma operação usada em FEM se refere a eliminação de linhas e colunas
de uma matriz, que corresponde na prática a aplicação de uma condição de
contorno ou restrição no sistema em estudo. Suponha uma matriz A dada
por:
A =
 a11 a12 a13a21 a22 a23
a31 a32 a33
 (2.17)
Se uma restrição for imposta de tal forma que a segunda linha e coluna
sejam eliminadas temos uma matriz M22 dada por:
21
M22 =
[
a11 a13
a31 a32
]
(2.18)
Este conceito também é usado para cálculo do cofator Cij:
Cij = (−1)i+j |Mij| (2.19)
Já o adjunto de uma matriz é Aij = CTij .
Uma aplicação comum em FEM é ter que resolver sistemas lineares do
tipo:
Ax = f (2.20)
Onde o vetor x representa as incógnitas do problemas que são os graus de
liberdade em cada nó de um elemento (por exemplo, deslocamento), a matriz
A os parâmetros conhecidos representando uma matriz de rigidez e o vetor f
representando as fontes ou forças atuantes. A solução deste problema é feita
a partir da inversão da matriz de rigidez:
A−1Ax = Ix = x = A−1f (2.21)
Porém este método é ineficiente para solucionar sistemas de grandes equa-
ções. Uma maneira mais efetiva e elegante é propor uma decomposição da
matriz de rigidez A, como por exemplo, o método de eliminação de Gauss.
Exemplo 2.1 Use o método de eliminação de Gauss para resolver o sistema
simultâneo de equações:
4x1 + 2x2 − 2x3 − 8x4 = 4
x1 + 2x2 + x3 = 2
0.5x1 − x2 + 4x3 + 4x4 = 10
−4x1 − 2x2 − x4 = 0
Este sistema de equações pode ser descrito na forma matricial como:
4 2 −2 −8
1 2 1 0
0.5 −1 4 4
−4 −2 0 1


x1
x2
x3
x4
 =

4
2
10
0
 (2.22)
Primeiro é dividido a 1.o linha por 4 e subtraindo esta nova linha pela 2.o
linha. Na sequência a nova linha 1 é dividida por 0.5 e subtraída da linha 3.
Por fim, a linha 1 é dividida por -4 e subtraída da linha 4. O resultado é:
22

1 0.5 −0.5 −2
0 1.5 1.5 2
0 −1.25 4.25 5
0 0 −2 −7


x1
x2
x3
x4
 =

1
1
9.5
4
 (2.23)
Agora neste novo sistema a linha 2 é dividida por 1.5, a nova linha 2 é
multiplicada por -1.25 e subtraída da linha 3. Como um zero já apareceu na
linha 4 nenhuma modificação é exigida. Este resultado é:
1 0.5 −0.5 −2
0 1 1 1.3333
0 0 5.5 6.6667
0 0 −2 −7


x1
x2
x3
x4
 =

1
0.6667
10.3333
4
 (2.24)
Por fim, a linha 3 é dividida por 5.5. Multiplicando esta nova linha por
3 por -2 é subtraindo da linha 4:
1 0.5 −0.5 −2
0 1 1 1.3333
0 0 1 1.2122
0 0 0 −4.578


x1
x2
x3
x4
 =

1
0.6667
1.8788
7.7576
 (2.25)
Agora a solução do sistema é trivial e é dada por: x1 = 0.0794, x2 =
−1.0066, x3 = 3.9338 e x4 = −1.6954.
2.3 Equações diferenciais
Como já discutido no capítulo 1 deste texto o FEM é uma formulação
para solucionar de forma numérica e com aproximações uma equação diferen-
cial. Sendo assim, é primordial que o engenheiro saiba modelar fisicamente
o seu problema com o conhecimento necessário para construir este sistema
de equações diferenciais.
Apesar desta vertente de utilização de FEM, as aplicações clássicas nor-
malmente envolvem estruturas civis e aeroespaciais baseadas em elementos
simples, como barra, viga e placa. Estes problemas podem ser descritos por
equações diferenciais parciais (problema de valor de contorno). Nestes exem-
plos, o FEM pode ser formulado a partir de métodos de energia envolvendo
funções de variação (método de Lagrange) sem necessariamente considerar
equações diferenciais. Esta representa uma das formas de abordagem de
FEM. Porém, neste texto inicialmente será dada uma abordagem da for-
mulação de FEM diretamente nas equações diferenciais juntamente com a
metodologia se empregando métodos de energia.
23
Neste capítulo vale apenas lembrar que a maioria dos problemas de enge-
nharia podem ser escritos através da equação (para o caso unidimensional5):
d
dx
α(x)A(x)
dΦ(x)
dx
+ C(x)A(x) = 0 (2.26)
Sendo α(x) um parâmetro do material, C(x) uma fonte externa e A(x) a
área da secção transversal. Se estes parâmetros forem variantes significa que
o sistema varia de elemento a elemento. A forma básica é assumir homoge-
neidade, assim a eq. (2.26) torna-se:
α
d2Φ
dx2
+ C = 0 (2.27)
Inúmeros métodos analíticos podem ser usados para solucionar este tipo
de problema, como separação de variáveis, coeficientes desconhecidos, trans-
formada de Laplace, etc. O capítulo 3 começa com uma revisão rápida dos
principais problemas de engenharia com suas respectivas equações.
2.4 Tensores cartesianos
Quando se trabalha com FEM envolvendo sistemas complexos, normal-
mente são encontradas equações de grande dimensão. Nestes casos a notação
de subscritos6 pode ser útil. Em primeiro lugar é preciso lembrar a definição
de tensor. Tensor é uma grandeza que precisa de nove elementos para poder
ser completamente conhecida. Em alguns casos com 6 elementos é possível
descrever um tensor, como por exemplo, no caso de um estado de tensões,
onde as tensões cisalhantes no mesmo plano são iguais. A notação tenso-
rial pode ser usada como forma de propor uma notação compacta para uma
notação vetorial. Um vetor descrito nesta notação é um tensor de primeira
ordem.
Imagine um vetor f escrito em função do sistema de coordenadas (x, y, z):
f = fxi + fyj + fzk (2.28)
Agora em vez do sistema de coordenadas (x, y, z) imagine um equivalente
(x1, x2, x3). Neste novo sistema de coordenadas este vetor é descrito como:
f = f1i + f2j + f3k (2.29)
5Estas equações também podem ser escritas de forma parcial quando envolvem mais
de duas variáveis
6Também conhecida como notação de Einstein ou notação tensorial.
24
Em uma notação tensorialeste vetor pode ser dado por:
f = fi, i = 1, 2, 3 (2.30)
2.5 Exercícios
Ex. 2.1 Considere um campo vetorial dado por a = (−y, xy, z). Calcule o
divergente e o rotacional deste campo vetorial.
Ex. 2.2 Sendo Φ e Ψ dois campos escalares mostre que ∇ (ΦΨ) = Φ∇Ψ +
Ψ∇Φ
Ex. 2.3 Resolva o seguinte sistema linear usando o método de eliminação
de Gauss. 
2 1 2 −3
2 −2 1 −4
1 0 2 −3
4 4 −4 1


x1
x2
x3
x4
 =

0
5
−4
−6
 (2.31)
Ex. 2.4 Reescreva a solução do exercício 2 usando notação de tensores car-
tesianos.
Ex. 2.5 Dada uma função escalar J(u) = k(Bu)2, onde k, B e u são
definidos como: k é uma matriz n × n, B é uma matriz n ×m e u é uma
matriz m× p. Escreva J(u) usando uma equação matricial.
Ex. 2.6 Calcule ∂J(u)/∂u para a equação matricial do problema anterior.
25
Capítulo 3
Elementos Finitos
Unidimensionais
FEM representa uma solução aproximada de um problema de valor de
contorno (PVC) descrito por uma equação diferencial. Equações que go-
vernam fenômenos de interesse em engenharia são geralmente obtidas por
equações de balanço e equações constitutivas. Neste capítulo todas as bases
matemáticas para FEM visando resolver este tipo de problema são apresen-
tadas usando casos unidimensionais. O caso a ser estudado em detalhes neste
capítulo é o problema de elasticidade em barras sujeitas a cargas axiais. O
capítulo começa apresentando a solução exata deste problema. Na sequência
o método de aproximação de Galerkin é mostrado em detalhes preparando
para a sua aplicação no método de elementos finitos visando o seu uso em pro-
blemas unidimensionais. Por fim, outros problemas clássicos de engenharia
são apresentados, como o problema de transferência de calor e massa, fluxo
de potência e eletricidade. Extensões para estes problemas podem ser pro-
postas a partir dos conceitos exemplificados ao longo deste capítulo. Sendo
assim, no final do capítulo exercícios são propostos ao estudante envolvendo
diferentes aplicações de FEM em problemas unidimensionais de interesse em
engenharia.
3.1 Solução exata de problemas de deformação
axial em barra uniforme
É possível integrar uma equação diferencial para obter a sua solução exata
uma vez conhecidas as condições de contorno. Nesta seção é formulado o
problema básico que será usado como exemplo benchmark no decorrer deste
capítulo.
26
Um problema unidimensional em elasticidade, que provavelmente já deve
ter sido estudado por todos em um curso básico de resistência dos materi-
ais, é descrito pelo balanço de forças em uma barra elástica sujeita a uma
deformação linear em termos de área A, tensão normal σ e força axial f . O
balanço de forças é dado por:
d[σ(x)A(x)]
dx
+ f(x)A(x) = 0 (3.1)
Já a equação constitutivas representa relações do material, neste caso a
Lei de Hooke envolvendo a deformação �(x) e o módulo de Young E(x), e é
dada por:
σ(x) = E(x)�(x) (3.2)
sendo a deformação �(x) relacionada ao deslocamento u(x) da barra:
�(x) =
du(x)
dx
(3.3)
Assim com a Lei de Hooke, a eq. (3.2) é escrita como:
σ(x) = E(x)
du(x)
dx
(3.4)
Considere agora o exemplo descrito pela barra engastada-livre, fig. (3.1).
Esta barra uniforme é sujeita a uma carga estática na extremidade livre. A
sua carga é linearmente variante q(x) = cx, sendo c uma constante. Esta es-
trutura é descrita pelo seguinte problema de valor de contorno (PVC) usando
as eqs. (3.1) e (3.4)1:
EA
d2u
dx2
+ cx = 0; 0 < x < L (3.5)
u(0) = 0; EA
du(L)
dx
= P
As condições de contorno em u(x) são essenciais ou geométricas e as
condições de contorno em σ são naturais.
A solução exata para este problema pode ser facilmente obtida integrando
a eq. (3.5) duas vezes e usando as condições de contorno para avaliar as
constantes de integração. Integrando ambos os lados da equação uma vez,
obtém-se:
1Assumindo que o módulo de Young e a área são constantes em todo o domínio.
27
Fig. 3.1: Barra uniforme carregada axialmente [2].
EA
du
dx
+
cx2
2
= C1 (3.6)
sendo C1 uma constante de integração. Integrando mais uma vez:
EAu(x) +
cx3
6
= C1x+ C2 (3.7)
sendo que C2 é outra constante de integração. Rearranjando estes termos:
u(x) =
1
EA
(
C1x+ C2 − cx
3
6
)
(3.8)
Usando as condições de contorno e a eq. (3.8):
u(0) = 0⇒ C2 = 0 (3.9)
e
EA
du(L)
dx
= P ⇒ 1
6
(−3cL2 + 6C1) = P (3.10)
que conduz a:
C1 =
1
2
(2P + cL2) (3.11)
Assim a solução exata para este problema é:
u(x) =
x(6P + 3cL2 − cx2)
6EA
(3.12)
Uma vez tendo esta solução pode-se usar a eq. (3.4) para calcular a
distribuição de tensão ao longo da barra. Para este problema a solução
28
exata é calculada facilmente. Infelizmente, nem todos os problemas reais de
engenharia podem ser descritos de forma tão simples assim. Com o propósito
de propor uma solução aproximada a próxima seção apresenta o método de
Galerkin.
3.2 Aproximação via método de Galerkin
No método de Galerkin assume-se uma forma geral de solução para um
PVC. Esta solução aproximada, denotada por u˜(x), deve ser tal que o erro
entre a solução exata e a aproximada seja o menor possível. Esta solução
pode ter qualquer forma. No geral, são mais usadas soluções aproximadas
em formas polinomiais do tipo:
u˜(x) = a0 + a1x+ a2x
2 + · · ·+ anxn (3.13)
sendo a0, a1,. . .,an parâmetros desconhecidos. É importante observar
que uma vez que u˜(x) é uma aproximação, esta solução não irá satisfazer a
equação diferencial para todos os valores de x. Substituindo está solução no
PVC da deformação axial da barra obtém-se um erro e(x):
e(x) =
d
dx
(
AE
du˜
dx
)
+ q(x) 6= 0 (3.14)
O erro total, chamado de resíduo, para o domínio inteiro da solução pode
ser obtido integrando o erro e(x) sobre todo o domínio. Entretanto, uma
vez que erros positivos e negativos pode ser cancelados nesta integração é
interessante o uso de funções de ponderação wi(x) para os i = 0, 1, . . . , n,
sendo n o número de parâmetros desconhecidos. Assim o resíduo pode ser
dado por: ∫ xl
x0
e(x)wi(x)dx = 0; i = 0, 1, . . . , n (3.15)
O desejo é que este resíduo seja nulo. Está forma é conhecida como forma
fraca. A equação diferencial é a forma forte pois exige que o erro e(x) seja
nulo em todo o domínio 0 < x < L.
O método mais popular nas aplicações de elementos finitos é o método de
Galerkin. Neste método a funções de ponderações wi(x) da função erro e(x)
são definidas como derivadas parciais das soluções assumidas em relação aos
parâmentros desconhecidos ai2:
2No caso descrito a seguir irá ser considerado a derivada parcial, uma vez que em um
caso mais geral u pode ser função de várias variáveis.
29
wi(x) =
∂u
∂ai
; i = 0, 1, . . . , n (3.16)
Assim o método de Galerkin define o seguinte resíduo para as n equações
para solução dos parâmetros desconhecidos:∫ xl
x0
e(x)
∂u
∂ai
dx = 0; i = 0, 1, . . . , n (3.17)
Para a maioria das aplicações de engenharia o método de Galerkin for-
nece a mesma solução que outro método popular de aproximação, método
de Rayleigh-Ritz. Considerando nosso problema benchmark envolvendo a
deformação axial da barra, podemos calcular o resíduo para esta situação,
substituindo a eq. (3.14) na eq. (3.15):
∫ xl
x0
e(x)wi(x)dx =
∫ xl
x0
(
d
dx
(
AE
du˜
dx
)
+ q(x)
)
wi(x)dx = 0; (3.18)
i = 0, 1, . . . , n
Vamos assumir que A e E possam ser funções. Escrevendo as integrais
separadas:
∫ xl
x0
d
dx
(
AE
du˜
dx
)
wi(x)dx+
∫ xl
x0
q(x)wi(x)dx = 0 (3.19)
i = 0, 1, . . . , n
A primeira integral da eq. (3.19) contém derivação de termos de segunda
ordem em u˜, sendo assim deve-se aplicar integração por partes. Lembrando
que se tivermos duas funções f(x) e g(x), a integração por partes é definida
como:
∫ xl
x0
[
d
dx
(f(x))]
g(x)dx = f(xl)g(xl)− f(x0)g(x0)−
∫ xl
x0
[
d
dx
(g(x))
]
f(x)dx
(3.20)
Aplicando integração por partes na eq. (3.19) fornece:
A(xl)E(xl)
du˜(xl)
dx
wi(xl)− A(x0)E(x0)du˜(x0)
dx
wi(x0) (3.21)
−
∫ xl
x0
AE
du˜
dx
dwi
dx
dx+
∫ xl
x0
q(x)wi(x)dx
30
Os primeiros dois termos desta integral incorporam força ou condições de
contorno derivativas dentro da forma fraca. Se uma força Px0 é aplicada em
x0 tem-se que:
A(xl)E(xl)
du˜(xl)
dx
wi(xl) = Px0wi(x0) (3.22)
e se uma força Pxl é aplicada na extremidade l tem-se:
A(x0)E(x0)
du˜(x0)
dx
wi(x0) = Pxlwi(x0) (3.23)
Aplicando as eqs. (3.22) e (3.23) na eq. (3.21) a forma fraca pode ser
escrita por:
∫ xl
x0
E
du˜
dx
dwi
dx
Adx =
∫ xl
x0
q(x)wi(x)dx+ Pxlwi(xl) + Px0wi(x0) (3.24)
Para problemas estruturais a eq. (3.24) pode ser interpretada como o
bem conhecido princípio do trabalho virtual3. Se a função de ponderação
wi(x) for vista como um deslocamento virtual, então o lado esquerdo da eq.
(3.24) é o trabalho interno virtual total, desde que E deu
dx
= σx é a tensão
axial na barra e dwi
dx
é a deformação axial virtual na barra. Assim, a forma
fraca no método de Galerkin implica que quando é dado um deslocamento
virtual na barra, o trabalho externo virtual é igual ao trabalho interno virtual.
Este princípio é altamente usado em mecânica estrutural, seja em aplicações
estáticas ou dinâmica. Sua ampla aplicação é uma das principais razões que
fazem com que o método de Galerkin seja o mais popular no desenvolvimento
das equações de elementos finitos.
3.2.1 Aplicação do método de Galerkin na solução do
problema da barra axial
Após apresentado os passos básicos do método de Galerkin, está seção
apresenta sua aplicação para a solução do PVC descrito pela eq. (3.5). Várias
soluções aproximadas podem ser propostas. Uma vez que consideramos que
EA é constante, q(x) = cx, domínio 0 < x < L e as condições de contorno
são u(0) = 04 e EAu′(L) = P , a forma fraca específica deste problema é dada
por:
3Que dever ter sido estudado no curso básico de mecânica geral.
4E consequentemente w(0) = 0.
31
Pwi(L) +
∫ L
0
(
−AEdu˜
dx
dwi
dx
+ cxwi(x)
)
dx = 0 (3.25)
3.2.2 Solução linear
A solução mais simples para este problema é assumir que o campo de
deslocamento para esta barra pode ser aproximado por5:
u˜(x) = a0 + a1x (3.26)
Para satisfazer as condições de contorno:
u(0) = 0⇒ a0 + 0a1 = 0⇒ a0 = 0 (3.27)
Assim, a solução admissível é dada por:
u(x) = xa1 ⇒ du
dx
= a1 (3.28)
Substituindo a eq. (3.28) na eq. (3.25) obtem-se:
Pwi(L) +
∫ L
0
(
−AE(a1)dwi
dx
+ cxwi(x)
)
dx = 0 (3.29)
Existe somente um parâmetro desconhecido, a1. A função de ponderação
de Galerkin é dada por:
w1 =
∂u
∂a1
= x;
∂w1
∂x
= 1; w1(0) = 0; w1(L) = L (3.30)
Substituindo os resultados da eq. (3.30) na forma fraca, eq. (3.29),
obtém-se:
PL+
∫ L
0
(−AE(a1) + cx2) dx = 0 (3.31)
Resolvendo esta integração:
cL3
3
+ PL− EAa1L = 0 (3.32)
Resolvendo esta equação para a1 tem-se:
5Ressalta-se que a solução aproximada u˜(x) será denotada daqui para frente simples-
mente por u(x) para facilitar a notação.
32
a1 = −−cL
2 − 3P
3EA
(3.33)
Assim a solução linear aproximada é dada pela eq. (3.28) considerando o
resultado da eq. (3.33):
u(x) = a1x =
cL2 + 3P
3EA
x (3.34)
3.2.3 Solução quadrática
Uma solução melhor do que a solução linear pode ser obtida se utilizarmos
uma polinômio com mais termos, por exemplo, uma forma quadrática:
u(x) = a2x
2 + a1x+ a0 (3.35)
Para satisfazer as condições de contorno essenciais devemos ter que:
u(0) = 0⇒ a20 + a10 + a0 = 0⇒ a0 = 0 (3.36)
Assim, a solução admissível neste caso é
u(x) = a2x
2 + a1x⇒ du
dx
= 2a2x+ a1 (3.37)
Substituindo o resultado da eq. (3.37) na forma fraca, eq. (3.25) tem-se:
Pwi(L) +
∫ L
0
(
−AE (2a2x+ a1) dwi
dx
+ cxwi(x)
)
dx = 0 (3.38)
Agora existem dois parâmetros desconhecidos, a1 e a2, e necessitaremos
de duas equações para encontrá-los. Estas equações são obtidas a partir das
funções de ponderação de Galerkin:
w1 =
∂u
∂a1
= x;
∂w1
∂x
= 1; w1(0) = 0; w1(L) = L (3.39)
w2 =
∂u
∂a2
= x2;
∂w2
∂x
= 2x; w2(0) = 0; w2(L) = L
2 (3.40)
Substituindo w1 na forma fraca eq. (3.38)tem-se:
PL+
∫ L
0
(−AE (2a2x+ a1) + cx2) dx = 0 (3.41)
⇒ cL
3
3
− EAa2L2 + PL− EAa1L = 0 (3.42)
33
e substituindo w2 na eq. (3.38) tem-se:
PL2 +
∫ L
0
(−AE (2a2x+ a1) (2x) + cx3) dx = 0 (3.43)
⇒ cL
4
4
− 4
3
EAa2L
3 + PL2 − EAa1L2 = 0 (3.44)
Resolvendo estas duas equações obtém-se as expressões para os dois pa-
râmetros desconhecidos a1 e a2:
a1 = −−7cL
2 − 12P
12EA
(3.45)
a2 = − cL
4EA
(3.46)
Por fim, a solução quadrática aproximada é:
u(x) = a2x
2 + a1x =
(12P + cL(7L− 3x))
12EA
x (3.47)
A solução exata pode ser obtida utilizando uma aproximação cubica. Isto
é dado como exercício proposto ao estudante no final do capítulo.
3.3 Forma de elementos finitos de soluções as-
sumidas
A qualidade da aproximação via método de Galerkin depende da esco-
lha da solução assumida. Entretanto, não há um procedimento para guiar
esta escolha, especialmente para problemas de duas ou três dimensões. Além
disto, para ser admissível , estas soluções devem satisfazer as condições essen-
ciais de contorno. Para sistemas complicados, pode ser impossível a obtenção
de soluções admissíveis apropriadas. O FEM contorna estas dificuldades in-
troduzindo dois conceitos fundamentais:
(i) Domínio da solução é discretizado dentro de elementos finitos. O domí-
nio de solução é dividido dentro de vários subdomínios mais simples
chamados de elementos. Cada elemento tem uma geometria simples
tal que a solução aproximada seja facilmente escrita para um elemento.
Cada elemento cobre apenas um pequena parte do domínio de solução,
um polinômio de grau baixo pode, no geral, ser usado para descrever
a solução para cada elemento. A integral na forma funcional ou fraca
34
pode ser avaliada separadamente em cada elemento e somada (mon-
tada) junta como:
∫ l
0
(· · ·) =
∫ x2
x1=0
(· · ·) dx+
∫ x3
x2
(· · ·) dx+ · · · (3.48)
(ii) Coeficientes nas soluções assumidas sobre um elemento representam
a solução em suas derivadas nos nós. Nos métodos clássicos os
coeficientes desconhecidos na solução assumida não tem significado
físico, por exemplo, os polinômios a0, a1,· · ·,an no método de Galer-
kin. Estes coeficientes são apenas quantidades matemáticas que quando
substituídas dentro da solução assumida dão uma solução aproximada.
No FEM, ao contrário, os coeficientes neste polinômio são descritos em
termos das soluções desconhecidas em alguns pontos selecionados do
elemento. Estes pontos são intersecção entre elementos e são chamados
de nós. As variáveis nestes nós são chamadas de graus de liberdade
nos nós. As condições de contorno necessárias dependem do tipo de
problema conforme a tabela (3.1).
Tab. 3.1: Graus de liberdade (DOF) necessário nos nós.
Problema DOFs necessários
segunda-ordem u
quarta-ordem u e u′
sexta-ordem u, u′ e u′′
Nas próximas seções são apresentados as ideias para se escrever elementos
de problemas unidimensionais envolvendo equações diferenciais de segunda
ordem e quarta ordem. O domínio de solução para estes elementos é uma
linha com dois nós.
3.3.1 Funções de interpolação linear para problemas de
segunda ordem
Um elemento linha com dois nós para um problema de segunda ordem
é mostrado na fig. (3.2). O elemento se estende de x1 até x2 com um
comprimento l = x2−x1. Os círculos representam os nós. Para um problema
de segunda ordem os graus de liberdade nos nós são as soluções desconhecidas
nestes nós e são dados simplesmente pelos deslocamento axiais u1 e u2.
35Fig. 3.2: Elemento simples com dois nós para um problema de segunda ordem
[2].
Uma vez que há dois nós em um elemento, uma solução linear pode ser es-
crita em todo o elemento começando com um polinômio com dois coeficientes
e então avaliando estes coeficientes em termos de nós desconhecidos como
segue-se:
u(x) = a0 + a1x (3.49)
x = x1; u(x1) = u1 ⇒ u1 = a0 + a1x1 (3.50)
x = x′2; u(x2) = u2 ⇒ u2 = a0 + a1x2 (3.51)
Resolvendo estas duas equações para a0 e a1 tem-se:
a0 = −u1x2 − u2x1
x1 − x2 (3.52)
a1 = −u2 − u1
x1 − x2 (3.53)
Assim a solução linear em um elemento em termos dos graus de liberdade
nodais é escrita da forma:
u(x) = −u1x2 − u2x1
x1 − x2 −
u2 − u1
x1 − x2x (3.54)
Arranjando a eq. (3.54) de tal forma que seja evidenciado os graus de
liberdade, deslocamento nos nós, os deslocamentos u1 e u2, pode-se escrever:
u(x) =
x− x2
x1 − x2u1 +
x− x1
x2 − x1u2 (3.55)
A eq. (3.55) representa um elemento finito de uma solução linear. Ele é
claramente equivalente a um elemento linear. Entretanto, ao invés de escreve-
lo em função dos coeficientes a0 e a1, que não tem nenhum significado físico,
são escritos em função dos deslocamentos u1 e u2 dos dois nós do elemento.
36
As equações de elementos finitos são geralmente escritas em notação ma-
tricial6. A solução do elemento finito linear visto na eq. (3.55) pode ser
escrita na forma:
u(x) =
2∑
i=1
Ni(x)ui = (N1 N2)
(
u1
u2
)
≡ NTd (3.56)
sendoNi(x) as funções de forma, ui os nós desconhecidos, u(x) a função de
interpolação e d um vetor contendo os nós desconhecidos. Um vetor coluna
2×1 contendo as funções de forma é representado por N. Note que a função
de forma N1 é 1 no nó 1 e 0 no nó 2, enquanto a função de forma N2 é 1 no
nó 2 e 0 no nó 1. Portanto, a função de forma N1 designa a influência de u1
na solução assumida u(x) e N2 a influência de u2. Por esta razão, as funções
de forma são também conhecidas como funções de influência.
3.3.2 Funções de ponderação de Galerkin na forma de
elementos finitos
Em uma forma geral a solução assumida em FEM se expressa com o uso
de funções de forma:
u(x) = (N1 N2 · · ·)
(
u1
u2
)
≡ NTd (3.57)
As funções de ponderação de Galerkin são derivadas das soluções assu-
midas u(x) com respeito aos coeficientes desconhecidos ui:
wi =
du
dui
≡ Ni (3.58)
Assim, as funções de forma Ni(x) são as funções de ponderação quando
se usa o método de aproximação de Galerkin.
3.3.3 Funções de interpolação Hermitianas para um ele-
mento com dois nós
Suponha um elemento linear como mostrado na fig. (3.3). Neste caso o
elemento contém dois nós e em cada nó são dois graus de liberdade (ui, u′i).
Assim ao todo temos quatro graus de liberdade por elemento. Neste caso
37
Fig. 3.3: Elemento simples com dois nós para um problema de quarta ordem
[2].
precisamos de quatro coeficientes na função de forma. Esta interpolação é
conhecida como função de interpolação Hermitiana, e muita usada em FEM.
Para este elemento uma aproximação cubica pode ser assumida:
u(x) = a0 + a1x+ a2x
2 + a3x
3 (3.59)
Derivando uma vez em relação a x obtém-se u′(x):
u′(x) = a1 + 2xa2 + 3x2a3 (3.60)
Tendo as condições de contorno do elemento:
x = 0; ⇒ u1 = a0; u′1 = a1 (3.61)
x = l; ⇒ u2 = a0 + la1 + l2a2 + l3a3; u′2 = a1 + 2la2 + 3l2a3 (3.62)
Com isto sabe-se que a0 = u1 e que a1 = u′1. Os coeficientes a2 e a3 são
dados por:
a2 = −3u1 − 3u2 + 2lu
′
1 + lu
′
2
l2
(3.63)
a3 = −−2u1 + 2u2 − lu
′
1 − lu′2
l3
(3.64)
Aplicando estes resultados na eq. (3.59) e escrevendo em termos dos graus
de liberdade nodais u1, u′1, u2 e u′2:
u(x) = u1 + xu
′
1 −
x3(−2u1 + 2u2 − lu′1 − lu′2)
l3
− x
2(3u1 − 3u2 + 2lu′1 + lu′2)
l2
(3.65)
6Para facilidade de implementação computacional.
38
Escrevendo na forma matricial para evidenciar as funções de forma Ni(x),
i = 1, 2, 3, 4, temos:
u(x) = (N1 N2 N3 N4)

u1
u′1
u2
u′2
 ≡ NTd (3.66)
As quatro funções de forma hermitianas para este elemento linear são
dadas por:
N1(x) = 1 +
2x3
l3
− 3x
2
l2
(3.67)
N2(x) = x+
x3
l2
− 2x
2
l
(3.68)
N3(x) = −2x
3
l3
+
3x2
l2
(3.69)
N4(x) =
x3
l2
− x
2
l
(3.70)
3.4 Solução de elementos finitos de problemas
de deformação axial
O problema da deformação axial de uma barra apresentado pelo PVC da
eq. (3.5) pode ser resolvido de forma simples com um elemento linear com
dois nós com um grau de liberdade por nó, totalizando 2 graus de liberdade
por elemento, u1 e u2, mostrado na fig. (3.4). Este elemento se estende de x1
até x2 e tem um comprimento L = x2 − x1. Cada nó pode ter uma carga Pi
sendo aplicada no nó i. A seguir são apresentados os passos para a solução
deste problema via FEM.
3.4.1 Solução linear assumida
Para este elemento, uma função de interpolação linear u(x) é suficiente,
conforme já discutido anteriormente. Assim a eq. (3.55) pode ser usada
para definir a interpolação do elemento. Evidenciando o comprimento L do
elemento, temos então:
u(x) = (N1 N2)
(
u1
u2
)
≡ NTd (3.71)
39
Fig. 3.4: Elemento linear de barra para o problema de deformação axial [2].
As funções de forma N1(x) e N2(x) são dadas por:
N1(x) =
x− x2
x1 − x2 = −
x− x2
L
(3.72)
N2(x) =
x− x1
x2 − x1 =
x− x1
L
(3.73)
Como será preciso usar a diferenciação da solução aproximada u(x) com
respeito a x quando da utilização do cálculo do resíduo ponderado do método
de Galerkin, devemos ter:
u′(x) =
du(x)
dx
= (N ′1 N
′
2)
(
u1
u2
)
≡ BTd (3.74)
Sendo as derivadas em relação a x das funções de forma Ni(x):
N ′1(x) =
1
x1 − x2 = −
1
L
(3.75)
N ′2(x) =
1
x2 − x1 =
1
L
(3.76)
As derivadas da função de forma podem ser representadas como uma
matriz coluna B.
3.4.2 Equações de elementos usando o método de Ga-
lerkin
A forma fraca do método de Galerkin pode ser usada com os passos já
vistos na seção 3.2 deste capítulo, a saber: escrever os resíduos ponderados
40
pela função wi(x), integração por partes e incorporação das condições de
contorno natural. Como já visto anteriormente, as funções de ponderação
em FEM são as funções de forma Ni(x). Como u(x) é uma solução assumida
o erro é escrito conforme a eq. (3.14), que é repetida a seguir:
e(x) = q +
d
dx
(AEu′) (3.77)
Multiplicando o erro pela função de forma Ni(x)7 e integrando no inter-
valo x1 até x28, tem-se o resíduo no método de Galerkin:∫ x2
x1
e(x)wi(x)dx =
∫ x2
x1
(
qNi +
d
dx
(AEu′)Ni
)
dx = 0 (3.78)
Usando integração por partes conforme a eq. (3.21):
AENi(x2)u
′(x2)− AENi(x1)u′(x1) +
∫ x2
x1
(qNi − AEu′N ′i)dx = 0 (3.79)
As condições de contorno natural no elemento conforme apresentado na
fig. (3.4) são:
−P1 − AEu′(x1) = 0 (3.80)
AEu′(x2)− P2 = 0 (3.81)
Rearranjando estas condições de contorno temos:
u′(x1)→ − P1
AE
(3.82)
u′(x2)→ P2
AE
(3.83)
Assim, aplicando as condições de contorno na forma fraca da eq. (3.79)
teremos:
P1Ni(x1) + P2Ni(x2) +
∫ x2
x1
(qNi − AEu′N ′i)dx = 0 (3.84)
Como temos duas funções de interpolação N1(x) e N2(x) teremos então
duas equações para o elemento a partir da eq. (3.84):
7Lembrando que a função de forma Ni(x) é o resíduo ponderado de Galerkin wi(x).
8Que corresponde ao domínio do elemento.
41
P1N1(x1) + P2N1(x2) +
∫ x2
x1
(qN1(x)− AEu′N ′1(x)) dx = 0 (3.85)
P1N2(x1) + P2N2(x2) +
∫ x2
x1
(qN2(x)− AEu′N ′2(x)) dx = 0 (3.86)
Lembrando que uma função de forma é uma função de influência temos
que N1(x1) = 1, N1(x2) = 0, N2(x1) = 0 e N2(x2) = 1. Escrevendo as
equações em uma forma matricial:
∫ x2
x1
((
N1
N2
)
q −
(
N ′1
N ′2
)
AEu′(x)
)
dx+
(
P1
P2
)
=
(
0
0
)
(3.87)
Uma vez que u′(x) já foi descrito na eq. (3.74) em função da matriz B
contendoas derivadas das funções de forma Ni(x) podemos escrever a eq.
(3.87) em uma forma mais compacta e elegante:∫ x2
x1
Nqdx−
∫ x2
x1
BAEBTddx+ rp = 0 (3.88)
sendo rp um vetor 2 × 1 com as cargas aplicadas nos nós dos elementos
(P1, P2)
T e 0 um vetor 2× 1 de zeros. A eq. (3.88) pode ser escrita na forma
de matriz de rigidez rearranjando-a:∫ x2
x1
BAEBTdxd =
∫ x2
x1
Nqdx+ rp ⇒ kd = rq + rp (3.89)
onde:
k =
∫ x2
x1
BAEBTdx (3.90)
rq =
∫ x2
x1
Nqdx (3.91)
rp =
(
P1
P2
)
(3.92)
d =
(
u1
u2
)
(3.93)
Substituindo as derivadas das funções de forma Ni(x) e conduzindo as
integrações considerando que EA e q são constantes em um elemento teremos
o seguinte resultado para a matriz de rigidez k do elemento local9:
9Lembrando que L = x2 − x1.
42
k =
∫ x2
x1
BAEBTdx =
( ∫ x2
x1
AE 1
L2
dx − ∫ x2
x1
AE 1
L2
dx
− ∫ x2
x1
AE 1
L2
dx
∫ x2
x1
AE 1
L2
dx
)
(3.94)
Resolvendo a integral anterior para o domínio do elemento encontra-se a
matriz de rigidez local k do elemento de barra:
k =
AE
L
(
1 −1
−1 1
)
(3.95)
Já o vetor rq pode ser resolvido por:
rq =
∫ x2
x1
Nqdx =
( ∫ x2
x1
−x−x2
L
qdx∫ x2
x1
x−x1
L
qdx
)
=
qL
2
(
1
1
)
(3.96)
Em resumo, para este problema de deformação axial no elemento temos
que:
kd = rq + rp (3.97)
sendo a matriz de rigidez do elemento dada pela eq. (3.95), o vetor rq
representado pela eq. (3.96) e o vetor rp = (P1 P2)T . A matriz de rigidez
global K pode ser montada a partir da matriz local k tomando cuidado
para escrever a contribuição de cada nó. O exemplo a seguir mostra como
construir a matriz global de rigidez K.
Exemplo 3.1 Considere a barra biengastada, vista na fig. (3.5). Esta barra
tem área da secção transversal A, módulo de elasticidade E e comprimento
total 3L com uma carga uniforme f sendo aplicada. Usando três elementos
lineares de comprimento L cada um calcule a matriz de rigidez global K.
Solução: A viga tem quatro nós, sendo que os nós 1 e 2 estão engas-
tados. A matriz de rigidez local é dada pela eq. (3.95). Os elementos (I) e
(II) dividem o mesmo nó (nó 2), assim como os elementos (II) e (III) que
dividem o nó 3. Aplicando o conceito de conectividade nodal, isto significa
que os termos da matriz de rigidez local k no nó 2 devem ser somados nos
elementos (I e II), o que corresponde a soma existente na linha 2 coluna 2
(que representa o que acontece no nó 2). O mesmo ocorre no nó 3 (linha 3,
coluna 3). Assim:
K =
AE
L

1 −1 0 0
−1 1 + 1 −1 0
0 −1 1 + 1 −1
0 0 −1 1
 (3.98)
43
Fig. 3.5: Barra biengastada discretizada com três elementos [3].
A matriz rigidez global tem ordem 4×4 para esta configuração. As incóg-
nitas neste caso são os deslocamentos axiais nos quatro nós, u1, u2, u3 e u4.
A solução deste conjunto de equações fornece os deslocamento que ocorrem
nestes nós:
Ku = f (3.99)
AE
L

1 0 0 0
0 2 −1 0
0 −1 2 0
0 0 0 1


u1
u2
u3
u4
 = fAL

0
1
1
0
 (3.100)
Deve ser observado que no vetor força f também são somados nos res-
pectivos nós as influências de cada elemento. Note também que a 1.o linha e
coluna e a quarta linha e coluna são escritas de forma diferente em virtude
da aplicação da condição de contorno biengastada, uma vez que u1 e u4 são
nulos. Isto é feito apenas por razões computacionais no momento de inverter
a matriz de rigidez global K visando encontrar o vetor de deslocamento em
toda a barra u.
3.5 Exemplos de problemas unidimensionais
em engenharia
Todo o desenvolvimento apresentado anteriormente, começando com o
método de aproximação de Galerkin até a construção de um elemento finito
linear, passando pelos conceitos de função de interpolação e função de forma,
44
foi realizado até aqui utilizando-se um exemplo clássico de elasticidade para
ilustrar o método. Porém, como já discutido neste capítulo, inúmeros pro-
blemas de engenharia são escritos na forma unidimensional e podemos obter
soluções similares para estes casos. Na sequência apresenta-se as equações
diferenciais básicas de alguns problemas unidimensionais de interesse em apli-
cações de engenharia. Todos estes problemas podem ser solucionados efeti-
vamente usando FEM.
3.5.1 Transferência de calor
Equações que descrevem a condução de calor unidimensional em regime
permanente também são derivadas de equações de balanço e relações consti-
tutivas, assim como problemas de elasticidade. O balanço de energia diz que
uma mudança no fluxo de calor q(x) é balanceada por uma fonte externa de
calor Q(x):
d[q(x)A(x)]
dx
= Q(x)A(x) (3.101)
sendo A(x) uma área. Um valor negativo de Q(x) indica que o calor está
sendo removido do sistema. A relação constitutiva neste caso é conhecida
como lei de Fourier é dada por:
q(x) = −k(x)dT (x)
dx
(3.102)
sendo T (x) o campo de temperatura e k(x) a condutividade térmica.
Combinando as eqs. (3.101) e (3.102) obtém-se a equação diferencial de
segunda ordem que governa a transferência de calor no sistema:
d
dx
(
k(x)A(x)
dT (x)
dx
)
+Q(x)A(x) = 0 (3.103)
É importante notar que uma equação de condução de calor é matemati-
camente idêntica ao problema de deformação axial em uma barra e todo o
desenvolvimento envolvendo o método de aproximação de Galerkin e FEM
pode ser aplicado de forma similar a este exemplo. As condições de contorno
em T são essenciais e podem também ser aplicadas no fluxo de calor q(x).
Um caso particular pode ser encontrado caso a área A, a fonte de calor Q e
a condutividade térmica k sejam assumidas constantes em x, assim:
k
d2T
dx2
+Q = 0 (3.104)
45
3.5.2 Fluxo de potência
O fluxo de potência é uma área especial em mecânica dos fluídos que pode
ser aplicada em problemas de movimento da camada limite. Nesta aplicação
se assume um fluído incompressível, e o problema é completamente descrito
pela equação da continuidade ou balanço de massa. Uma função potencial é
postulada em uma dimensão, assumindo área constante:
φ(x) = −K(x)h(x) = −K(x)
(
z + p
γ
)
(3.105)
e
u(x) =
dφ
dx
(3.106)
sendo u a velocidade do fluído, h é a cabeça piezométrica, z é a cabeça
de elevação, γ é uma ponderação específica em problemas de camada limite,
p é a pressão, e K é um coeficiente de permeabilidade ou condutividade
hidráulica. A equação de condutividade é dada pela lei de Darcy:
u(x) = −Kdh(x)
dx
(3.107)
e segue que a lei de Darcy é relacionada a definição do potencial eq.
(3.105). A equação governante é obtida combinando as eqs. (3.105) e (3.107),
e para um fluído incompressível, du
dx
= 0:
d2φ
dx2
= 0 (3.108)
A solução da equação anterior é uma função linear, e em uma dimensão
será uma constante. Entretanto, a sua contrapartida em duas dimensões é um
desafio e será estudada em detalhes em capítulos seguintes nesta apostila. As
condições de contorno essenciais podem ser aplicadas no fuxo φ e as condições
de contorno naturais na velocidade u(x).
3.5.3 Transferência de massa
A difusão resulta em uma equação similar a vista anteriormente em fluxo
de potência. Neste exemplo, a equação de balanço de massa será descrita tal
que exista uma mistura diluída. Esta teoria é aplicada para uma variedade
de aplicações físicas. Considerando uma área constante, o balanço de massa
pode ser escrito como:
u(x)
dC(x)
dx
+
dj(x)
dx
+KrC(x) = m (3.109)
46
sendo que u(x) é a velocidade da mistura, C(x) é a sua concentração,
j(x) é o fluxo dos componentes da mistura, Kr é a taxa da reação, e m é
uma fonte externa de massa. A equação constitutiva neste caso é a lei de
Fick dada por:
j(x) = −D(x)dC(x)
dx
(3.110)
sendo D(x) a difusidade. Combinando as eqs. (3.109) e (3.110) leva a
equação governante para o fenômeno de transferência de massa:
u(x)dC(x)
dx
− d
dx
[
D(x)
dC(x)
dx
]
+KrC(x) = m (3.111)
A velocidade u(x) é assumida ser conhecida. As condições de contorno
essenciais estão em C(x) e as naturais no fluxo j(x).
3.5.4 Eletricidade
A equações envolvidas em problemas de eletroestática são similares as de
condução de calor. Neste caso, o balanço de carga fornece uma relação entre
o deslocamento elétrico D(x) e a densidade de carga ρ(x) tal que:
d [A(x)D(x)]
dx
= ρ(x)A(x) (3.112)
sendo A(x) a área da secção tranversal perpendicular ao eixo x. O campo
elétrico E(x) é relacionado ao potencial elétrico φ(x) por:
E(x) = −dφ(x)
dx
(3.113)
A relação constitutiva é calculada por:
D(x) = �(x)E(x) = −�(x)dφ(x)
dx
(3.114)
Combinando as eqs. (3.112) e (3.114) dão a seguinte equação:
d
dx
[
�(x)A(x)
dφ(x)
dx
]
+ ρ(x)A(x) = 0 (3.115)
sendo as condições de contorno essenciais aplicadas a φ(x) e as condições
naturais em D(x).
47
3.6 Função no Matlab para resolver um PVC
unidimensional
Todos os problemas de valor de contorno unidimensional vistos neste ca-
pítulo tem a seguinte forma geral:
d
dx
(
k(x)
du(x)
dx
)
+ p(x)u(x) + q(x) = 0; xl < x < xn (3.116)
Por exemplo, no caso do problema elástico k(x) = EA, p(x) = 0 e q(x) =
q. Já para um caso de transferência de calor com parâmetros constantes
k(x) = kA, p(x) = 0 e q(x) = Q, como descrito na eq. (3.104).
Abaixo é transcrita uma função no Matlab contida no livro [2] que pode
ser usada para definir um elemento linear em problemas de valor de contorno
unidimensional. Esta função é geral e pode ser usada em qualquer tipo de
problema. Os parâmetros de entrada na função são: k, p e q além das
coordenadas iniciais e finais do elemento linear, variável coord.
function[ke,re]=BVP1DLinElement(k,p,q,coord)
\% [ke,re]=BVP1DLinElement(k,p,q,coord)
\% Generates equations for a linear element for 1D BVP
\% k,p,q = parameters defining the BVP
\% coord = coordinates at the elements ends
\% ke = local stiffness matrix
\% re = local load
L = coord(2)-coord(1)
ke = [k/L-(L*p)/3, -(k/L)-(L*p)/6; -(k/L)-(L*p)/6, k/L-(L*p)/3];
re = [(L*q)/2;(L*q)/2];
3.7 Exercícios
Ex. 3.1 Para o exemplo da barra uniforme da fig. (3.1) proponha e avalie
uma solução cubica usando o método de Galerkin. Compare graficamente
usando algum software esta solução aproximada com a solução exata e com a
solução aproximada linear e quadrática. Para simulação utilize os parâmetros
E = 70 GPa, A = 600 mm2, c = 200 e L = 600 mm.
Ex. 3.2 Para a barra da fig. (3.6) descreva o problema de valor de contorno
48
resultante10. Escreva a solução exata e usando o método de Galerkin encontre
as soluções aproximadas usando aproximações lineares, quadráticas e cúbicas.
Compare graficamente as respostas usando algum software. Para simulação
utilize os valores E = 1, c = 0, P = 1, A0 = 2, AL = 1 e L = 1.
Fig. 3.6: Barra não-uniforme carregada axialmente [2].
Ex. 3.3 Considere o problema de transferência de calor em uma barra de
comprimento L, com área constante A, condutividade térmica k constante,
e fonte de calor Q também constante. Assuma as seguintes condições de
contorno T (0) = T (L) = 0. Obtenha a solução exata deste problema. Usando
o método de Galerkin obtenha uma solução aproximada empregando uma
função de interpolação quadrática. Compare as soluções encontradas.
Ex. 3.4 Resolva novamente os problemas (3.1) e (3.2), porém, discretizando
as barras em cinco elementos iguais e usando FEM. Calcule a distribuição
de tensões resultantes nas barras para estes dois casos. Compare com as
soluções exatas e avalie a qualidade das respostas. Se possível, após implantar
uma solução em rotina própria utilize algum software comercial de FEM para
realizar a mesma tarefa.
Ex. 3.5 Resolva novamente o exercício (3.3), porém, discretizando a barra
em cinco elementos e usando FEM. Implemente esta solução em algum soft-
ware, como Matlab ou Scilab, e mostre graficamente a distribuição de tempe-
ratura ao longo da barra. Para efeitos de cálculo considere o comprimento da
barra com L = 0.2 m, T (0) = 100o C, a fonte de calor constante é Q = 3×106
W/m3 ao longo da barra, fluxo de calor q = 1.8×106 W/m2 removendo calor
em x = 0.2 m, ou seja, q(0.2) = 1.8 × 106 W/m2, e condutividade térmica
10Cuidado: Para este caso a área não é constante! Uma dica é escrever A(x) =
A0 − A0−ALL x.
49
k = 6000 W/(mK). Assuma que a área da barra é 0.4 × 10−3 m2. Calcule
o campo de temperatura e o fluxo de distribuição de calor em cada nó dos
elementos.11
Ex. 3.6 Uma barra com área da secção transversal constante A e módulo
de elasticidade E é fixada de ambos os lados em suportes rígidos e carre-
gada axialmente com uma força P como ilustrado na fig. (3.7). Encontre o
deslocamento axial e a distribuição de carga usando somente dois elementos
finitos lineares com dois nós. Utilize os valores numéricos: a = 300 mm,
b = 600 mm, E = 200 GPa, A = 200 mm2 e P = 10 kN.
Fig. 3.7: Barra uniforme carregada axialmente [2].
Ex. 3.7 Um coluna no interior de uma construção com quatro pisos é su-
jeita a cargas axiais nos diferentes andares, como mostrado na fig. (3.8).
Determine o deslocamento axial em cada um dos andares. Usando estes des-
locamentos, calcule as distribuições de deformação e tensão na coluna. Uti-
lize os valores numéricos12: altura de cada piso h = 15 ft, Carga nos andares
P1 = 50000 lb, P2 = P3 = 40000 lb, P4 = 35000 lb, módulo de elasticidade
E = 29× 106 lb/in2, área da secção transversal A = 21.8 in2.
Ex. 3.8 Considere a condução de calor através de uma placa com 5 mm
de espessura com uma fonte de calor de 900 W. A área da base é 150 m2
e a condutividade térmica é k = 20 W/m.oC. A fase interna da placa é
sujeita à um fluxo de calor gerado por uma resistência, como mostrado na
fig. (3.9). Durante o regime permanente a temperatura na face externa da
placa é T = 84 oC. Determine a temperatura interna da placa usando FEM?
Este problema pode ser modelado com a seguinte equação:
11Dica: A matriz de "rigidez"local para este exemplo é idêntica a eq. (3.95), porém,
substituindo AE/L por Ak/L e a força f pela fonte de calor Q.
12As unidades estão no sistema inglês.
50
Fig. 3.8: Coluna de um prédio modelada com quatro elementos lineares [2].
kA
(
d2T
dx2
)
= 0; 0 < x < L (3.117)
k = 20 W/m.oC; A =
150
1002
m2; L =
5
1000
m
T (L) = 84; −kAdT (0)
dx
= 900 (3.118)
Utilize um elemento finito linear com dois nós apenas para facilitar13.
Ex. 3.9 Uma placa com 2 metros de altura e 3 metros de largura está a uma
temperatura de 100 oC. Ela é resfriada usando uma aleta de alumínio com
20 cm de largura e 0.3 cm de espessura (com k = 237 W/m.oC). A aleta se
estende por todo o comprimento da placa como mostrado na fig. (3.10). Ao
longo da altura da placa há um espaçamento de 0.4 cm. A temperatura do ar
ambiente é 25 oC e o coeficiente de convecção médio é 30 W/m2.C Determine
13Perceba que no caso geral para usar a função Matlab fornecida k(x) = kA = 310 ,
p(x) = 0, q(x) = 0 e coord = [0 0.005].
51
Fig. 3.9: Placa com fonte de calor [2].
a distribuição de temperatura através da aleta usando FEM e discretizando-a
em 4 elementos lineares. A equação governante deste sistema é:
d
dx
(
kA
dT
dx
)
− hPT + hPT∞ = 0; 0 < x < L (3.119)
P = 2(W + t)/; A = Wt (3.120)
T (0) = 100; kA
dT (L)
dx
= −hA(T − T∞) (3.121)
(3.122)
Os coeficientes na forma geral são:
k(x) = kA; p(x) = −hP ; q(x) = hPT∞ (3.123)
NBC x = L; α = −hA
kA
; β =
hAT∞
kA
(3.124)
Use os seguintes dados numéricos para simulação W = 3 m, t = 0.3/100
m; L = 20/100 m; k = 237; h = 30; T∞ = 30; P = 6.006; A = 0.009;
kA = 2.133; hP = 180.18; hPT∞ = 4504.5; α = −0.126582 e β = 3.16456.
Ex. 3.10 Determine o perfil de velocidades de um fluído escoando entre duas
placas fixas como mostrado na fig.