Analise-Estruturas1 (1)

7/18/2019 Analise-Estruturas1 (1)

http://slidepdf.com/reader/full/analise-estruturas1-1 1/196

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 aoMétodo dos Elementos Finitos (FEM) do curso de graduação em EngenhariaMecâ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 excelenteslivros 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 queesta apostila diminua esta deficiência, porém, sem a pretensão de substituirum 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 aomáximo transcrever o tema de maneira natural com rigor matemático, masnã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 computacionalpara 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 umnível introdutório e acessível aos alunos da graduação e com boa qualidadegráfica1. Espero contar com o apoio dos alunos e demais colaboradores paramelhorar este texto constantemente, visto que ele foi feito de forma rápidae 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 Silvafevereiro de 2009.

1O texto foi redigido com o LATEX 2ε.

2



Sumário

Lista de Figuras 6

Lista de Tabelas 9

1 Introdução 101.1 Exemplos de aplicação . . . . . . . . . . . . . . . . . . . . . . 111.2 Etapas na solução de um problema via FEM . . . . . . . . . . 131.3 Discretização por elementos finitos . . . . . . . . . . . . . . . 13

2 Fundamentos Matemáticos Básicos de FEM 172.1 Análise vetorial . . . . . . . . . . . . . . . . . . . . . . . . . . 17

2.1.1 Produto escalar . . . . . . . . . . . . . . . . . . . . . . 18

2.1.2 Produto vetorial . . . . . . . . . . . . . . . . . . . . . 182.1.3 Gradiente . . . . . . . . . . . . . . . . . . . . . . . . . 192.1.4 Divergente . . . . . . . . . . . . . . . . . . . . . . . . . 192.1.5 Rotacional . . . . . . . . . . . . . . . . . . . . . . . . . 192.1.6 Teorema da divergência . . . . . . . . . . . . . . . . . 202.1.7 Teorema de Green-Gauss . . . . . . . . . . . . . . . . . 20

2.2 Análise matricial . . . . . . . . . . . . . . . . . . . . . . . . . 212.3 Equações diferenciais . . . . . . . . . . . . . . . . . . . . . . . 232.4 Tensores cartesianos . . . . . . . . . . . . . . . . . . . . . . . 242.5 Exercícios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3 Elementos Finitos Unidimensionais 263.1 Solução exata de problemas de deformação axial em barra uni-

forme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263.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 . . . . . . . . . . . . . . . . . . . . . . . 323.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 393.4.1 Solução linear assumida . . . . . . . . . . . . . . . . . 393.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 . . . . . . . . . . . . . . . . . . 453.5.2 Fluxo de potência . . . . . . . . . . . . . . . . . . . . . 463.5.3 Transferência de massa . . . . . . . . . . . . . . . . . . 463.5.4 Eletricidade . . . . . . . . . . . . . . . . . . . . . . . . 47

3.6 Função no Matlab para resolver um PVC unidimensional . . . 483.7 Exercícios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

4 Elementos de Treliça, Viga e Frame 554.1 Treliças planas . . . . . . . . . . . . . . . . . . . . . . . . . . 55

4.1.1 Exemplo de implementação prática . . . . . . . . . . . 604.2 Deformação transversal em vigas . . . . . . . . . . . . . . . . 71

4.2.1 Equação diferencial para deflexão em vigas . . . . . . . 724.2.2 Condições de contorno em vigas . . . . . . . . . . . . . 764.2.3 Exemplo de solução exata em vigas . . . . . . . . . . . 784.2.4 Diagramas de cortante e momento usando o método

das seções . . . . . . . . . . . . . . . . . . . . . . . . . 794.2.5 Elemento de viga com dois nós . . . . . . . . . . . . . 824.2.6 Solução cúbica assumida . . . . . . . . . . . . . . . . . 834.2.7 Matriz de rigidez local em vigas . . . . . . . . . . . . . 844.2.8 Exemplo de aplicação com carregamento concentrado . 874.2.9 Exemplo de aplicação com carregamento distribuído . . 98

4.2.10 Viga com carregamento trapezoidal . . . . . . . . . . . 1054.3 Frames . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

4.3.1 Exemplo de aplicação com carregamento concentrado . 1134.3.2 Exemplo de aplicação com carregamento distribuído . . 122

4.4 Considerações finais . . . . . . . . . . . . . . . . . . . . . . . . 1294.5 Exercícios . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131

4



5 Análise Dinâmica de Estruturas via FEM 142

5.1 Matriz de massa para elementos estruturais . . . . . . . . . . 1425.2 Matriz de massa para barra . . . . . . . . . . . . . . . . . . . 1445.3 Matriz de massa para elemento de treliça plana . . . . . . . . 1455.4 Matriz de massa para elemento de viga . . . . . . . . . . . . . 1465.5 Matriz de massa para elemento de frame . . . . . . . . . . . . 1475.6 Análise modal analítica - resposta livre . . . . . . . . . . . . . 1475.7 Algoritmo de Newmark . . . . . . . . . . . . . . . . . . . . . . 1515.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 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1656.1.1 Teorema da divergência de Gauss . . . . . . . . . . . . 1656.1.2 Teorema de Green-Gauss . . . . . . . . . . . . . . . . . 166

6.2 Equações de elementos finitos usando o método de Galerkin . 1676.3 Elementos finitos retangulares . . . . . . . . . . . . . . . . . . 1726.4 Exemplo: equação de Laplace em um domínio quadrado . . . 1786.5 Elementos finitos triangulares . . . . . . . . . . . . . . . . . . 179

7 Introdução à Mecânica dos Sólidos Computacional 1847.1 Revisão básica de elasticidade . . . . . . . . . . . . . . . . . . 185

7.1.1 O conceito de tensão . . . . . . . . . . . . . . . . . . . 1857.1.2 Tensão principal e direção principal . . . . . . . . . . . 1877.1.3 Critérios de falhas . . . . . . . . . . . . . . . . . . . . 1887.1.4 Relações constitutivas . . . . . . . . . . . . . . . . . . 190

7.2 Equação de elementos finitos em sólidos elásticos . . . . . . . 1927.2.1 Problema no estado plano de tensões . . . . . . . . . . 1957.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 bloco de motor [4]. . . . . . . . . . . . . . 141.2 Tipos e formas de elementos [2]. . . . . . . . . . . . . . . . . . 15

1.3 Diferença entre o contorno físico e a geometria do contorno domodelo 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]. . . . . . . . . . . . . 283.2 Elemento simples com dois nós para um problema de segunda

ordem [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 363.3 Elemento simples com dois nós para um problema de quarta

ordem [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 383.4 Elemento linear de barra para o problema de deformação axial

[2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403.5 Barra biengastada discretizada com três elementos [3]. . . . . 443.6 Barra não-uniforme carregada axialmente [2]. . . . . . . . . . . 493.7 Barra uniforme carregada axialmente [2]. . . . . . . . . . . . . 503.8 Coluna de um prédio modelada com quatro elementos lineares

[2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.9 Placa com fonte de calor [2]. . . . . . . . . . . . . . . . . . . . 523.10 Transferência de calor através de uma aleta [2]. . . . . . . . . 533.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 noplano [2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

4.2 Treliça plana a ser analisada [4]. . . . . . . . . . . . . . . . . . 604.3 Visualização da malha implementada com o Matlab R. . . . . 634.4 Flexão em vigas [2]. . . . . . . . . . . . . . . . . . . . . . . . . 734.5 Forças e momentos agindo em um elemento diferencial da viga

[2]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 734.6 Condições de contorno comuns em vigas [2]. . . . . . . . . . . 764.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. . . . . . . . . . . . . . . . 814.10 Elemento de viga com flexão. . . . . . . . . . . . . . . . . . . 824.11 Carregamento distribuído concentrado nos nós. . . . . . . . . . 854.12 Viga com carregamento concentrado aplicado. . . . . . . . . . 884.13 Deslocamento transversal da viga solucionado via FEM com

funções de interpolação Hermitiana. . . . . . . . . . . . . . . . 964.14 Momento fletor M (x) calculado usando FEM. . . . . . . . . . 974.15 Cortante V (x) calculada usando FEM. . . . . . . . . . . . . . 974.16 Tensão normal σ(x) devido a flexão calculada usando FEM. . 984.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ídocalculado 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 deforma trapezoidal em um elemento. . . . . . . . . . . . . . . . 107

4.22 Elemento de frame no plano. . . . . . . . . . . . . . . . . . . . 109

4.23 Frame no plano com carregamento concentrado. . . . . . . . . 1134.24 Malha do frame. . . . . . . . . . . . . . . . . . . . . . . . . . . 1144.25 Frame no plano com carregamento distribuído. . . . . . . . . . 1234.26 Treliça 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1314.27 Treliça 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1324.28 Treliça 3. Considere d = 2.0 m e L = 1000 N. . . . . . . . . . 1324.29 Treliça 4. Considere a = 3.0 m . . . . . . . . . . . . . . . . . . 1324.30 Treliça 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1334.31 Treliça 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1344.32 Treliça 7. α = 30o e P = 500 N . . . . . . . . . . . . . . . . . . 1344.33 Treliça 8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1354.34 Viga 1 com w0 = 10 kN/m, L = 5 m, a = 2 m, b = 4 m e

M = 20 kNm. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1354.35 Viga 2 com L = 4 m e w0 = 10 kN/m. . . . . . . . . . . . . . 1354.36 Viga 3 com a = 1 m, L = 4 m e w0 = 10 kN/m. . . . . . . . . 1364.37 Viga 4 com a = 2 m, b = 3 m, L = 4 m, P = 10 kN e w0 = 10

kN/m. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1364.38 Viga 5 com a = 1 m e w0 = 10 kN/m. . . . . . . . . . . . . . . 1374.39 Viga 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1374.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−2m2 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. . . . . . . . . . . . . . . . . . . . . . . 1394.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. . . . . . . . . . . . . . 1394.45 Frame 4 com barras quadradas de 200 mm × 200 mm, E = 10

GPa e P = 20 kN. . . . . . . . . . . . . . . . . . . . . . . . . 1404.46 Frame 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1404.47 Frame 6. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141

4.48 Frame 7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1414.49 Barragem de concreto. . . . . . . . . . . . . . . . . . . . . . . 141

5.1 Esquema de aceleração média constante de Newmark. . . . . . 1535.2 Resposta para o sistema. . . . . . . . . . . . . . . . . . . . . . 1605.3 Resposta para o sistema. . . . . . . . . . . . . . . . . . . . . . 1615.4 1.o e 2.o Modo de vibrar da viga calculado analiticamente. . . 1625.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. . . . . . . . . . . . . 1656.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. . . . . . . . . . 353.2 Viscosidade µ do fluído. . . . . . . . . . . . . . . . . . . . . . 53

4.1 Deslocamentos nodais (mm). . . . . . . . . . . . . . . . . . . . 704.2 Forças externas aplicadas nos nós (kgf). . . . . . . . . . . . . . 704.3 Deformações , tensões σ e forças internas axiais F nas barras. 71774.5 Deslocamentos transversais e rotações na viga. . . . . . . . . . 964.6 Esforços atuantes nos nós do frame. . . . . . . . . . . . . . . . 1224.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 ede 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 formaexata só é possível para sistemas muito simples. Assim, para sistemas maiscomplexos envolvendo geometrias e condições de contorno mais sofisticadasnão é possível se obter uma solução exata. Nestes casos deve-se optar porprocedimentos de aproximação com precisão aceitável para a aplicação deengenharia 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 decontorno, 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 tornandoo mais popular de todos, sobretudo pelo aparecimento de diferentes pacotesde 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 paradiscreto. Esta divisão do domínio é conhecida como malha ou grid, que nadamais é do que o conjunto de elementos finitos resultante da discretização. Amalha é 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 paratodo o domínio, o que pode ser praticamente impossível em um problemacomplexo, e sim buscar estas soluções em cada elemento separadamente.Suponha que o funcional para um elemento seja Ψi, sua soma sobre a malhacom 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 ordem m descrita em função dos nós dos elementos (parâmetrosnodais α j) e por funções de forma (φ). A função de interpolação é descritacomo:

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 valoresdos parâmetros nodais α j . Os parâmetros nodais podem estar associados adeslocamentos, forças internas, tensões, temperaturas, pressão, etc. e dep-dende da formulação do elemento usado. Todos estes pontos serão melhordiscutidos no decorrer deste curso. Porém de antemão o que se espera é que oaluno compreenda que o FEM é uma busca por uma solução local que possaser generalizada para todo o domínio. Várias abordagens do método FEMsão usadas.

1.1 Exemplos de aplicação

O FEM têm inúmeras aplicações nos diferentes ramos da ciência, emespecial em aplicações estruturais. Historicamente, as primeiras utilizaçõesde FEM em engenharia foram em aplicações aeronáuticas e de estruturascivis, 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 projetosde aeronaves sem o uso consistente de ferramentas envolvendo elementosfinitos.

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 deseu 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 umproblema via FEM sem o computador não saberá como proceder tendo umamá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 deprojetos basta "decorar"meia dúzia de comandos para se dizer especialistaem FEM. Porém, isto é um conceito errado. O autor do livro [4] cita umexemplo interessante: Imagine que você está muito doente e procura um mé-dico que não é um grande especialista na sua enfermidade. O médico dizpara não se preocupar, pois ele tem um programa onde "basta"digitar naentrada os sintomas que ele fornece na saída os diagnósticos com a profilaxiaadequada. 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 umaaeronave 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ásicosdo método. Isto permite que ele use pacotes comerciais com maior rigor deaná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 softwarese nunca estudou FEM de forma convencional não se pode dizer que saiba oque é 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 geometriacomplexa. A meta foi calcular o panorama de tensões, que é mostrado namesma 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 emprojetos 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 elementosfinitos.

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 comcalma. 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 elementosde 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 exemplosdestes elementos. Os programas comerciais de FEM possuem bibliotecas comcentenas 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 elementodeve 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 cursopretende-se estudar alguns elementos básicos e clássicos como os elementosde 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 envolvefunçõ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 elementosfinitos no domínio global conduz em um sistema de equações algébricas degrande 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 paraencontrar soluções de equações diferenciais.

A qualidade (acurácia) da aproximação é diretamente proporcional aquantidade de elementos usados. O custo computacional também é ligadoao 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 podemusar malhas mais grosseiras. Os contornos curvilíneos são exemplos ondemalhas finas devem ser usadas, conforme a fig. 1.3.

A escolha de um elemento, numero de elementos, etc. deve ser pautada notipo 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 obtidoe depois refinar a malhar conforme desejo, obtendo assim uma economiacomputacional 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 aconectividade dos nós. Uma malha deve ser formada por elementos que seconectam 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. Afig. (1.4) apresenta um exemplo de malha válida e inválida para um elementocom quatro nós. A malha fica inválida quando usando uma configuração comtrês elementos, uma vez que o nó 4 dos elementos (2) e (3) não está conectadocom 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 formasimé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ênticosde treliças. Porém deve-se ter cuidado, pois se a estrutura for simétrica aalguma linha de centro, a escolha de um elemento não adequado pode fazercom 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 procedimentosmatemáticos que vão dos mais simples, como manipulação de matrizes, até osmais avançados, envolvendo por exemplo teoremas de cálculo vetorial. Estecapítulo tem como meta revisar alguns destes conceitos que serão utilizadosno decorrer do curso. Não será dado nenhum rigor matemático, muito menosprovas dos teoremas mostrados. O aluno interessado pode consultar livrosbá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çãode vetores unitários (i, j, k):

a = axi + ay j + 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çãovetorial, o que forma um campo vetorial. Um campo vetorial pode ser umdeslocamento, 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 emum 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:

∇ ≡ ∂

∂xi +

∂

∂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

∂ 2xi +

∂ 2

∂ 2y j +

∂ 2

∂ 2z k (2.3)

O operador ∇ é usado para definir três operações básicas envolvendocampos escalares e vetoriais:

• Gradiente

• Divergente

• Rotacional

Estas operações são usadas na definição de teoremas fundamentais deintegrais de vetores tais como o Teorema da Divergência e o Teorema deGreen-Gauss. Estes dois teoremas são a base matemática para compreendero método de Galerkin, que por sua vez é uma das bases fundamentais de

FEM. Antes de aprofundar nestas questões é interessante revisar operaçõesbá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 daoperaçã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ãoortonormais (θ = 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 quei × j = k. Importante observar que i × j = j × i.

2.1.3 Gradiente

O gradiente de uma função escalar Φ(x,y,z )1 é dado por:

∇Φ =

∂

∂xi +

∂

∂y j +

∂

∂z k

Φ =

∂ Φ

∂xi +

∂ Φ

∂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 deuma grandeza escalar por unidade de espaço. O significado físico pode serinterpretado como a diferença de temperatura nas faces de um bloco, paraeste tipo de aplicação.

2.1.4 Divergente

O divergente já é uma operação envolvendo um campo vetorial dado poruma função vetorial do tipo a(x,y,z )2e calculado por:

∇ · a = ∂

∂x i +∂

∂y j +∂

∂z k · (axi + ay j + azk) (2.7)

o que leva a seguinte expressão:

∇ · a =∂ax

∂x+

∂ay

∂y+

∂az

∂z (2.8)

note que ∇ · a = 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. noespaç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 ) quevisa observar suas características nestas novas coordenadas. A representaçãoacima 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 vetor ortonormalà esta superfície S 3. O teorema da divergência relaciona o divergente totalde um campo vetorial a em um volume V com o fluxo total deste campovetorial 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 umvalor 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 ba

d

dx

dΦ

dxΨ

dx = k ba

d2Φ

dx2Ψdx + k

ba

dΦ

dx

dΨ

dxdx (2.12)

Notanto que o lado esquerdo da eq. (2.12) forma uma integral perfeitatem-se que:

3a · n pode representar um fluxo.

4O próximo capítulo irá revisar alguns destes problemas.

20



k ba

ddΦdx

Ψ = k dΦdx

Ψba

(2.13)

Substituindo a eq. (2.13) em (2.12) e rearranjando tem-se:

k

ba

d2Φ

dx2Ψdx = k

dΦ

dxΨ

ba

− k

ba

dΦ

dx

dΨ

dxdx (2.14)

Considerando que a = β b o teorema da divergência, eq. (2.10), pode serreescrito 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. EmFEM a eq. (2.14) é uma extensão da eq. (2.16) sendo que Φ e Ψ sãomatrizes representando funções de interpolação (funções aproximadores) doselementos empregados em uma discretização. Em um momento oportunonos 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çãode 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 dematrizes 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 decontorno ou restrição no sistema em estudo. Suponha uma matriz A dadapor:

A =

a11 a12 a13a21 a22 a23a31 a32 a33

(2.17)

Se uma restrição for imposta de tal forma que a segunda linha e colunasejam eliminadas temos uma matriz M22 dada por:

21



M22 = a11 a13a31 a32

(2.18)

Este conceito também é usado para cálculo do cofator C ij:

C ij = (−1)i+ j |M ij| (2.19)

Já o adjunto de uma matriz é Aij = C T ij .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 deliberdade em cada nó de um elemento (por exemplo, deslocamento), a matrizA 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 é feitaa 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 −81 2 1 0

0.5 −1 4 4−4 −2 0 1

x1

x2

x3

x4

=

42

100

(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 −20 1.5 1.5 20 −1.25 4.25 50 0 −2 −7

x1

x2

x3

x4

=

11

9.54

(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 −20 1 1 1.3333

0 0 5.5 6.66670 0 −2 −7

x1

x2

x3x4

=

10.6667

10.33334

(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 −20 1 1 1.33330 0 1 1.21220 0 0 −4.578

x1

x2

x3

x4

=

1

0.66671.87887.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çãopara solucionar de forma numérica e com aproximações uma equação diferen-cial. Sendo assim, é primordial que o engenheiro saiba modelar fisicamenteo seu problema com o conhecimento necessário para construir este sistemade 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 elementossimples, como barra, viga e placa. Estes problemas podem ser descritos porequações diferenciais parciais (problema de valor de contorno). Nestes exem-plos, o FEM pode ser formulado a partir de métodos de energia envolvendofunções de variação (método de Lagrange) sem necessariamente considerarequações diferenciais. Esta representa uma das formas de abordagem deFEM. Porém, neste texto inicialmente será dada uma abordagem da for-mulação de FEM diretamente nas equações diferenciais juntamente com ametodologia 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 queo 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 tipode 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 dosprincipais 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çãode tensor. Tensor é uma grandeza que precisa de nove elementos para poderser completamente conhecida. Em alguns casos com 6 elementos é possíveldescrever 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 umanotação vetorial. Um vetor descrito nesta notação é um tensor de primeiraordem.

Imagine um vetor f escrito em função do sistema de coordenadas (x,y,z ):

f = f xi + f y j + f zk (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 = f 1i + f 2 j + f 3k (2.29)

5Estas equações também podem ser escritas de forma parcial quando envolvem maisde duas variáveis

6Também conhecida como notação de Einstein ou notação tensorial.

24



Em uma notação tensorial este vetor pode ser dado por:

f = f i, i = 1, 2, 3 (2.30)

2.5 Exercícios

Ex. 2.1 Considere um campo vetorial dado por a = (−y,xy,z ). Calcule odivergente 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çãode Gauss.

2 1 2 −32 −2 1 −41 0 2 −34 4 −4 1

x1

x2

x3

x4

=

05

−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ãodefinidos 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 decontorno (PVC) descrito por uma equação diferencial. Equações que go-vernam fenômenos de interesse em engenharia são geralmente obtidas porequações de balanço e equações constitutivas. Neste capítulo todas as basesmatemáticas para FEM visando resolver este tipo de problema são apresen-tadas usando casos unidimensionais. O caso a ser estudado em detalhes nestecapítulo é o problema de elasticidade em barras sujeitas a cargas axiais. Ocapítulo começa apresentando a solução exata deste problema. Na sequênciao método de aproximação de Galerkin é mostrado em detalhes preparandopara 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 engenhariasão apresentados, como o problema de transferência de calor e massa, fluxode potência e eletricidade. Extensões para estes problemas podem ser pro-postas a partir dos conceitos exemplificados ao longo deste capítulo. Sendoassim, no final do capítulo exercícios são propostos ao estudante envolvendodiferentes aplicações de FEM em problemas unidimensionais de interesse emengenharia.

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 exatauma vez conhecidas as condições de contorno. Nesta seção é formulado oproblema básico que será usado como exemplo benchmark no decorrer destecapí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 umadeformação linear em termos de área A, tensão normal σ e força axial f . Obalanç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 aLei 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. Asua carga é linearmente variante q (x) = cx, sendo c uma constante. Esta es-trutura é descrita pelo seguinte problema de valor de contorno (PVC) usandoas eqs. (3.1) e (3.4)1:

EAd2u

dx2+ cx = 0; 0 < x < L (3.5)

u(0) = 0; EAdu(L)

dx= P

As condições de contorno em u(x) são essenciais ou geométricas e ascondições de contorno em σ são naturais.

A solução exata para este problema pode ser facilmente obtida integrandoa eq. (3.5) duas vezes e usando as condições de contorno para avaliar asconstantes 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





exata é calculada facilmente. Infelizmente, nem todos os problemas reais de

engenharia podem ser descritos de forma tão simples assim. Com o propósitode propor uma solução aproximada a próxima seção apresenta o método deGalerkin.

3.2 Aproximação via método de Galerkin

No método de Galerkin assume-se uma forma geral de solução para umPVC. 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çãopode ter qualquer forma. No geral, são mais usadas soluções aproximadasem formas polinomiais do tipo:

u(x) = a0 + a1x + a2x2 + · · · + anxn (3.13)

sendo a0, a1,. . .,an parâmetros desconhecidos. É importante observarque uma vez que u(x) é uma aproximação, esta solução não irá satisfazer aequação diferencial para todos os valores de x. Substituindo está solução noPVC da deformação axial da barra obtém-se um erro e(x):

e(x) =d

dx AE d

u

dx+ q (x) = 0 (3.14)

O erro total, chamado de resíduo, para o domínio inteiro da solução podeser obtido integrando o erro e(x) sobre todo o domínio. Entretanto, umavez 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 serdado 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 formafraca. A equação diferencial é a forma forte pois exige que o erro e(x) sejanulo em todo o domínio 0 < x < L.

O método mais popular nas aplicações de elementos finitos é o método deGalerkin. 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 aosparâmentros desconhecidos ai

2:

2No caso descrito a seguir irá ser considerado a derivada parcial, uma vez que em umcaso 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çõespara 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étodode 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):

xlx0

e(x)wi(x)dx =

xlx0

d

dx

AE

dudx

+ q (x)

wi(x)dx = 0; (3.18)

i = 0, 1, . . . , n

Vamos assumir que A e E possam ser funções. Escrevendo as integraisseparadas:

xlx0

d

dx

AE

dudx

wi(x)dx +

xlx0

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 segundaordem em u, sendo assim deve-se aplicar integração por partes. Lembrandoque se tivermos duas funções f (x) e g(x), a integração por partes é definidacomo:

xlx0

ddx

(f (x)) g(x)dx = f (xl)g(xl) − f (x0)g(x0) − xlx0

ddx

(g(x)) f (x)dx

(3.20)Aplicando integração por partes na eq. (3.19) fornece:

A(xl)E (xl)du(xl)

dxwi(xl) − A(x0)E (x0)

du(x0)

dxwi(x0) (3.21)

− xlx0

AE dudx

dwi

dxdx +

xlx0

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 P x0 é aplicada emx0 tem-se que:

A(xl)E (xl)du(xl)

dxwi(xl) = P x0wi(x0) (3.22)

e se uma força P xl é aplicada na extremidade l tem-se:

A(x0)E (x0)du(x0)

dxwi(x0) = P xlwi(x0) (3.23)

Aplicando as eqs. (3.22) e (3.23) na eq. (3.21) a forma fraca pode ser

escrita por: xlx0

E dudx

dwi

dxAdx =

xlx0

q (x)wi(x)dx + P xlwi(xl) + P x0wi(x0) (3.24)

Para problemas estruturais a eq. (3.24) pode ser interpretada como obem conhecido princípio do trabalho virtual 3. Se a função de ponderaçãowi(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 deslocamentovirtual na barra, o trabalho externo virtual é igual ao trabalho interno virtual.Este princípio é altamente usado em mecânica estrutural, seja em aplicaçõesestáticas ou dinâmica. Sua ampla aplicação é uma das principais razões quefazem com que o método de Galerkin seja o mais popular no desenvolvimentodas equações de elementos finitos.

3.2.1 Aplicação do método de Galerkin na solução doproblema da barra axial

Após apresentado os passos básicos do método de Galerkin, está seçãoapresenta sua aplicação para a solução do PVC descrito pela eq. (3.5). Váriassoluções aproximadas podem ser propostas. Uma vez que consideramos queEA é constante, q (x) = cx, domínio 0 < x < L e as condições de contornosão u(0) = 04 e EAu(L) = P , a forma fraca específica deste problema é dadapor:

3Que dever ter sido estudado no curso básico de mecânica geral.4E consequentemente w(0) = 0.

31



P wi(L) + L0

−AE dudx

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 dedeslocamento 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:

P wi(L) +

L0

−AE (a1)

dwi

dx+ cxwi(x)

dx = 0 (3.29)

Existe somente um parâmetro desconhecido, a1. A função de ponderaçãode 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:

P L + L

0 −AE (a1) + cx

2 dx = 0 (3.31)Resolvendo esta integração:

cL3

3+ P L − 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 = −−cL2 − 3P 3EA

(3.33)

Assim a solução linear aproximada é dada pela eq. (3.28) considerando oresultado da eq. (3.33):

u(x) = a1x =cL2 + 3P

3EAx (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) = a2x2 + 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) = a2x2 + a1x ⇒ du

dx= 2a2x + a1 (3.37)

Substituindo o resultado da eq. (3.37) na forma fraca, eq. (3.25) tem-se:

P wi(L) +

L0

−AE (2a2x + a1)

dwi

dx+ cxwi(x)

dx = 0 (3.38)

Agora existem dois parâmetros desconhecidos, a1 e a2, e necessitaremosde duas equações para encontrá-los. Estas equações são obtidas a partir dasfunçõ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) = L2 (3.40)

Substituindo w1 na forma fraca eq. (3.38)tem-se:

P L +

L0

−AE (2a2x + a1) + cx2

dx = 0 (3.41)

⇒ cL3

3− EAa2L2 + P L − EAa1L = 0 (3.42)

33



e substituindo w2 na eq. (3.38) tem-se:

P L2 +

L0

−AE (2a2x + a1) (2x) + cx3

dx = 0 (3.43)

⇒ cL4

4− 4

3EAa2L3 + P L2 − 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 = −−7cL2

− 12P 12EA (3.45)

a2 = − cL

4EA(3.46)

Por fim, a solução quadrática aproximada é:

u(x) = a2x2 + a1x =(12P + cL(7L − 3x))

12EAx (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 guiaresta escolha, especialmente para problemas de duas ou três dimensões. Alémdisto, 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çãode 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 simpleschamados de elementos . Cada elemento tem uma geometria simplestal 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 descrevera 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: l0

(· · ·) =

x2x1=0

(· · ·) dx +

x3x2

(· · ·) 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 oscoeficientes desconhecidos na solução assumida não tem significadofí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 emtermos das soluções desconhecidas em alguns pontos selecionados doelemento. Estes pontos são intersecção entre elementos e são chamadosde 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 deproblema conforme a tabela (3.1).

Tab. 3.1: Graus de liberdade (DOF) necessário nos nós.

Problema DOFs necessáriossegunda-ordem uquarta-ordem u e u

sexta-ordem u, u e u

Nas próximas seções são apresentados as ideias para se escrever elementosde problemas unidimensionais envolvendo equações diferenciais de segundaordem e quarta ordem. O domínio de solução para estes elementos é umalinha com dois nós.

3.3.1 Funções de interpolação linear para problemas desegunda 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 umcomprimento l = x2−x1. Os círculos representam os nós. Para um problemade segunda ordem os graus de liberdade nos nós são as soluções desconhecidasnestes nós e são dados simplesmente pelos deslocamento axiais u1 e u2.

35



Fig. 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 coeficientese então avaliando estes coeficientes em termos de nós desconhecidos comosegue-se:

u(x) = a0 + a1x (3.49)

x = x1; u(x1) = u1 ⇒ u1 = a0 + a1x1 (3.50)

x = x2; 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 liberdadenodais é escrita da forma:

u(x) = −u1x2 − u2x1

x1

−x2

− u2 − u1

x1

−x2

x (3.54)

Arranjando a eq. (3.54) de tal forma que seja evidenciado os graus deliberdade, deslocamento nos nós, os deslocamentos u1 e u2, pode-se escrever:

u(x) =x − x2

x1 − x2

u1 +x − x1

x2 − x1

u2 (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 serescrita na forma:

u(x) =2

i=1

N i(x)ui = (N 1 N 2)

u1

u2

≡ NT d (3.56)

sendo N i(x) as funções de forma, ui os nós desconhecidos, u(x) a função deinterpolação e d um vetor contendo os nós desconhecidos. Um vetor coluna2 × 1 contendo as funções de forma é representado por N. Note que a função

de forma N 1 é 1 no nó 1 e 0 no nó 2, enquanto a função de forma N 2 é 1 nonó 2 e 0 no nó 1. Portanto, a função de forma N 1 designa a influência de u1

na solução assumida u(x) e N 2 a influência de u2. Por esta razão, as funçõesde forma são também conhecidas como funções de influência .

3.3.2 Funções de ponderação de Galerkin na forma deelementos finitos

Em uma forma geral a solução assumida em FEM se expressa com o usode funções de forma:

u(x) = (N 1 N 2 · · ·)

u1

u2

≡ NT d (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

≡ N i (3.58)

Assim, as funções de forma N i(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 oelemento contém dois nós e em cada nó são dois graus de liberdade (ui, ui).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 + a2x2 + a3x3 (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; u1 = a1 (3.61)

x = l; ⇒ u2 = a0 + la1 + l2a2 + l3a3; u2 = a1 + 2la2 + 3l2a3 (3.62)

Com isto sabe-se que a0 = u1 e que a1 = u1. Os coeficientes a2 e a3 sãodados por:

a2 = −3u1 − 3u2 + 2lu1 + lu2l2

(3.63)

a3 = −−2u1 + 2u2

−lu1

−lu2

l3 (3.64)Aplicando estes resultados na eq. (3.59) e escrevendo em termos dos graus

de liberdade nodais u1, u1, u2 e u2:

u(x) = u1 + xu1 − x3(−2u1 + 2u2 − lu1 − lu2)

l3− x2(3u1 − 3u2 + 2lu1 + lu2)

l2(3.65)

6Para facilidade de implementação computacional.

38



Escrevendo na forma matricial para evidenciar as funções de forma N i(x),

i = 1, 2, 3, 4, temos:

u(x) = (N 1 N 2 N 3 N 4)

u1

u1u2

u2

≡ NT d (3.66)

As quatro funções de forma hermitianas para este elemento linear sãodadas por:

N 1(x) = 1 + 2x

3

l3 − 3x

2

l2 (3.67)

N 2(x) = x +x3

l2− 2x2

l(3.68)

N 3(x) = −2x3

l3+

3x2

l2(3.69)

N 4(x) =x3

l2− x2

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 daeq. (3.5) pode ser resolvido de forma simples com um elemento linear comdois nós com um grau de liberdade por nó, totalizando 2 graus de liberdadepor 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 P isendo aplicada no nó i. A seguir são apresentados os passos para a soluçãodeste 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 usadapara definir a interpolação do elemento. Evidenciando o comprimento L doelemento, temos então:

u(x) = (N 1 N 2)

u1

u2

≡ NT d (3.71)

39



Fig. 3.4: Elemento linear de barra para o problema de deformação axial [2].

As funções de forma N 1(x) e N 2(x) são dadas por:

N 1(x) =x − x2

x1 − x2

= −x − x2

L(3.72)

N 2(x) =x − x1

x2 − x1

=x − x1

L(3.73)

Como será preciso usar a diferenciação da solução aproximada u(x) comrespeito 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

≡ BT d (3.74)

Sendo as derivadas em relação a x das funções de forma N i(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 umamatriz 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çãoem FEM são as funções de forma N i(x). Como u(x) é uma solução assumidao 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 N i(x)7 e integrando no inter-valo x1 até x2

8, tem-se o resíduo no método de Galerkin:

x2

x1

e(x)wi(x)dx = x2

x1 qN i +d

dx

(AEu)N i dx = 0 (3.78)

Usando integração por partes conforme a eq. (3.21):

AEN i(x2)u(x2) − AEN i(x1)u(x1) +

x2x1

(qN i − AEuN i)dx = 0 (3.79)

As condições de contorno natural no elemento conforme apresentado nafig. (3.4) são:

−P 1 − AEu

(x1) = 0 (3.80)AEu(x2) − P 2 = 0 (3.81)

Rearranjando estas condições de contorno temos:

u(x1) → − P 1AE

(3.82)

u(x2) → P 2AE

(3.83)

Assim, aplicando as condições de contorno na forma fraca da eq. (3.79)teremos:

P 1N i(x1) + P 2N i(x2) +

x2x1

(qN i − AEuN i)dx = 0 (3.84)

Como temos duas funções de interpolação N 1(x) e N 2(x) teremos entãoduas equações para o elemento a partir da eq. (3.84):

7Lembrando que a função de forma N i(x) é o resíduo ponderado de Galerkin wi(x).8Que corresponde ao domínio do elemento.

41



P 1N 1(x1) + P 2N 1(x2) + x2x1

(qN 1(x) − AEuN 1(x)) dx = 0 (3.85)

P 1N 2(x1) + P 2N 2(x2) +

x2x1

(qN 2(x) − AEuN 2(x)) dx = 0 (3.86)

Lembrando que uma função de forma é uma função de influência temosque N 1(x1) = 1, N 1(x2) = 0, N 2(x1) = 0 e N 2(x2) = 1. Escrevendo asequações em uma forma matricial:

x2x1 N 1

N 2 q − N 1N 2 AEu

(x) dx + P 1P 2 = 0

0 (3.87)

Uma vez que u(x) já foi descrito na eq. (3.74) em função da matriz B

contendo as derivadas das funções de forma N i(x) podemos escrever a eq.(3.87) em uma forma mais compacta e elegante: x2

x1

Nqdx − x2x1

BAE BT ddx + r p = 0 (3.88)

sendo r p um vetor 2 × 1 com as cargas aplicadas nos nós dos elementos(P 1, P 2)T e 0 um vetor 2 × 1 de zeros. A eq. (3.88) pode ser escrita na forma

de matriz de rigidez rearranjando-a: x2x1

BAE BT dxd =

x2x1

Nqdx + r p ⇒ kd = rq + r p (3.89)

onde:

k =

x2x1

BAE BT dx (3.90)

rq =

x2x1

Nqdx (3.91)

r p = P 1P 2 (3.92)

d =

u1

u2

(3.93)

Substituindo as derivadas das funções de forma N i(x) e conduzindo asintegrações considerando que EA e q são constantes em um elemento teremoso seguinte resultado para a matriz de rigidez k do elemento local9:

9Lembrando que L = x2 − x1.

42



k = x2x1

BAE BT dx = x2

x1AE 1

L2dx − x2

x1AE 1

L2dx

− x2x1

AE 1L2

dx x2x1

AE 1L2

dx

(3.94)

Resolvendo a integral anterior para o domínio do elemento encontra-se amatriz 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 = x2x1

Nqdx = x2

x1−x−x2

L qdx x2x1

x−x1L

qdx

=

qL

2

11

(3.96)

Em resumo, para este problema de deformação axial no elemento temosque:

kd = rq + r p (3.97)

sendo a matriz de rigidez do elemento dada pela eq. (3.95), o vetor rqrepresentado pela eq. (3.96) e o vetor r p = (P 1 P 2)T . A matriz de rigidezglobal 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 comoconstruir 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 comprimentototal 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 00 −1 1 + 1 −10 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 00 2 −1 00 −1 2 00 0 0 1

u1u2

u3

u4

= f AL 0

110

(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ãonulos. 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 ométodo de aproximação de Galerkin até a construção de um elemento finitolinear, 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 obtersoluções similares para estes casos. Na sequência apresenta-se as equaçõesdiferenciais 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 queuma mudança no fluxo de calor q (x) é balanceada por uma fonte externa decalor 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 é conhecidacomo 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 desegunda 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 odesenvolvimento envolvendo o método de aproximação de Galerkin e FEMpode ser aplicado de forma similar a este exemplo. As condições de contornoem 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 ea condutividade térmica k sejam assumidas constantes em x, assim:

kd2T

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 podeser aplicada em problemas de movimento da camada limite. Nesta aplicaçãose assume um fluído incompressível, e o problema é completamente descritopela 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çade elevação, γ é uma ponderação específica em problemas de camada limite, p é a pressão, e K é um coeficiente de permeabilidade ou condutividadehidráulica. A equação de condutividade é dada pela lei de Darcy:

u(x) = −K dh(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ãoserá uma constante. Entretanto, a sua contrapartida em duas dimensões é umdesafio e será estudada em detalhes em capítulos seguintes nesta apostila. Ascondições de contorno essenciais podem ser aplicadas no fuxo φ e as condiçõesde 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 fluxode potência. Neste exemplo, a equação de balanço de massa será descrita talque exista uma mistura diluída. Esta teoria é aplicada para uma variedadede aplicações físicas. Considerando uma área constante, o balanço de massapode ser escrito como:

u(x)dC (x)

dx+

dj(x)

dx+ K rC (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, K r é a taxa da reação, e m éuma fonte externa de massa. A equação constitutiva neste caso é a lei deFick 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 aequação governante para o fenômeno de transferência de massa:

u(x)dC (x)

dx −

d

dx D(x)dC (x)

dx + K rC (x) = m (3.111)

A velocidade u(x) é assumida ser conhecida. As condições de contornoessenciais 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 decondução de calor. Neste caso, o balanço de carga fornece uma relação entreo 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 campoelé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çõesnaturais 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 constantesk(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 contornounidimensional. Esta função é geral e pode ser usada em qualquer tipo deproblema. Os parâmetros de entrada na função são: k, p e q além dascoordenadas 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íciosEx. 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 mm 2, c = 200 e L = 600 mm.

Ex. 3.2 Para a barra da fig. (3.6) descreva o problema de valor de contorno

48



resultante 10. 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çãoutilize 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. Usandoo 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, discretizandoas barras em cinco elementos iguais e usando FEM. Calcule a distribuiçãode 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/m 3 ao longo da barra, fluxo de calor q = 1.8×106 W/m 2 removendo calor em x = 0.2 m, ou seja, q (0.2) = 1.8 × 106 W/m 2, e condutividade térmica

10Cuidado: Para este caso a área não é constante! Uma dica é escrever A(x) =

A0 − A0−AL

Lx.

49



k = 6000 W/(mK). Assuma que a área da barra é 0.4

×10−3 m 2. 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ódulode 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 odeslocamento 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 mm 2 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 é sujeita 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éricos 12: altura de cada piso h = 15 ft, Carga nos andares P 1 = 50000 lb, P 2 = P 3 = 40000 lb, P 4 = 35000 lb, módulo de elasticidade E = 29 × 106 lb/in 2, área da secção transversal A = 21.8 in 2.

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 m 2

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

1002m 2; L =

5

1000m

T (L) = 84; −kAdT (0)

dx= 900 (3.118)

Utilize um elemento finito linear com dois nós apenas para facilitar 13.

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). Aolongo 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/m 2.C Determine

13Perceba que no caso geral para usar a função Matlab fornecida k(x) = kA = 3

10,

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− hP T + hP T ∞

= 0; 0 < x < L (3.119)P = 2(W + t)/; A = W t (3.120)

T (0) = 100; kAdT (L)

dx= −hA(T − T ∞) (3.121)

(3.122)

Os coeficientes na forma geral são:

k(x) = kA; p(x) = −hP ; q (x) = hP T ∞ (3.123)

NBC x = L; α = −hAkA

; β = hAT ∞kA

(3.124)

Use os seguintes dados numéricos para simulação W = 3 m, t = 0.3/100m; L = 20/100 m; k = 237; h = 30; T ∞ = 30; P = 6.006; A = 0.009;kA = 2.133; hP = 180.18; hP T ∞ = 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. (3.11). Utilize FEM empregando quatroelementos lineares com dois nós cada um. A temperatura do fluído varia

52



Fig. 3.10: Transferência de calor através de uma aleta [2].

linearmente de 80 oF na base até 200 oF no topo. A viscosidade do fluído em diferentes temperaturas é mostrada na tabela a seguir 14:

Tab. 3.2: Viscosidade µ do fluído.

Temperatura µ × 10−6

80 16120 14160 11200 7

A equação diferencial que governa este problema é descrita por:

d

dy

µ(y)

du

dy

− dP

dx= 0; 0 < y < h (3.125)

u(0) = 0; u(h) = 0 (3.126)

14Dica: Utilize algum método de interpolação para encontrar para encontrar uma fun-

ção do tipo µ(y).

53



Fig. 3.11: Perfil de velocidades de um fluído viscoso escoando entre duasplacas [2].

Use na simulação os seguintes dados: h = 0.3 ft e dP dx

= −2 × 10−5 lb/ft 3.

Ex. 3.11 Em coordenadas esféricas o potencial eletromagnético V é escritocomo uma função de θ pelo seguinte problema de valor de contorno:

d

dθ

sinθ

dV

dθ

= 0; α < θ < π/2 (3.127)

V (α) = V 0; V (π/2) = 0 (3.128)

sendo α um ângulo dado menor do que π e V 0 um potencial especificadopara este ângulo. Assumindo α = π/4 e V 0 = 1 e usando dois elementos finitos lineares, determine o potencial V em θ = 3π/8.

54



Capítulo 4

Elementos de Treliça, Viga e

Frame

A modelagem de sistemas mecânicos e estruturais complexos normal-mente envolvem elementos mais simples, como treliças, vigas e pórticos1.Exemplos são a modelagem de telhados, torres de antenas, pontes, máquinasde elevação, estruturas de automóveis, etc.

A análise de sistemas estruturais envolvendo estes elementos pode ser feitade maneira eficiente, usando o método dos elementos finitos. Neste sentido,a meta deste capítulo será fornecer todo o desenvolvimento para analisarum problema mecânico envolvendo treliças, vigas e frames. Em especial, oobjetivo final é encontrar qual a matriz de rigidez local ke para cada umdestes elementos, ilustrando exemplos de aplicação prática acompanhadosde procedimentos de programação numérica em ambiente Matlab R.

4.1 Treliças planas

Uma treliça é um sistema estrutural composto por elementos de barraarranjados de forma a suportarem apenas cargas axiais. A forma com que

estes elementos são montados é tal que o resultado final envolva sempre aformação de triângulos básicos entre as barras e as conexões pinadas. Estamontagem pode ser discretizada tal que as barras formem elementos e asconexões representem os nós em uma formulação FEM. Como as barras su-portam apenas carregamentos axiais, o elemento linear estudado no capítulo3 para o problema de deformação em barras axiais pode ser empregado nestaformulação.

1Em inglês frames .

55



A questão é que estes elementos são dispostos formando estruturas no

plano ou no espaço. Sendo assim, é necessário realizar uma transformaçãode coordenadas visando descrever todos os elementos em um mesmo sistemade referência, para assim ser possível realizar a montagem nas equações emum sistema global.

Neste contexto, um elemento de treliça plana é um elemento com defor-mação axial orientado arbitrariamente em um espaço bidimensional. Comomostrado na fig. (4.1), em um sistema s ao longo do eixo da barra o ele-mento é idêntico ao problema de deformação axial em barras. Em termosdo deslocamento assumido na coordenada s no elemento, tem-se a seguintefunção de interpolação:

u(s) =L − s

L

s

L

d1d2

; 0 ≤ s ≤ L (4.1)

Fig. 4.1: Coordenadas global e local para uma barra axial orientada no plano[2].

As equações do elemento neste sistema de coordenada local s são escritas

como:

EA

L

1 −1

−1 1

d1d2

=

P 1P 2

(4.2)

sendo que esta expressão pode ser escrita em uma forma matricial emfunção da matriz de rigidez do elemento no sistema de coordenadas locals kl, graus de liberdade local dl = (d1 d2)T e forças aplicadas nos nósrl = (P 1 P 2)T :

56



kldl = rl (4.3)

Porém uma vez que pode existir diferentes elementos de barra em umatreliça com orientações diferentes se faz necessário a definição de um sistemaglobal de coordenadas, como mostrado na fig. (4.1). Neste exemplo, define-se o sistema (x, y) e localiza-se todos os elementos em função deste. Assim,o que se faz é decompor os deslocamentos no sistema s na direção x, queaqui é denotado por u, e decompor na direção y, denotado por v. Portanto,no sistema global de coordenadas os graus de liberdade do elemento serãoquatro:

d =

u1

v1u2

v2

(4.4)

O mesmo é feito com as forças nodais:

r =

F x1F y1F x2F y2

(4.5)

A transformação de coordenadas é realizada em função do ângulo α entreo sistema local e global. Os deslocamentos u1 e v1 são escritos como:

u1 = d1cosα = d1ls (4.6)

v1 = d1cos(90 − α) = d1senα = d1ms (4.7)

sendo ls e ms o cosseno e seno do ângulo α, respectivamente. A relaçãode d1 com u1 e v1 é feita de maneira matricial:

u1v1 = lsms

d1 (4.8)

O comprimento do elemento de barra pode ser expressado em função dascoordenadas dos seus nós, que são: ponto 1 (x1, y1) e ponto 2 (x2, y2):

L =

(x2 − x1)2 + (y2 − y1)2 (4.9)

Por razões computacionais é interessante definir os valores de ls e ms emfunção da geometria da malha e não do angulo α. Assim:

57



ls = cosα =x2 − x1

L(4.10)

ms = senα =y2 − y1

L(4.11)

Multiplicando u1 por ls e v1 por ms e adicionando as duas equações obtémpor fim a transformação de global para local:

u1ls + v1ms = d1cos2α + d1sen2α = d1

cos2α + sen2α

≡ d1 (4.12)

Relações idêntica são obtidas tomando o nó 2. Com isto a transformaçãoentre os graus de liberdade no sistema global e local pode ser implementadapelas matrizes de transformação T. Para transformação global para local amatriz T é:

T =

ls ms 0 00 0 ls ms

(4.13)

Assim a transformação é feita por:

dl = Td (4.14)

Já a transformação inversa do sistema global para o sistema local podeser feita pela transposta de T, portanto:

d = TT dl (4.15)

Subsitituindo a eq. (4.14) na eq. (4.3) podemos obter a realização emfunção das coordenadas globais:

kldl = rl ⇒ klTd = rl (4.16)

Multiplicando ambos os lados da equação anterior por TT chega-se a:

TT klTd = TT rl (4.17)

Uma vez que r = TT rl e considerando que a matriz de rigidez do elementono sistema global de coordenadas (x, y) por ser escrita em função da matrizde transformação T e da matriz de rigidez do elemento no sistema local s,kl:

k = TT klT (4.18)

58



Com isto a eq. (4.17) pode ser escrita totalmente no sistema global de

coordenadas como:

kd = r (4.19)

Uma vez conhecidas as matrizes kl e T e os vetores dl e rl, a eq. (4.19)para o elemento torna-se:

EA

L

l2s lsms −l2s −lsms

lsms m2s −lsms −m2

s

−l2s −lsms l2s lsms

−lsms

−m2

s

lsms m2

s

u1

v1u2

v2

=

F x1F y1F x2F y

2

(4.20)

Por conveniência de implementação computacional é interessante adicio-nar o carregamento concentrado nos nós somente após a montagem globaldas equações. Assim, nas equações dos elementos o lado direito da equaçãoanterior fica zero:

EA

L



s


−lsms −m2s lsms m2

s

u1

v1u2

v2

=

0000

(4.21)

Com isto a matriz de rigidez do elemento de treliça plana no sistemaglobal de coordenadas (x, y) é descrita por:

k =EA

L



s


−lsms −m2s lsms m2

s

(4.22)

Estas equações são usadas para analisar qualquer montagem de barrasno plano. Depois de calculado os deslocamento nodais para elemento, asolução final é programada pela transformação dos deslocamentos nodais

para o sistema de coordenadas locais conforme a eq. (4.14). O deslocamentoaxial em qualquer ponto do elemento é calculado pela função de interpolaçãoda eq. (4.1).

O objetivo em uma análise estrutural é determinar a distribuição de ten-sões na treliça, assim é necessário calcular a deformação axial da barra. Noelemento de treliça plana a deformação axial é assumida constante em todoo elemento e calculada por2:

2Lembrando do capítulo anterior que u(s) = BT d, sendo B um vetor contendo as

derivadas das funções de forma.

59



= duds

= − 1L

1L d1

d2 = 1

L(−d1 + d2) (4.23)

A tensão axial no elemento é calculada com o conhecimento do módulode elasticidade E como σ = E. Já as forças axiais em cada barra podem sercomputadas por F = σA, sendo A a área da secção transversal em cada barra.A convenção de sinal usada nestas equações assume que tensão tracionandoé positiva e comprimindo é negativa.

4.1.1 Exemplo de implementação prática

Esta seção tem como objetivo apresentar uma aplicação prática de FEMpara análise estrutural de uma treliça. A fig. (4.2) mostra a treliça a ser usadaneste teste. Este exemplo é resolvido manualmente no livro [4]. Porém, nopresente texto irá ser apresentada uma solução implementada em ambienteMatlab R utilizando funções para o elemento treliça 2D descritas no livro [5].

Fig. 4.2: Treliça plana a ser analisada [4].

A estrutura é composta por seis barras conectadas em quatro nós. Todasas barras tem módulo de elasticidade de 21000 kgf/mm2. As áreas das seçõestransversais e os comprimentos das barras são variáveis em função de A = 300

60



mm2 e L = 1000 mm. As barras a e c têm área de 0.6A e comprimento 0.6L.

As barras b e d têm área de 0.8A e comprimento 0.8L. Já as barras e e f têm área A e comprimento L. A carga externa aplicada é P = 500 kgf.

A meta deste teste é calcular as reações de apoio, as forças internas emtodas as barras, assim como a distribuição de tensão em cada um dos mem-bros estruturais da treliça. O primeiro passo é construir a malha a ser usada,onde cada nó da treliça será equivalente a um nó da malha. O programaem Matlab R malha.m transcrito a seguir mostra como construir esta malharealizando a numeração dos nós e dos elementos. A fig. (4.3) apresenta oresultado gráfico desta rotina.

%--------------------------------------------------------%%% Descrição do Problema:% Construção de uma Malha 2D para uma estrutura do tipo treliça%%--------------------------------------------------------%clc; clear%--------------------------------------------------------%% Entrada de Dados%--------------------------------------------------------%

nel=6; % numero de elementosnnel=2; % numero de nós por elementondof=2; % numero de dofs por nónnode=4; % numero total de nós no sistemasdof=nnode*ndof; % numero total de dofs no sistema

%--------------------------------------------------------%% Coordenadas Nodais%--------------------------------------------------------%L=1000; % comprimento em mm

% gcoord(nó,x)= posição;% gcoord(nó,y)=posição % coordenadas dos nós

gcoord(1,1)=0.0; gcoord(1,2)=0.0; % Coordenada do nó Agcoord(2,1)=0.6*L; gcoord(2,2)=0.0; % Coordenada do nó Bgcoord(3,1)=0.6*L; gcoord(3,2)=0.8*L; % Coordenada do nó Cgcoord(4,1)=0.0; gcoord(4,2)=0.8*L; % Coordenada do nó D

61



%--------------------------------------------------------%

% Conectividade Nodal%--------------------------------------------------------%

% nodes(elemento,1)=nó

nodes(1,1)=1; nodes(1,2)=2; % Elemento (a) entre nós A e Bnodes(2,1)=2; nodes(2,2)=3; % Elemento (b) entre nós B e Cnodes(3,1)=3; nodes(3,2)=4; % Elemento (c) entre nós C e Dnodes(4,1)=4; nodes(4,2)=1; % Elemento (d) entre nós D e Anodes(5,1)=1; nodes(5,2)=3; % Elemento (e) entre nós A e C

nodes(6,1)=2; nodes(6,2)=4; % Elemento (f) entre nós B e D

%--------------------------------------------------------%% Plotar Malha%--------------------------------------------------------%% visualização da malha

figure(1)set(gcf,’name’,’Treliça - UNIOESTE’,’Numbertitle’,’off’)hold on

% Plotando os nós com círculos vermelhosset(plot(gcoord(:,1),gcoord(:,2),’ro’),’linewidth’,’linewidth’,3,3);

for i=1:6

aux1=nodes(i,1);aux2=nodes(i,2);x1=[gcoord(aux1,1),gcoord(aux2,1)];x2=[gcoord(aux1,2),gcoord(aux2,2)];set(plot(x1,x2),’linewidth’,’linewidth’,2,2);

axis off

end

% salvar figura com a malha no formato epssaveas(1,’malhatruss.eps’,’psc2’)

%--------------------------------------------------------%

62



Nó A Nó B

Nó C

elementoa

elemento

b

elemento

c

elemento

d

Nó D

elemento

e

elemento

f

Fig. 4.3: Visualização da malha implementada com o Matlab R.

O programa para análise irá usar três funções Matlab R feitas por [5]. Noinício de cada função existe uma descrição do comando, que pode ser obtidacom help mais o nome da função no prompt do Matlab R

. A primeira rotina

é chamada de feeldof e é usada para enumerar os elementos e os graus deliberdade em cada elemento.

function [index]=feeldof(nd,nnel,ndof)%----------------------------------------------------------% Objetivo:% Calcular os dofs associados com o elemento%% Comando:% [index]=feeldof(nd,nnel,ndof)%

% Descrição das variáveis:% index - vetor com os dofs do sistema associados com o iel-th elemento% iel - numero do elemento que os dofs do sistema irão ser calculados% nnel - numero de nós por elemento% ndof - numero de dofs por nó%---------------------------------

edof = nnel*ndof;k=0;for i=1:nnel

63



start = (nd(i)-1)*ndof;

for j=1:ndofk=k+1;index(k)=start+j;

endend

A segunda rotina fetruss2 contém as matrizes locais de rigidez e massade cada barra. Apesar da rotina retornar a matriz de massa para o elemento,

esta informação não será usada, uma vez que este teste tem como objetivofornecer apenas uma análise estática do problema.

function [k,m]=fetruss2(el,leng,area,rho,ls,ms,ipt)

%---------------------------------------------------% Objetivo:% Calcular as matrizes de massa e rigidez para o elemento de treliça 2D% nodal dof u_1 v_1 u_2 v_2%% Comando:% [k,m]=fetruss2(el,leng,area,rho,ls,ms,ipt)%% Descrição das Variaveis:% k - matriz de rigidez do elemento (4 x 4)% m - matriz de massa do elemento (4 x 4)% el - modulo de elasticidade% leng - tamanho do elemento% area - area da secao transversal% rho - densidade% ls = cos(beta); ms = sin(beta)

% ipt = 1 - matriz de massa distribuida% = 2 - matriz de massa concentrada%---------------------------------------------------

% matriz de rigidez

% c=cos(beta); s=sin(beta);c = ls; s = ms;k= (area*el/leng)*[ c*c c*s -c*c -c*s;...

64



c*s s*s -c*s -s*s;...

-c*c -c*s c*c c*s;...-c*s -s*s c*s s*s];

% consistent mass matrix

if ipt==1

m=(rho*area*leng/6)*[ 2*c*c+2*s*s 0 c*c+s*s 0;...0 2*c*c+2*s*s 0 c*c+s*s;...

c*c+s*s 0 2*c*c+2*s*s 0;...0 c*c+s*s 0 2*c*c+2*s*s];

% lumped mass matrix

else

m=(rho*area*leng/2)*[ c*c+s*s 0 0 0;...0 c*c+s*s 0 0;...0 0 c*c+s*s 0;...

0 0 0 c*c+s*s];

end

Por fim a terceira rotina feasmbl1 realiza a montagem global das matrizesde rigidez e massa a partir do conhecimento das matrizes para cada elemento.

function [kk]=feasmbl1(kk,k,index)%-------------------------------------------% Objetivo:

% Montar as matrizes globais do sistema%% Comando:% [kk]=feasmbl1(kk,k,index)%% Descrição das variáveis:% kk - matriz do sistema todo% k - matriz do elemento% index - vetor com os dofs associados com o elemento

65



%-------------------------------------------

edof = length(index);for i=1:edof

ii=index(i);for j=1:edof

jj=index(j);kk(ii,jj)=kk(ii,jj)+k(i,j);

endend

É importante observar que estas rotinas são gerais e podem ser usadaspara qualquer problema envolvendo estruturas do tipo treliça 2D, sendo queas funções feeldof e feasmbl1 podem inclusive ser adaptadas para outros tiposelemento, como será visto no decorrer deste capítulo.

Para efetuar os cálculos das tensões e forças atuantes o programa modelofoi implementado. A seguir é feita a transcrição desta rotina:

%------------------------------------------------%% Descrição do Problema% Analise Estática de uma Estrutura Tipo Treliça (2D)%% Descrição das variáveis% k = matriz de rigidez do elemento% m = matriz de massa do elemento% kk = matriz de rigidez do sistema% mm = matriz massa do sistema% index = vetor contendo dofs do sistema associados com cada elemento% gcoord = matriz de coordenadas global% prop = matriz de propriedades do elemento% propA = áreas das barras

% nodes = matriz de conectividade para cada elemento% bcdof = vetor contendo dofs associados com as condições de restrição% bcval = vetor contendo os valores das condições de restrições associadas% com os dofs no vetor ’bcdof’% x = vetor de deslocamentos nodais% F = forças externas atuantes nos nós% FR = reações de apoio% d = deslocamentos axiais em cada elemento% Eps = deformação axial nas barras

66



% Sigma = tensão nas barras

%------------------------------------------------

clc; clear; close all

%------------------------------------------------% Abertura do Arquivo de dados com a malha%-----------------------------------------------

malha

%------------------------------------------% Propriedades Geometricas e dos Materiais%------------------------------------------

prop(1)=21000; % modulo de elasticidade [kgf]

A=300; % area da seção transversal [mm^2] - Barras tem áeas diferentespropA(1) = 0.6*A; % barra (a) tem ·rea 0.6*A;propA(2) = 0.8*A; % barra (b) tem ·rea 0.8*A;propA(3) = 0.6*A; % barra (c) tem ·rea 0.6*A;

propA(4) = 0.8*A; % barra (d) tem ·rea 0.8*A;propA(5) = A; % barra (e) tem ·rea A;propA(6) = A; % barra (f) tem ·rea A;

prop(2)=1; % densidade (seria usado este valor para cáculo da matriz de massa)

%----------------------------% Inicialização com Zero%----------------------------

kk=zeros(sdof,sdof); % matriz de ridez do sistema

mm=zeros(sdof,sdof); % matriz de massa do sistemaindex=zeros(nnel*ndof,1); % vetor de indices

%--------------------------% Loop para os Elementos%--------------------------

for iel=1:nel % loop para o numero total de elementos

67



nd(1)=nodes(iel,1); % 1st no conectado para o (iel)-th elemento

nd(2)=nodes(iel,2); % 2nd no conectado para o (iel)-th elemento

x1=gcoord(nd(1),1); y1=gcoord(nd(1),2); % coordenada do 1.o nóx2=gcoord(nd(2),1); y2=gcoord(nd(2),2); % coordenada do 2.o nó

leng=sqrt((x2-x1)^2+(y2-y1)^2); % tamanho do elemento

% Duas formas possíveis para calcular os angulos de orientação das barras% if (x2-x1)==0;% beta=2*atan(1);% angulo entre o sistema local e global

% else% beta=atan((y2-y1)/(x2-x1));% end

ls = (x2-x1)/leng;ms = (y2-y1)/leng;

el=prop(1); % extraindo modulo de elasticidadearea=propA(iel); % extraindo area da seção transversal% (diferente para cada elemento)

rho=prop(2); % extraindo densidade (para massa)

index=feeldof(nd,nnel,ndof); % extraindo os dofs do sistema para o elemento

ipt=1; % escolha da opção de matriz de massa continua

[k,m]=fetruss2(el,leng,area,rho,ls,ms,ipt); % matrizes dos elementos

kk=feasmbl1(kk,k,index); % montagem da matriz de rigidez do sistema% mm=feasmbl1(mm,m,index); % montagem da matriz de massa do sistema% Pois não estamos interessados na matriz de massa

% Isto seria interessante para análise dinâmica - cálculo de forças de inércia

end

%-------------------------------------------% Aplicando as Restrições%-------------------------------------------% Os quatro primeiros graus de liberdade são restritos

68



K = kk(5:length(kk),5:length(kk));

% Vetor com forças externas atuantes no nó C (dof 5, 6)F = [1000; 500; 0;0];

%-------------------------------------------% Solução do sistema linear: F=Kx%-------------------------------------------

x=K\F; % x: vetor de deslocamentos nodais nos nós C e D

x = [zeros(4,1); x]; % vetor de deslocamentos nodais na estrutura toda% notando que pelas condições de contorno x(1),x(2),x(3) e x(4) são nulos

%-------------------------------------------% As forças atuantes nos nós (A) e (B) são as reações de apoio% que agem nos dofs 1,2,3 e 4% Observe que no nó (D) não atuam nenhuma força (dof 7,8)%-------------------------------------------

FR = kk*x; % calculo das reações de apoio

%-------------------------------------------% Calculo dos deslocamentos axiais -% Deformações e tensões nas barras

for iel=1:nel % loop para o numero total de elementos

nd(1)=nodes(iel,1); % 1st no conectado para o (iel)-th elementond(2)=nodes(iel,2); % 2nd no conectado para o (iel)-th elemento

x1=gcoord(nd(1),1); y1=gcoord(nd(1),2); % coordenada do 1.o nó

x2=gcoord(nd(2),1); y2=gcoord(nd(2),2); % coordenada do 2.o nó

leng=sqrt((x2-x1)^2+(y2-y1)^2); % tamanho do elemento

index=feeldof(nd,nnel,ndof); % extraindo os dofs do sistema para o elemento

ls = (x2-x1)/leng;ms = (y2-y1)/leng;

69



% Transformação de coordenadas globais para locais

d(:,iel) = [ls ms 0 0;0 0 ls ms]*x(index);

teste(:,iel)=index;Eps(iel) = 1/leng*(-d(1,iel)+d(2,iel)); % Cálculo da deformação

end

Sigma = el*Eps;Fax(1) = Sigma(1)*propA(1);Fax(2) = Sigma(2)*propA(2);

Fax(3) = Sigma(3)*propA(3);Fax(4) = Sigma(4)*propA(4);Fax(5) = Sigma(5)*propA(5);Fax(6) = Sigma(6)*propA(6);%--------------------------------

As tabelas (4.1), (4.2) e (4.3) apresentam os resultado final alcançadoscom com a análise estrutural via FEM utilizando as rotinas anteriores.

Tab. 4.1: Deslocamentos nodais (mm).

Nó Coordenada x Coordenada yA 0.0 0.0B 0.0 0.0C 0.3392 −0.0509D 0.2781 0.0814

Tab. 4.2: Forças externas aplicadas nos nós (kgf).

Nó Força em x Força em yA −615.4 −1333.36B −384.6 833.3C 1000.0 500.0D 0.0 0.0

70



Tab. 4.3: Deformações , tensões σ e forças internas axiais F nas barras.Barra (mm/mm) σ (kgf/mm2) F (kgf) Estado

a 0.0 0.0 0.0 Não é solicitadab −6.36 × 10−5 −1.335 −320.5 Compressãoc 1.017 × 10−4 2.137 384.6 Traçãod 1.017 × 10−4 2.137 512.8 Traçãoe 1.628 × 10−4 3.42 1025.6 Traçãof −1.017 × 10−4 −2.137 −641.03 Compressão

Estas rotinas são muito próximas ao que existe no interior de um pro-grama comercial de FEM, como Ansys R ou Nastran R. Este exemplo, ape-sar de simples, mostra quais os passos básicos necessários para implementarcomputacionalmente qualquer aplicação envolvendo FEM. Os estudantes sãoencorajados a resolver o exemplo acima usando algum software comercialcomo Ansys R. Um bom tutorial para iniciantes no ANSYS R é encontradono site http://www.mece.ualberta.ca/Tutorials/ansys/.

4.2 Deformação transversal em vigas

Muitos sistemas reais podem ser modelados como sendo vigas, por exem-plo, colunas, asas de aviões, eixos de rotores3, etc. Uma viga é um elementoestrutural que pode suportar tanto carga axial, como em uma barra, quantocargas transversais (cortante) e torção. Assim, uma viga pode apresentar aolongo de seu eixo longitudinal uma distribuição de tensões normais σ(x), quepode ser causada tanto por eventuais forças axiais, quanto por momentosfletores4, e uma distribuicão de tensões de cisalhamento τ (x), provocada oupor torção e/ou por forças cortantes. Do ponto de vista geométrico, uma vigaé definida como um sistema onde o comprimento é bem mais representativo

do que a largura e espessura5

.Em cursos básicos de mecânica geral e mecânica dos sólidos o objetivoera determinar a distribuição de momentos fletores M (x) e cortantes V (x) aolongo de vigas e depois calcular quais as tensões máximas atuantes no sistema

3Ressalta-se que "90%"dos equipamentos industriais envolvem componentes rotativos,como por exemplo, compressores, turbinas, ventiladores, geradores, motores, etc., portantosaber analisar um sistema destes é importante para um engenheiro projetista.

4Provocados pela possível existência de cargas transversais5Um elemento linear pode ser usado em aplicação FEM.

71



envolvendo expressões relacionando a geometria da viga, propriedades físicas

do material e os esforços de flexão e cortante estimados. Destaca-se que oproblema estástico e/ou dinâmico do movimento e deformação de uma vigapode ser solucionado de forma exata com o uso de uma equação diferencialparcial envolvendo um problema de valor de contorno e condição inicial. Nocaso de envolver somente tensões normais este modelo é conhecido como Vigade Euler-Bernoulli, já no caso de envolver também tensões de cisalhamentoo modelo é conhecido como Viga de Timoshenko. Porém, como constatadono decorrer desta seção esta formulação é muito trabalhosa e pouco prática, justificando o uso de um procedimento numérico.

Estas metodologias podem ser adaptadas para análise e projeto de vigas e

eixos. Entre os métodos clássicos para calcular M (x) e V (x), além da soluçãoexata da equação diferencial parcial que descreve o problema, se destacamo método das seções e o método das funções singulares. Nesta seção seráapresentado um procedimento para análise estática de viga usando FEM. Agrande vantagem de FEM é ser capaz de trabalhar com vigas envolvendo umnúmero maior de restrições, apoios e descontinuidade, nos casos em que osprocedimentos convencionais ficam muito exaustivos. A meta do capítulo éapresentar a matriz de rigidez k para o elemento de viga de Euler-Bernoullie exemplificar o uso de FEM em problemas de análise estrutural em rotinascomputacionais programadas no ambiente Matlab R.

4.2.1 Equação diferencial para deflexão em vigas

A área da secção transversal de uma viga pode ser assimétrica, porém,a formulação apresentada aqui se restringirá aos casos onde a secção trans-versal de uma viga é simétrica. A fig. (4.4) mostra uma viga sujeita a umcarregamento q (x, t) transversal ao longo do eixo x. O momento de inérciade área6 é dado por I (x). Considera-se que esta viga tenha uma distribuiçãode massa por unidade de comprimento dada por m(x).

A meta aqui é obter uma equação diferencial que seja governante domovimento transversal v(x, t) desta viga causado por eventuais condições de

valor inicial e valor de contorno. Para isto será realizado um diagrama decorpo livre de um elemento diferencial posicionado a uma distância x e comcomprimento dx nesta viga, como mostrado na fig. (4.5). A cortante V e omomento M na posição x é mostrado nesta figura. Usando uma expansãocom a série de Taylor a cortante na posição x + dx será V + (∂V/∂x) dx e omomento será M + (∂M/∂x) dx.

Como a viga tem massa m(x) e se movimenta, pela segunda lei de New-

6Também conhecido como segundo momento de área

72



Fig. 4.4: Flexão em vigas [2].

Fig. 4.5: Forças e momentos agindo em um elemento diferencial da viga [2].

ton existe uma força de inércia agindo em sentido contrário a tendência domovimento. Esta força é dada em função da aceleração v(x, t) = ∂ 2v/∂t2.Assim a força de inércia é m(x)dxv(x, t).

Para garantir o equilíbrio dinâmico, o somatório de forças no elementodiferencial dx da viga deve atender o princípio de D’Alembert:

q (x, t)dx + V −V +

∂V

∂x dx− m(x)dxv(x, t) = 0 (4.24)o que integrando em x leva a expressão:

mv(x, t) +∂V

∂x= q (x, t) (4.25)

Já para garantir o equilíbrio da soma de momentos do lado direito doelemento diferencial dx, têm-se que7:

7Adotando como convenção que o momento no sentido anti-horário é positivo

73



−q (x, t)dxdx2

− M − V dx + M + ∂M ∂x

dx = 0 (4.26)

desprezando os termos de segunda ordem, chega-se a:

V =∂M

∂x(4.27)

Neste ponto é importante observar que no elemento diferencial dx a fibrasda face superior da viga são comprimidas e as fibras inferiores tracionadas,sendo assim a aplicação de uma carga transversal em uma viga provoca tensãonormal que varia ao longo da secção transversal. Observa-se também que em

uma posição ao longo da secção transversal y da viga, as fibras permanecemsem nenhuma deformação. Esta posição é conhecida como linha neutra edefine o raio de curvatura de uma viga. Este conceito é a base para afirmarque a flexão causa o aparecimento de uma deformação axial u(x) na vigaque é relacionada a posição transversal y e ao raio de curvatura desta8. Estadeformação axial u(x) é escrita como:

u(x) = −ydv

dx(4.28)

A eq. (4.28) é muito útil, pois permitirá obter uma relação entre tensãonormal é momento fletor. A deformação axial x pode ser calculada em

função do resultado da eq. (4.28)9:

x =du

dx= −y

d2v

dx2(4.29)

Para um material elástico e linear a relação tensão/deformação é definidapela Lei de Hooke:

σx = Ex = −Eyd2v

dx2(4.30)

A eq. (4.30) mostra claramente que a tensão normal varia linearmente

com relação a distância y da secção transversal em cada ponto x ao longo docomprimento de uma viga.

Considerando que um força qualquer em uma posição x é dF e atua naseção transversal da viga na mesma direção da carga q (x, t), sendo dada por:

8Para relembrar este conceito pesquise em um livro de resistência dos materiais umcapítulo sobre flexão em vigas.

9O sinal negativo na expressão mostrada se refere ao fato de adotar que na face superioras fibras da viga são comprimidas e na face inferior tracionadas

74



dF = σxdA (4.31)

sendo dA um elemento diferencial da área da secção transversal. O mo-mento desta força na secção transversal é dado por dM :

dM = −ydF = −yσxdA (4.32)

Integrando em dA a equação anterior e considerando a eq. (4.30) obtém-se:

M (x) =

− A σxydA = A Eyd2v

dx2

dA = E d2v

dx2 A y2dA (4.33)

Lembrando que o momento de inércia ou momento de segunda ordem deárea é I :

I =

A

y2dA (4.34)

Assim:

M (x) = EI d2v

dx2(4.35)

Combinando as eqs. (4.35) e (4.30):

σx = −My

I (4.36)

A eq. (4.36) pode ser usada para projeto de vigas e eixos considerandocalcular a geometria do sistema10 em função das tensões a serem solicitadas.

Note que a cortante V na eq. (4.27) pode ser escrita da forma:

V =∂M

∂x=

∂

∂x

EI

∂ 2v

∂x2

(4.37)

Substituindo a expressão anterior na equação de equilíbrio dinâmicoobtém-se a equação de viga de Euler-Bernoulli:

mv(x, t) +∂ 2

∂x2

EI

∂ 2v

∂x2

= q (x, t) (4.38)

A eq. (4.38) é uma equação diferencial parcial representado a deflexãov(x, t) de um sistema de parâmetros concentrados. Esta equação é de quartaordem e envolve uma derivada de segunda ordem no tempo t e uma derivada

10Informação está contida dentro do momento de inércia de área.

75



de quarta ordem em x. Para solucioná-la de forma exata é necessário o

conhecimento de duas condições iniciais em t e quatro condições de contornoem x. A solução da equação de viga de Euler Bernoulli fornece os modos devibrar da estrutura φ(x) em todo o domínio envolvendo infinitas frequenciasnaturais.

Aqui o interesse é analisar a condição estática e irá se admitir que asforças de inércia são nulas, assim o problema fica apenas do tipo PVC:

d2

dx2

EI

d2v

dx2

= q (x) (4.39)

4.2.2 Condições de contorno em vigasAs condições de contorno em viga envolvem especificações do desloca-

mento v e suas primeiras tês derivadas em x. A fig. (4.6) mostra algumascondições de contorno comuns em vigas.

Fig. 4.6: Condições de contorno comuns em vigas [2].

Relembrando que a derivada segunda em x do deslocamento v é relacio-

nada diretamente ao momento fletor e a cortante, conforme visto na seçãoanterior, em um ponto x0 aparecerá as condições de contorno vistas na tab.(4.4).

Assim as condições de contorno em uma viga engastada-livre são:

M = EI ∂ 2v(x0)

∂x2= 0 (4.40)

V =∂

∂x

EI

∂ 2v(x0)

∂x2

= 0 (4.41)

76



Tab. 4.4: Condições de contorno11

em viga em uma posição arbitrária x0.Variável Lado esquerdo da seção Lado direito da seção

Deslocamento v(x0) = vx0 v(x0) = vx0Rotação θ(x0) = ∂v(x0)

∂x= θx0 θ(x0) = ∂v(x0)

∂x= θx0

Momento −EI ∂ 2v(x0)∂x2

= M x0 EI ∂ 2v(x0)∂x2

= M x0Cortante EI ∂

3v(x0)∂x3

= F x0 −EI ∂ 3v(x0)∂x3

= F x0

As condições em posições fixa:

v = 0 (4.42)

θ =∂v

∂x= 0 (4.43)

As condições de contorno em um suporte em suporte (pino) simples são:

v = 0 (4.44)

M = EI

∂ 2v(x0)

∂x2 = 0 (4.45)

Já em um suporte deslizante são:

θ =∂v

∂x= 0 (4.46)

V =∂

∂x

EI

∂ 2v(x0)

∂x2

= 0 (4.47)

Condições de contorno mais complicadas podem ser impostas, principal-

mente quando uma viga é suportada por outros membros elásticos, comobarras ou molas. O uso de molas para suportar vigas é muito comum namodelagem de máquinas rotativas. No caso a mola com rigidez k representao mancal no modelo matemático. A força F e em uma mola posicionada emx é:

F e(x0) = kv(x0) (4.48)

Esta força deve ser balanceada pela força cortante neste suporte ou pelomomento fletor, caso a mola seja torsional.

77



Para o caso da viga ser apoiada em uma mola k a condição de contorno

neste ponto fica sendo:

∂

∂x

EI

∂ 2v(x0)

∂x2

− kv = 0 (4.49)

No caso de uma viga apoiada em uma mola torsional tem-se:

EI ∂ 2v(x0)

∂x2+ kv = 0 (4.50)

4.2.3 Exemplo de solução exata em vigas

É possível obter uma solução exata para obter o deslocamento transversalv de uma viga, porém o procedimento não é prático. Apenas para ilustraresta tarefa, considere a viga da fig. (4.7). A viga é sujeita a uma carga linearcrescente do tipo q (x) = ax/L, sendo a uma constante. Este problema devalor de contorno é escrito como:

EI d4v

dx4=

ax

L; 0 < x < L (4.51)

v(0) = 0 EI d2v(0)

dx2

= 0 (4.52)

EI d3v(L)

dx3− kv(L) = 0 EI

d2v(L)

dx2= 0 (4.53)

Fig. 4.7: Viga pinada e apoiada em uma mola.

Integrando quatro vezes a equação diferencial que descreve o problema:

78



EI d3v

dx3− ax2

2L= c1 (4.54)

EI d2v

dx2− ax3

6L= c1x + c2 (4.55)

EI dv

dx− ax4

24L= c1

x2

2+ c2x + c3 (4.56)

EI v − ax5

120L= c1

x3

6+ c2

x2

2+ c3x + c4 (4.57)

sendo c1, c2, · · · , c4 são constantes de integração. Usando as condições de

contorno:

v(0) = 0 ⇒ c4 = 0 (4.58)

EI d2v(0)

dx2= 0 ⇒ c2 = 0 (4.59)

EI d2v(L)

dx2= 0 ⇒ −20aL3 − 120L2c1

120L= 0 ⇒ c1 = −aL

6(4.60)

EI d3v(L)

dx3− kv(L) = 0

⇒ −−60aL2

−120Lc1

120L +

k (

−aL5

−20L4c1

−120L2c3)

120EI L = 0

⇒ c3 = −−120aEI − 7akL3

360k(4.61)

Assim a solução exata do problema é:

v(x) =ax (120EI L + k (7L4 − 10L2x2 + 3x4))

360EIkL(4.62)

Apenas para ilustrar a fig. (4.8) mostra o deslocamento v(x) assumindoEI = 1, L = 1, a = 1 e com vários k para as molas.

Com o conhecimento da função v(x) é possível determinar o momento

fletor M (x) em toda a extensão de x, a cortante V (x) e a distribuição detensões de forma exata. Porém, para aplicações práticas obter a soluçãoexata do problema de valor de contorno associado a viga é inviável.

4.2.4 Diagramas de cortante e momento usando o mé-todo das seções

Apenas para relembrar um método usado para estimar as funções demomento fletor M (x) e cortante V (x) esta seção mostra a partir de um

79



0 0.2 0.4 0.6 0.8 10

1

2

3

4

5

6

7

8

9

x 10−3

Comprimento

v ( x )

k=100

k=200

k=500

k=1000k=10000

Fig. 4.8: Soluções exatas para a viga com várias constantes.

exemplo resolvido o uso do método das seções. Com este método levanta-seV (x) e M (x) e consequentemente σx(x) sem a necessidade da solução exatade v(x).

Considere a viga da fig. (4.9). A primeira etapa é a construção de umdiagrama de corpo livre para cálculo das reações de apoio nos vínculos12. Asreações de apoio neste exemplo são:

M A = 0 ∴ BL − woL

4

L

3+

L

2

⇒ B =

5w0L

24(4.63)

F y = 0 ∴ A + B = woL4 ⇒ A = woL24 (4.64)

Como a viga tem uma descontínuidade em x = L/2 irá se encontrar duasfunções V (x) e M (x) para dois trechos distintos. O primeiro trecho é entre0 < x < L/2, neste caso a cortante e o momento são:

12Aconselha-se que os estudantes relembrem o cálculo de carregamentos concentradosa partir de carregamentos distribuídos em vigas através da consulta em livros básicos deestática.

80



Fig. 4.9: Exemplo de viga com carregamento.

F y = 0 ∴ V (x) =

woL

24(4.65)

M = 0 ∴ M (x) − woL

24x ⇒ M (x) =

woL

24x (4.66)

Já o segundo trecho da viga é L/2 < x < L. A cortante V (x) neste trechoé:

F y = 0 ∴

woL

24− wo

L

x − L

2

2

− V (x) = 0

V (x) =woL

24− wo

L

x − L

2

2

(4.67)

Já o momento fletor M (x) neste trecho é:

M = 0 ∴ M (x) − woL

24x +

wo

L

x − L

2

2 1

3

x − L

2

= 0

M (x) =woL

24x − wo

3L

x − L

2

3

(4.68)

Uma vez tendo estas funções pode-se calcular os valores máximos de V (x)e M (x) e conhecendo a geometria e propriedades da viga calcular as tensões

81



máximas atuantes. Este procedimento também pode ser dificil de ser usado

na prática, principalmente se a viga contiver muitas descontinuidades, forestaticamente indeterminada para cálculo das reações de apoio, etc. Assim,é plenamente justificável o uso de uma formulação de vigas via FEM.

4.2.5 Elemento de viga com dois nós

As condições essenciais para o problema de viga de Euler-Bernoulli en-volvem o deslocamento tranversal v e a rotação θ = dv/dx. Portanto paramontar a solução via FEM serão necessários dois graus de liberdade, v e θpor nó em um elemento linear. Com isto um elemento de viga tem ao todo

quatro graus de liberdade. Visando a simplicidade na integração o momentode inércia de área e a carga distribuída são assumidas constantes no ele-mento. Além disto, as cargas concentradas F 1 e F 2 e os momentos M 1 e M 2são permitidos serem aplicados apenas nos nós.

O elemento de viga a ser usado é mostrado na fig. (4.10). Mesmo sendobem simples, este elemento pode ser usado para a grande maioria dos exem-plos práticos envolvendo vigas. Para vigas uniformes apenas um elemento játem capacidade de fornecer a solução exata. Para vigas com muitas desconti-nuidades, para cada variação seja de carga, presença de vínculo ou mudançade propriedade física e geométrica.

Fig. 4.10: Elemento de viga com flexão.

82



4.2.6 Solução cúbica assumida

Uma vez que os graus de liberdade em cada nó são o deslocamento trans-versal v e a rotação θ, uma função de interpolação apropriada para o elementodeverá ter quatro constantes (ideal para solução de problemas de quarta or-dem, como é o caso do problema de Viga de Euler-Bernoulli). Conformeestudado no capítulo anterior a função de interpolação Hermitiana é apropri-ada para este caso. No capítulo anterior foi mostrado que é possível escreverfunções de forma N i(x) em função das coordenadas dos nós x1 e x2. Porém,partir por este caminho para obter a solução do elemento se mostra bastantetrabalhosa e monótono. Uma expressão mais conveniente para as funções deforma podem ser obtidas introduzindo mudança de coordenadas:

s = x − x1; x1 ≤ x ≤ x2; 0 ≤ s ≤ L (4.69)

ds = dx ⇒ dv

dx=

dv

ds(4.70)

Aplicando esta transformação de coordenadas nas funções de forma dainterpolação Hermitiana obtém-se:

N 1(s) =

2s3

L3 −3s2

L2 + 1 (4.71)

N 2(s) =s3

L2− 2s2

L+ s (4.72)

N 3(s) =3s2

L2− 2s3

L3(4.73)

N 4(s) =s3

L2− s2

L(4.74)

O vetor de funções de forma no elemento é então dado por:

N

T

= N 1(s) N 2(s) N 3(s) N 4(s) (4.75)assim a solução assumida v(s) é:

v(s) = NT

v1θ1v2θ2

= NT d (4.76)

a solução assumida para a rotação θ(s) no elemento é:

83



θ(s) = dv(s)ds

= BT d (4.77)

sendo BT = B1(s) B2(s) B3(s) B4(s) onde Bi(s) = dN i/ds.O objetivo no método dos elementos finitos é obter uma solução aproxi-

mada na forma algébrica da eq. (4.39):

kd = rq + rp (4.78)

Assim, em vez de ter que resolver a equação diferencial de Viga de Euler-Bernoulli, pode-se obter a solução do problema para cada elemento caso seconheça a matriz de rigidez k deste elemento.

4.2.7 Matriz de rigidez local em vigas

O método de Galerkin envolvendo a forma fraca e o procedimento de in-tegração por partes pode ser usado aqui para obter r p, rq e k. Porém, para oproblema de viga é mais fácil obter a solução usando o método de Rayleigh-Ritz13, que utiliza conceitos variacionais relacionados a energia. Felizmente,independente do método a ser empregado a matriz de rigidez do elementosempre será da forma já apresentada no capítulo 3, considerando que o ele-mento14 se estende de s = 0 até s = L:

k = L0

BEI BT = EI L0

BBT (4.79)

Assim para encontrar uma expressão para a matriz de rigidez local doelemento de viga é necessário primeiramente derivar em s as funções de formada interpolação Hermitiana. Assim:

BBT =

36(L2s)2

L612(2L−3s)(L−2s)

L5−36(L−2s)2

L612(L−3s)(L−2s)

L512(L−3s)(L−2s)

L54(2L−3s)2

L4−12(2L−3s)(L−2s)

L54(L−3s)(2L−3s)

L4

−36(L−2s)2

L6 −12(2L−3s)(L−2s)

L536(L−2s)2

L6 −12(L−3s)(L−2s)

L512(L−3s)(L−2s)

L54(L−3s)(2L−3s)

L4−12(L−3s)(L−2s)

L54(L−3s)2

L4

(4.80)Subsituindo a eq. (4.80) na eq. (4.79) obtém-se:

13Que deverá ser estudado nos próximos capítulos.14Além de se assumir que EI é constante no elemento.

84



k = EI

L0

36(L2s)2

L6

ds

L0

12(2L−3s)(L−2s)L5

ds · · · L

0

12(L−3s)(L−2s)

L5

ds

L0

4(2L−3s)2

L4

ds · · ·

......

...

(4.81)

Após realizar a integração da eq. (4.81) obtém-se finalmente a matriz derigidez do elemento de viga de Euler-Bernoulli:

k =

EI

L3 12 6L −12 6L6L 4L2

−6L 2L2

−12 −6L 12 −6L6L 2L2 −6L 4L2

(4.82)

Também conforme foi visto no capítulo 3 o vetor de carga15 rq é descritoutilizando o vetor de funções de forma N:

rq =

L0

Nqds = q

L0

2s3

L3− 3s2

L2+ 1)ds L

0

s3

L2− 2s2

L+ s)ds

L

0 3s2

L2− 2s3

L3)ds L

0 s3L2

− s2L

)ds

=qL

2

1L6

1

−L

6

(4.83)

Com a eq. (4.83) observa-se que o carregamento distribuído ao longode um elemento de viga é considerado dividido e concentrados em cada nó,conforme visto na fig. (4.11).

Fig. 4.11: Carregamento distribuído concentrado nos nós.

Já o vetor r p representa o carregamento concentrado no nós do elementode viga. Para o caso geral ele é dado por:

15Assumindo ser constante no elemento.

85



r p =

F 1M 1F 2M 2

(4.84)

Assim a equação algébrica do elemento de viga de Euler-Bernoulli, eq.(4.78) é escrita em função dos resultados das eqs. (4.82), (4.83) e (4.84):

EI

L3

12 6L −12 6L6L 4L2 −6L 2L2

−12 −6L 12 −6L6L 2L2 −6L 4L2

v1θ1

v2θ2

=

qL

2

1L6

1−L6

+

F 1M 1

F 2M 2

(4.85)Assim como no caso da análise em treliças, normalmente as cargas ex-

ternas nos nós, vetor r p, são adicionadas somente no final da montagem dosistema global. Assim a eq. (4.85) tornasse:

EI

L3

12 6L −12 6L6L 4L2 −6L 2L2

−12

−6L 12

−6L

6L 2L2 −6L 4L2

v1θ1v2

θ2

=

qL

2

1L6

1

−L6

(4.86)

Para um problema qualquer envolvendo uma viga, a montagem da matrizglobal de rigidez K pode ser feita de forma convencional. A partir da soluçãodo sistema linear resultante pode se encontrar o vetor d = (v1 θ1 v2 θ2)T

de deslocamentos nos nós da viga. A solução final v(s) é encontrada pelafunção de interpolação Hermitiana para cada elemento, repetida abaixo:

v(s) = NT d (4.87)

Já a rotação para cada elemento é obtida por:

θ(s) = BT d (4.88)

O momento fletor M (s) é calculado por:

M (s) = EI d2v(s)

ds2= EI HT d (4.89)

sendo HT = (dB1(s)/ds dB2(s)/ds dB3(s)/ds dB4(s)/ds ) .Já a cortante V (s) é:

86



V (s) = EI d3v(s)ds3

= EI DT d (4.90)

sendo DT = (dH 1(s)/ds dH 2(s)/ds dH 3(s)/ds dH 4(s)/ds ) . Se fordo interesse obter a tensão normal em x da viga pode-se usar o valor encon-trado para o momento fletor M (s) e aplicar a eq. (4.36). Para uma viga comseção simétrica e com altura h a tensão normal máxima devido ao momentofletor será dada por:

σx,max =M (h/2)

I (4.91)

Apesar da teoria de vigas de Euler-Bernoulli não prever na sua formulaçãodeformações devido a forças cortantes, as tensões de cisalhamento τ provoca-das por elas podem ser calculadas utilizando conceitos básicos de resistênciados materiais que levam a equação:

τ =V Q

It(4.92)

sendo V a cortante estimada da eq. (4.90), I o momento de inércia deárea e Q e t são relacionados a pontos na seção transversal onde a tensãode cisalhamento está sendo calculada. t é a espessura da seção transversal eQ o momento de área cisalhante sobre a linha neutra. A tensão máxima decisalhamento para um viga com seção simétrica pode ser calculada por:

τ max =V Qmax

It(4.93)

Nas seções a seguir são apresentados alguns exemplos de solução numéricade deflexão em vigas, cálculos de momento, cortante e tensões usando FEM.Os exemplos são simulados com o Matlab R.

4.2.8 Exemplo de aplicação com carregamento concen-

tradoA fig. (4.12) mostra uma viga com carregamentos concentrados aplicados

em alguns pontos. O objetivo é calcular o deslocamento transversal v(x)da viga, as reações de apoio, o momento fletor M (x), a cortante V (x) ea distribuição de tensão normal σ(x). A viga é assumida ser feita de açocomo módulo de elasticidade E = 210 GPa. A seção transversal da viga éretangular com b = 0.144 m e h = 0.399 m.

A solução deste problema será realizada com o auxílio de rotinas imple-mentadas em Matlab R. As funções feeldof e feasmbl1 são as mesmas usadas

87



Fig. 4.12: Viga com carregamento concentrado aplicado.

no exemplo da treliça 2D mostrado anteriormente. A matriz de rigidez parao elemento de viga é calculada a partir da função febeam1 transcrita a seguir.

function [k,m]=febeam1(el,xi,leng,area,rho,ipt)

%-------------------------------------------------------------------% Objetivo:% Stiffness and mass matrices for Hermitian beam element% nodal dof v_1 theta_1 v_2 theta_2%% Synopsis:% [k,m]=febeam1(el,xi,leng,area,rho,ipt)%% Variable Description:% k - matriz de rigidez do elemento(size of 4x4)% m - matriz massa do elemento (size of 4x4)

% el -modulo de elasticidade% xi - second moment of inertia of cross-section% leng - tamanho do elemento% area - area da seçãao transversal da viga% rho - densidade% ipt = 1: consistent mass matrix% 2: lumped mass matrix% otherwise: diagonal mass matrix%-------------------------------------------------------------------

88



% stiffness matrix

c=el*xi/(leng^3);k=c*[12 6*leng -12 6*leng;...

6*leng 4*leng^2 -6*leng 2*leng^2;...-12 -6*leng 12 -6*leng;...

6*leng 2*leng^2 -6*leng 4*leng^2];


if ipt==1

mm=rho*area*leng/420;m=mm*[156 22*leng 54 -13*leng;...

22*leng 4*leng^2 13*leng -3*leng^2;...54 13*leng 156 -22*leng;...

-13*leng -3*leng^2 -22*leng 4*leng^2];


elseif ipt==2

m=zeros(4,4);mass=rho*area*leng;m=diag([mass/2 0 mass/2 0]);

% diagonal mass matrix

else

m=zeros(4,4);

mass=rho*area*leng;m=mass*diag([1/2 leng^2/78 1/2 leng^2/78]);

end

A rotina modelobeam abaixo apresenta os passos necessários. A viga édiscretizada em três elementos finitos em função das descontinuidades exis-tentes.

89



%-------------------------------------------------------------------------%

%% MODELO FEM - VIGA% Este programa realiza uma análise estrutural de uma viga% de Euler-Bernoulli%% Variáveis de entrada:% nel = número de elementos% E = modulo de elasticidade% rho = densidade% b = lárgura da viga

% h = espessura da viga% L = comprimento total da viga%% Descrição da variaveis internas da função% k = matriz de rigidez do elemento% m = matriz de massa do elemento% K = matriz de rigidez do sistema% M = matriz de massa do sistema% index = vetor contendo os graus de liberdade associado ao elemento% bcdof = vetor contendo graus de liberdade associados as restrições

% bcval = vetor contendo as condiÁoes de restriçãao (valores)% %%-------------------------------------------------------------------------%% Rotina preparado por: SAMUEL - NOVEMBRO de 2006%-------------------------------------------------------------------------%

clc; clear; close all;

%-------------------------------------------------------------------------%% Dados FEM para a viga %%-------------------------------------------------------------------------%

ndof=2; % numero de graus de liberdade por nonnel=2; % numero de nos por elementonel = 3; % numero de elementos usados para discretizar a vigannode=(nnel-1)*nel+1; % numero total de nos no sistemasdof=nnode*ndof; % total de graus de liberdade no sistema

%----------------------------------------------------------------------------%% Propriedades geometricas da viga %

90



%----------------------------------------------------------------------------%

b = 0.144; % base da seção transversal da vigah = 0.399; % altura da seção transversal da vigaI=(b*h^3)/12; % momento de inercia da seção transversal

E = 210e9; % Módulo de elasticidade da viga

% comprimento dos nel elementos (viga não-uniforme)Lb(1)=2; Lb(2)=2; Lb(3)=1;

%rho e area da seção transversal (Seria Usada se o objetivo fosse% obter a matriz de Massa - Interesse neste caso seria realizar uma% análise dinâmicaarea=b*h;rho=1;

%----------------------------------------------------------------------------%

kk=zeros(sdof,sdof); % inicialização da matriz de rigidez da vigamm=zeros(sdof,sdof); % inicialização da matriz de massa da viga

index=zeros(nel*ndof,1); % inicialização do vetor de indices

%----------------------------------------------------------------------------%% Montagem das matrizes de massa e rigidez da viga %%----------------------------------------------------------------------------%

for iel=1:nel % loop para numero total de elementos% extraindo os graus de liberdade associados ao elemento

index=feeldof1(iel,nnel,ndof);% calculo das matrizes de massa e rigidez do elemento

[k,m]=febeam1(E,I,Lb(iel),area,rho,1);

% montando a matriz de rigidez do sistema globalkk=feasmbl1(kk,k,index);

end

%----------------------------------------------------------------------------%% Aplicando as condições de contorno na Viga %%----------------------------------------------------------------------------%

% Restringe dof x1 e x7 (elimine linhas 1 e 7 e colunas 1 e 7)

91



K=kk([2 3 4 5 6 8],[2 3 4 5 6 8]);

% Forças externas atuantes nas coordenadas x3 e x5:% F = [M2 F3 M4 F5 M6 M8];F=[0 -10e3 0 -5e3 0 0]’;

%----------------------------------------------------------------------------%% Calculo do Problema inverso F=Kx --> x=F\K %%----------------------------------------------------------------------------%

% deslocamento e rotações nos nós da viga:

x=K\F;

% vetor de deslocamento e rotação completo na vigaxc =[0 x(1) x(2) x(3) x(4) x(5) 0 x(6)]’;

%--------------------------------------------------------------------------% Calculo das reações de apoio --> FR=kk*xc -->Problema direto:%----------------------------------------------------------------------------%

FR = kk*xc;

%--------------------------------------------------------------------------% Reações de apoio na viga:%--------------------------------------------------------------------------

RA = FR(1);RB = FR(7);

% Verifique que FR(2)=0 e FR(8)=0 - momentos nos vínculos é nulo!

%--------------------------------------------------------------------------

% Calculo do Momento Fletor M(x)=EId^2v/dv=dB/dsxc% Calculo da Cortante V(x)=EId^2v/dv=dB/dsxc%--------------------------------------------------------------------------

% pontos para interpolar em cada elementos = [0:.1:2;

2:.1:4;4:.05:5];

92



x1 = [0 2 4]; % primeira coordenada do primeior nó em cada elemento.

% Lembrando que s=x-x1!

[el,order]=size(s);

disp = []; % deslocamento trasnversal na vigaMf = []; % momento fletor na viga;Cort = []; % cortante na viga;

for iel=1:nel

% extraindo os graus de liberdade associados ao elementoindex=feeldof1(iel,nnel,ndof);d(iel,:)=xc(index); % dofs no elemento

for i = 1:order;

% calculo das funções de forma e de suas derivadas[N,H,D]=funforma(Lb(iel),s(iel,i)-x1(iel));

% deslocamento transversal no elemento

v(iel,i) = N*d(iel,:)’;

% Momento fletor;M(iel,i) = E*I*H*d(iel,:)’;

% CortanteV(iel,i) = E*I*D*d(iel,:)’;

end

disp = [disp v(iel,:)];

Mf = [Mf M(iel,:)];Cort = [Cort V(iel,:)];

end

% Distribuição de Tensão Normal na viga devido a flexão:

Sigma = Mf*(h/2)/I;

93



comp = [s(1,:) s(2,:) s(3,:)];

% comprimento da viga (com interpolação);

figure(1)plot(comp,1e3*disp,’o-’,’linewidth’,1); gridylabel(’Deslocamento transversal da viga v(x) [mm]’,’fontsize’,12,’fontweight’,’xlabel(’Comprimento [m]’,’fontsize’,12,’fontweight’,’bold’);set(gca,’fontsize’,12,’fontweight’,’bold’);saveas(1,’dispb1.eps’,’psc2’) % salvar figura para arquivo eps

% Momento obtido via método das seções

Msc = [7*s(1,:) -3*s(2,:)+20 -8*s(3,:)+40];

figure(2)plot(comp,1e-3*Mf,’o-’,’linewidth’,1); grid ; hold onplot(comp,Msc,’r’)ylabel(’Momento Fletor M(x) [kN.m]’,’fontsize’,12,’fontweight’,’bold’);xlabel(’Comprimento [m]’,’fontsize’,12,’fontweight’,’bold’);legend(’FEM’,’Método das Seções’);set(gca,’fontsize’,12,’fontweight’,’bold’);saveas(2,’momentob1.eps’,’psc2’) % salvar figura para arquivo eps

figure(3)plot(comp,1e-3*Cort,’o-’,’linewidth’,1); gridylabel(’ForÁa Cortante V(x) [kN]’,’fontsize’,12,’fontweight’,’bold’);xlabel(’Comprimento [m]’,’fontsize’,12,’fontweight’,’bold’);set(gca,’fontsize’,12,’fontweight’,’bold’);saveas(3,’cortanteb1.eps’,’psc2’) % salvar figura para arquivo eps

figure(4)plot(comp,1e-6*Sigma,’o-’,’linewidth’,1); gridylabel(’Tens„o Normal \sigma [MPa]’,’fontsize’,12,’fontweight’,’bold’);

xlabel(’Comprimento [m]’,’fontsize’,12,’fontweight’,’bold’);set(gca,’fontsize’,12,’fontweight’,’bold’);saveas(4,’Tensaob1.eps’,’psc2’) % salvar figura para arquivo eps

%--------------------------------------------------------------------------

A rotina modelobeam também usa a rotina funforma para cálculo dasfunções de forma e de suas derivadas que são usadas ao longo do programa16.

16Para cálculo do momento fletor, tensão e cortante.

94



A função funforma é transcrita a seguir.

function [N,H,D]=funforma(L,s);

%-------------------------------------------------------------------% Objetivo:% Calculo das funções de forma Hermitiana e de suas derivadas no% elemento%% Comando% [N,H,D]=funforma(s,L);

%% Vari·veis% N - vetor de funções de forma Hermitiana% H - vetor da derivada segunda do vetor de função de forma N% D - vetor da derivada terceira do vetor de função de forma% s = x-x1 sendo x1 a coordenada do nó 1 do elemento em questão% L - comprimento do elemento%-------------------------------------------------------------------

% Funções de forma Hermitianas Ni(s) no elemento ielN1 = 1-(3*s^2)/L^2 + 2*s^3/L^3;N2 = s - 2*(s^2)/L + s^3/L^2;N3 = 3*s^2/L^2 - (2*s^3)/L^3;N4 = -(s^2)/L + s^3/L^2;

N = [N1 N2 N3 N4];

% d^2Ni(s)/ds^2H1 = -6/L^2 + 12*s/L^3;H2 = -4/L + 6*s/L^2;H3 = 6/L^2 - 12*s/L^3;

H4 = -2/L + 6*s/L^2;

H = [H1 H2 H3 H4];

% d^3Ni(s)/ds^3D1 = 12/L^3;D2 = 6/L^2;D3 = -12/L^3;D4 = 6/L^2;

95



D = [D1 D2 D3 D4];%-------------------------------------------------------------------

As reações de apoio calculadas via FEM foram RA = 7kN e RB = 8kN .Já os deslocamentos e rotações no nós da viga são mostrados na tabela (4.5).

Tab. 4.5: Deslocamentos transversais e rotações na viga.

Nó v (mm) θ (rad/s)1 0 -0.1249

2 -0.1916 -0.03753 -0.1166 0.1004 0.0 0.1249

A fig. (4.13) apresenta o deslocamento transversal v(x) calculado comFEM assumindo a interpolação usando a função Hermitiana e os valoresconhecidos nos 4 nós da viga. Com o uso das derivadas da função de forma(comando funforma fornecido) e do conhecimento do módulo de elasticidadeE e do momento de inércia de área I , o momento fletor M (x), a cortanteV (x) e as tensões normais são calculadas usando o deslocamento v(x). As

figs. (4.14), (4.15) e (4.16) mostram os resultados gráficos desta análise.

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5−0.2

−0.18

−0.16

−0.14

−0.12

−0.1

−0.08

−0.06

−0.04

−0.02

0

D e s l o c a m e n t o t r a n s v e r s a l d a v i g a v ( x ) [ m m ]

Comprimento [m]

Fig. 4.13: Deslocamento transversal da viga solucionado via FEM com fun-ções de interpolação Hermitiana.

Apenas para ilustração a seguir se apresentam os resultados da cortantee momento usando o método das seções para o trecho 0 < x < 2:

96



0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5−2

0

2

4

6

8

10

12

14

M o m e n t o F l e t o r M ( x ) [ k N . m

]

Comprimento [m]

FEM

Método das Seções

Fig. 4.14: Momento fletor M (x) calculado usando FEM.

0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5−8

−6

−4

−2

0

2

4

6

8

F o r ç a C o r t a n t e V ( x ) [ k

N ]

Comprimento [m]

Fig. 4.15: Cortante V (x) calculada usando FEM.

V (x) = 7 kN (4.94)

M (x) = 7x kN.m (4.95)

Já para o trecho 2 < x < 4 tem-se:

97



0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5−0.5

0

0.5

1

1.5

2

2.5

3

3.5

4

T e n s ã o N o r m a l σ [ M

P a ]

Comprimento [m]

Fig. 4.16: Tensão normal σ(x) devido a flexão calculada usando FEM.

V (x) = −3 kN (4.96)

M (x) =

−3x + 20 kN.m (4.97)

(4.98)

Por fim o trecho 4 < x < 5 tem:

V (x) = −8 kN (4.99)

M (x) = −8x + 40 kN.m (4.100)

Da análise deste resultado se observa que a solução via FEM é igual asolução exata deste problema, não sendo necessário nenhum refinamento namalha.

4.2.9 Exemplo de aplicação com carregamento distri-buído

A fig. (4.17) mostra uma viga não-uniforme com um carregamento dis-tribuído em um trecho. Neste exemplo L = 2 m, F = 18 kN, q = 10 kN/m,E = 210 GPa e I = 4 × 10−4 m4.

A ideia é substituir o carregamento distribuído por um carregamentoconcentrado nos nós do elemento 3, conforme sugerido na fig. (4.11). A

98



Fig. 4.17: Viga não-uniforme com carregamento distribuido.

rotina modelbeam2 transcrita a seguir realiza e gera os gráficos desta análise,além do cálculo das reações de apoio na viga.

%-------------------------------------------------------------------------%%% MODELO FEM - VIGA% Este programa realiza uma análise estrutural de uma viga% de Euler-Bernoulli%% Vari·veis de entrada:% nel = número de elementos% E = modulo de elasticidade% rho = densidade% b = lárgura da viga

% h = espessura da viga% L = comprimento total da viga%% Descrição da variaveis internas da função% k = matriz de rigidez do elemento% m = matriz de massa do elemento% K = matriz de rigidez do sistema% M = matriz de massa do sistema% D = matriz de amortecimento do sistema (proporcional)

99



% index = vetor contendo os graus de liberdade associado ao elemento

% bcdof = vetor contendo graus de liberdade associados as restriçõoes% bcval = vetor contendo as condiçõoes de restrição (valores)%%-------------------------------------------------------------------------%% Rotina preparado por: SAMUEL - MARÇO DE 2009%-------------------------------------------------------------------------%

clc; clear;

%-------------------------------------------------------------------------%

% Dados FEM para a viga %%-------------------------------------------------------------------------%

ndof=2; % numero de graus de liberdade por nonnel=2; % numero de nos por elementonel = 3; % numero de elementos usados para discretizar a vigannode=(nnel-1)*nel+1; % numero total de nos no sistemasdof=nnode*ndof; % total de graus de liberdade no sistema

%----------------------------------------------------------------------------%

% Propriedades geometricas da viga%----------------------------------------------------------------------------%

q = -10e3; % carregamento (N/m)I=4e-4; % momento de inercia da seção transversal

E = 210e9; % Módulo de elasticidade da viga

% comprimento dos nel elementos (viga não-uniforme)Lb(1)=2; Lb(2)=2; Lb(3)=2;E(1:2) = 2*210e9;

E(3) = 210e9;

%rho e area da seçãao transversal (Seria Usada se o objetivo fosse% obter a matriz de Massa - Interesse neste caso seria a realização de uma% análise dinâmica.area=1;rho=1;

%----------------------------------------------------------------------------%

100



kk=zeros(sdof,sdof); % inicialização da matriz de rigidez da vigamm=zeros(sdof,sdof); % inicialização da matriz de massa da vigaindex=zeros(nel*ndof,1); % inicialização do vetor de indices

%----------------------------------------------------------------------------%% Montagem das matrizes de massa e rigidez da viga %%----------------------------------------------------------------------------%

for iel=1:nel% loop para numero total de elementos

index=feeldof1(iel,nnel,ndof);% extraindo os graus de liberdade associados ao elemento[k,m]=febeam1(E(iel),I,Lb(iel),area,rho,1);% calculo das matrizes de massa e rigidez do elementokk=feasmbl1(kk,k,index);% montando a matriz de rigidez do sistema global

end

%----------------------------------------------------------------------------%

% Aplicando as condições de contorno na Viga %%----------------------------------------------------------------------------%

% Restringe dof x1, x2 e x7 (elimine linhas 1,2 e 7 e colunas 1,2 e 7)K=kk([3 4 5 6 8],[3 4 5 6 8]);

% Forças externas atuantes nas coordenadas:% Com 3 3lementosF=[-18e3 0 q*Lb(3)/2 q*Lb(3)^2/12 -q*Lb(3)^2/12]’;

%----------------------------------------------------------------------------%

% Calculo do Problema inverso F=Kx --> x=F\K %%----------------------------------------------------------------------------%

% deslocamento e rotações na viga:x=K\F;

% vetor de deslocamento e rotação completo na vigaxc =[0 0 x(1) x(2) x(3) x(4) 0 x(5)]’;

101



%--------------------------------------------------------------------------

% Calculo das reações de apoio --> FR=kk*xc -->Problema direto:%----------------------------------------------------------------------------%

FR = kk*xc;

%--------------------------------------------------------------------------% Reações de apoio na viga:%--------------------------------------------------------------------------

RA = FR(1);

MA = FR(2);RB = FR(7);

%--------------------------------------------------------------------------% Calculo do Momento Fletor M(x)=EId^2v/dv3=dB/dsxc% Calculo da Cortante V(x)=EId^3v/dv3=d^2B/ds^2xc%--------------------------------------------------------------------------

% pontos para interpolar em cada elementos = [0:.1:2;

2:.1:4;4:.1:6];

x1 = [0 2 4]; % primeira coordenada do primeiro nó em cada elemento.% Lembrando que s=x-x1!

[el,order]=size(s);

disp = []; % deslocamento trasnversal na vigaMf = []; % momento fletor na viga;Cort = []; % cortante na viga;

for iel=1:nel

% extraindo os graus de liberdade associados ao elementoindex=feeldof1(iel,nnel,ndof);d(iel,:)=xc(index); % dofs no elemento

for i = 1:order;

102



% calculo das funçõeses de forma e de suas derivadas[N,H,D]=funforma(Lb(iel),s(iel,i)-x1(iel));

aux = 0;if iel==3 % correção com o principio da superposição!

% deslocamento transversal no elementoaux = s(iel,i)-x1(iel);vf(iel,i) = (q/(24*E(iel)*I))*(Lb(iel)-aux)^2*aux^2;v(iel,i) = N*d(iel,:)’+vf(iel,i);

% Momento fletor;

d2vfds2 = (2*Lb(iel)^2-12*aux*Lb(iel)+12*aux^2)*q/(24*E(iel)*I);M(iel,i) = E(iel)*I*H*d(iel,:)’+E(iel)*I*d2vfds2;

% Cortanted3vfds3 = (-12*Lb(iel)+24*aux)*q/(24*E(iel)*I);V(iel,i) = E(iel)*I*D*d(iel,:)’+E(iel)*I*d3vfds3;

elsev(iel,i) = N*d(iel,:)’;

% Momento fletor;M(iel,i) = E(iel)*I*H*d(iel,:)’;

% CortanteV(iel,i) = E(iel)*I*D*d(iel,:)’;

end

end

disp = [disp v(iel,:)];

Mf = [Mf M(iel,:)];Cort = [Cort V(iel,:)];

end

comp = [s(1,:) s(2,:) s(3,:)];% comprimento da viga (com interpolação);

figure(1)

103



plot(comp,1e3*disp,’ro-’,’linewidth’,2); grid on ; hold on

ylabel(’Deslocamento transversal da viga v(x) [mm]’,’fontsize’,12,’fontweight’,’xlabel(’Comprimento [m]’,’fontsize’,12,’fontweight’,’bold’);set(gca,’fontsize’,12,’fontweight’,’bold’);saveas(1,’dispb1.eps’,’psc2’) % salvar figura para arquivo eps

figure(2)plot(comp,1e-3*Mf,’ro-’,’linewidth’,2); grid on; hold onylabel(’Momento Fletor M(x) [kN.m]’,’fontsize’,12,’fontweight’,’bold’);xlabel(’Comprimento [m]’,’fontsize’,12,’fontweight’,’bold’);%legend(’3 elementos’,’7 elementos’);

set(gca,’fontsize’,12,’fontweight’,’bold’);saveas(2,’momentob2.eps’,’psc2’) % salvar figura para arquivo eps

figure(3)plot(comp,1e-3*Cort,’ro-’,’linewidth’,2); grid on; hold onylabel(’ForÁa Cortante V(x) [kN]’,’fontsize’,12,’fontweight’,’bold’);xlabel(’Comprimento [m]’,’fontsize’,12,’fontweight’,’bold’);set(gca,’fontsize’,12,’fontweight’,’bold’);saveas(3,’cortanteb2.eps’,’psc2’) % salvar figura para arquivo eps

%--------------------------------------------------------------------------

As figs. (4.18), (4.19) e (4.20) mostram os gráficos do deslocamentov(x), do momento M (x) e da cortante V (x) ao longo da viga. Porém, umavez que o carregamento distribuído no 3.o elemento foi aproximado por umcarregamento concentrado nos nós deste elemento, há um erro no trecho4 < x < 6. Uma forma de contornar este erro é aproximar um número maiorde elementos neste trecho, por exemplo 5 elementos neste trecho. Porém, umaforma melhor é usando um fator de correção vf (s). As figs. (4.18), (4.19) e(4.20) também apresentam a estimativa assumindo que o trecho 4 < x < 6 é

corrigido pelo fator de carga uniforme17:

vf (s) =q (L − s)2s2

24EI (4.101)

Assim no trecho do elemento o deslocamento transversal será interpoladopor:

17Não será dada a prova aqui desta expressão, mas ela é válida para qualquer carrega-mento contínuo constante ao longo de uma viga e é obtida pela princípio da superposição.

104



v(s) = NT d + vf (s) (4.102)

Já o momento fletor neste trecho será dado por:

M (s) = EI

HT d +

d2vf ds2

(4.103)

e a cortante por:

V (s) = EI

DT d +

d3vf ds3

(4.104)

Na transcrição do programa fornecido é apresentada também estámodificação feita na rotina modelbeam2 para encontrar a solução exata.

0 1 2 3 4 5 6−0.35

−0.3

−0.25

−0.2

−0.15

−0.1

−0.05

0

D e s l o c a m e n t o t r a n s v e r s a l d a v i g a v ( x ) [ m m ]

Comprimento [m]

3 elementos

3 elementos com correção

Fig. 4.18: Deslocamento transversal da viga com carregamento distribuídosolucionado via FEM com funções de interpolação Hermitiana.

4.2.10 Viga com carregamento trapezoidal

Qualquer carregamento pode ser tratado de forma exata da mesma ma-neira feita nos exemplos anteriores. Esta seção exemplifica como realizar esteprocedimento e como encontrar o fator de correção usando a superposição docarregamento no elemento.

105



0 1 2 3 4 5 6−40

−30

−20

−10

0

10

20

M o m e n t o F l e t o r M ( x ) [ k N . m

]

Comprimento [m]

3 elementos


Fig. 4.19: Momento fletor M (x) na viga com carregamento distribuído cal-culado usando FEM.

0 1 2 3 4 5 6−20

−15

−10

−5

0

5

10

15

20

25

F o r ç a C o r t a n t e V ( x ) [ k N ]

Comprimento [m]

3 elementos


Fig. 4.20: Cortante V (x) na viga com carregamento distribuído calculadausando FEM.

Considere um elemento de viga com carregamento linearmente variávelde um valor q 1 no primeiro nó até um valor q 2 no segundo nó, como visto nafig. (4.21). Usando interpolação esta carga pode ser descrita pela função:

106



q (s) = sq 2L

− (s − L)q 1L

(4.105)

Fig. 4.21: Carregamento equivalente devido a uma carga distribuída de formatrapezoidal em um elemento.

Usando a função de interpolação para a viga e integrando q (s) obtém-seo vetor de carregamento equivalente rq neste elemento:

rq = L0

Nqds =

L0

q (s)

1 − 3s2

L2+ 2s3

L3

ds

L

0q (s)s

−2s2

L+ s3

L2 ds L0

q (s)3s2

L2− 2s3

L3 ds L

0q (s)

− s2

L+ s3

L2

ds

(4.106)

que gera o resultado:

rq =

120

L (7q 1 + 3q 2)160

L2 (3q 1 + 2q 2)120

L (3q 1 + 7q 2)− 1

60L2 (2q 1 + 3q 2)

(4.107)

Para usar a superposição, a solução de correção com carga trapezoidaldistribuída q (s) pode ser obtida resolvendo a seguinte equação diferencialgovernante:

EI d4v

ds4− q (s) = 0, 0 < s < L (4.108)

Com as seguintes condições de contorno:

v(0) = 0;dv(0)

dx= 0; v(L) = 0;

dv(L)

dx= 0 (4.109)

107



Integrando quatro vezes e introduzindo as constantes de integração

c1, · · · , c4 temos:

v(s) = c4s3 + c3s2 + c2s + c1 − (1/120)s5(q 1 − q 2) − (1/24)Ls4q 1EI L

(4.110)

Com as condições de contorno nós temos:

c1 = 0; c2 = 0; (4.111)

c4L3 + c3L2 + c2L + c1−

(1/120)L5(q 1 − q 2) − (L4q 1/6)

EI L= 0 (4.112)

3c4L2 + 2c3L + c2 − (1/24)(q 1 − q 2) − (L4q 1/6)

EI L= 0 (4.113)

Solucionando este sistema as constantes de integração são dada por:

c1 = 0; (4.114)

c2 = 0; (4.115)

c3 = −−3q 1L2 − 2q 2L2

120EI (4.116)

c4 = −L(7q 1 + 3q 2)120EI

(4.117)

Assim, a solução corretora para o carregamento trapezoidal distribuídona viga é:

vf (s) =(L − s)2s2((3L − s)q 1 + (2L + s)q 2)

120EI L(4.118)

4.3 Frames

Os elementos de um frame (pórtico) são projetados para resistir deforma-ções axiais e flexão. O elemento de viga com dois nós estudado anteriormentepode ser combinado com o elemento de barra para formar um elemento bas-tante versatil que pode ser utilizado para analisar montagens no plano ouespaço de frames. Este elemento pode ser útil para modelar de forma sim-ples o comportamento dinâmico e estático de mecanismos, máquinas, robôs,rotores, pás de turbinas, paletas e até mesmo de uma aeronave ou uma bar-ragem de concreto. É assumido a priori que os efeitos axiais e de flexãosão desacoplados um do outro, o que é uma condição razoável de se assumir

108



dentro de estruturas sofrendo deformações pequenas. Neste capítulo apenas

a formulação para frames no plano será apresentada, porém a expansão paraelementos no espaço é direta a partir da teoria apresentada a seguir.

Como mostrado na fig. (4.22) um sistema local de coordenadas é esta-belecido para cada elemento. Neste sistema de coordenadas locais o eixo sestá ao longo do elemento axial, com direção positiva do nó 1 para o nó 2. Oeixo local t está a 90o no sentido anti-horário do eixo s. Cada nó tem então3 graus de liberdade, duas translações (em s e t) e uma rotação no plano. Omomento de inércia de área e carregamento distribuído transversalmente sãoconsiderados constantes em cada elemento. Além disto, cargas concentradassão permitidas apenas em cima dos nós de um elemento. A carga distribuída

q s está na direção axial e q t é a carga na direção transversal. Estes carrega-mentos aplicados em um elemento são positivos se eles atuam nas direçõespositivas dos eixos locais.

Fig. 4.22: Elemento de frame no plano.

No sistema de coordenadas local, as equações do elemento são simples-mente uma combinação do elemento de barra e de viga de Euler-Bernoulli.Com a ordem dos graus de liberdade como mostrado na fig. (4.22) as equa-ções locais do elemento de frame são da ordem 6 × 6. Assim:

EAL

0 0 −EAL

0 00 12EI

L36EI L2

0 −12EI L3

6EI L2

0 6EI L2

4EI L

0 −6EI L2

2EI L

−EAL

0 0 EAL

0 00 −12EI

L3−6EI

L20 12EI

L3−6EI

L2

0 6EI L2

2EI L

0 −6EI L2

4EI L

d1d2d3d4d5d6

=

12

q sL12

q tL112

q tL2

12

q sL12

q tL− 1

12q tL

2

(4.119)

A equação anterior pode ser escrita na forma matricial como:

109



kldl = rl (4.120)Assumindo um referência comum global (x, y) para todos os elementos18

no frame, os graus de liberdade nas direções globais são descritos por u1, v1,θ1, u2, v2 e θ2. Para desenvolver a transformação do sistema local (s, t) parao sistema global (x, y) deve-se observar que:

u1 = d1cosα = d1ls (4.121)

v1 = d1cos(90 − α) = d1senα = d1ms (4.122)

sendo ls e ms os ângulos entre os sistemas de coordenadas. Estes ângulospodem ser mais facilmente programados em função das coordenadas de cadaelemento, informação está que está contida na malha da estrutura. Deno-tando as coordenadas dos nós de um elemento como (x1, y1) para o nó 1 e(x2, y2) para o nó 2 tem-se facilmente que:

L =

(x2 − x1)2 + (y2 − y1)2 (4.123)

ls =x2 − x1

L(4.124)

ms =y2

−y1

L (4.125)

Isto é exatamente a mesma transformação feita com o elemento de barrano plano (elemento de treliça). Deve ser observado que a rotação não éafetada pela transformação, assim:

θ1 = d3 (4.126)

θ2 = d6 (4.127)

Assim a transformação dos graus de liberdade do elemento no sistema

(x, t) para o sistema (s, t) pode ser feita por:

d1d2d3d4d5d6

=

ls ms 0 0 0 0−ms ls 0 0 0 0

0 0 1 0 0 00 0 0 ls ms 00 0 0 −ms ls 00 0 0 0 0 1

u1

v1θ1u2

v2θ2

(4.128)

18Uma vez que os elementos de frame podem ser montados com geometria diversas.

110



ou na forma matricial:

dl = Td (4.129)

sendo T a matriz de transformação e d o vetor de graus de liberdade nosistema (x, y). Também é possível fazer uma transformação inversa de (s, t)para (x, y) a partir de:

d = TT ds (4.130)

Usando estas transformações, as equações do elemento no sistema decoordenadas locais pode ser relacionada com o sistema global como segue-se:

kldl = rl ⇒ klTd = rl (4.131)

Multipliando ambos os lados por TT tem-se:

TT klTd = TT rl (4.132)

Note que r = TT rl, assim:

kd = r (4.133)

sendo:

k = TT klT (4.134)

r = TT rl (4.135)

É possível realizar a multiplicação das matrizes de transformação e obtera matriz de rigidez no sistema (x, y) escrevendo a matriz explicitamente nestesistema. Porém isto pode ser feito numericamente no momento de implemen-tação em algum software. A montagem e a solução é feita na forma padrãode FEM, solucionando primeiro o problema inverso (restringindo as coorde-nadas nas condições de contorno) e obtendo as soluções nodais, resolvendo

o problema direto e encontrando as reações nos apoios, e por fim fazendo assoluções em cada elemento.

Uma vez conhecida a solução nodal no sistema (x, y) obtém-se a soluçãono sistema (s, t) a partir da transformação da eq. (4.129). O deslocamentoaxial será dado na coordenada s por:

u(s) =

−s − L

L

s

L

d1d4

, 0 ≤ s ≤ L (4.136)

já o deslocamento transversal em função de s será dado por:

111



v(s) = (N 1(s) N 2(s) N 3(s) N 4(s))

d2d3d5d6

+ vf (s) (4.137)

sendo N i(s) a função de forma para o elemento do tipo Hermitiano:

N 1(s) =2s3

L3− 3s2

L2+ 1 (4.138)

N 2(s) = s3L2

− 2s2L

+ s (4.139)

N 3(s) =3s2

L2− 2s3

L3(4.140)

N 4(s) =s3

L2− s2

L(4.141)

o termo vf (s) é um parâmetro corretor do fato de assumir que um carre-gamento distribuído é aplicado apenas nos nós do elemento. Para o caso docarregamento constante no elemento ele é:

vf (s) = q t(L − s)2s2

24EI ; 0 ≤ s ≤ L (4.142)

Note que é possível este termo ser referente a um carregamento trapezoi-dal, eq. (4.118), ou mesmo um carregamento qualquer. Usando estes termoscorretores a solução exata para um frame necessita apenas um elemento entrecada descontinuidade da estrutura toda. Por fim as forças axiais nas barrasF (s), o momento fletor M (s) e as forças cortantes V (s) são calculadas nodomínio 0 ≤ s ≤ L por:

F (s) = EABT

da (4.143)M (s) = EI

HT db +

d2vf ds2

(4.144)

V (s) = EI

DT db +

d3vf ds3

(4.145)

sendo da os graus de liberdade axiais, B a derivada primeira da função deforma da interpolação linear, db os graus de liberdade transversal e rotação,H e D as derivadas segunda e terceira do vetor função de formas N da

112



interpolação Hermitiana. Note também que para o caso dos dofs transversais

s = x − x1, sendo x1 a coordenada do primeiro nó do elemento em questão.Nas próximas seções são solucionados problemas simples usando o

Matlab R para ilustrar todo o procedimento de análise estrutural de umframe no plano.

4.3.1 Exemplo de aplicação com carregamento concen-trado

O primeiro exemplo consiste em determinar os deslocamentos, momentosfletores e forças cortantes no frame mostrado na fig. (4.23). Os parâmetros

que definem este problema são: M = 20 kN.m, P = 10 kN, L = 1 m,E = 210 GPa, A = 4 × 10−2 m2 e I = 4 × 10−4 m4.

Fig. 4.23: Frame no plano com carregamento concentrado.

A forma mais simples de resolver este problema é usando três elementoscom dois nós e com cada nó contendo três graus de liberdade. Após montar

uma malha, fig. (4.24), calculamos as matrizes de rigidez para cada elementoem coordenadas locais, transformamos para coordenadas globais e realizamosuma montagem do sistema completo de equações aplicando as cargas nosrespectivos nós e dofs que estas atuam.

Por fim, ajustando as condições de contorno essenciais conhecidas, ouseja, restringindo valores dos dofs conhecidos, obtemos a solução nodal comos deslocamentos axiais, tranversais e rotações.

Os 12 dofs nodais neste problema são:

113



−0.5 0 0.5 1 1.5 2 2.5−2

−1.5

−1

−0.5

0

0.5

2−D Frame

Width [m]

H e i g h t [ m ]

Fig. 4.24: Malha do frame.

u1 = 0 v1 = 0 θ1 = 0.0785 × 10−3 (4.146)

u2 = 0 v2 = 0.0686 × 10−3

θ2 = 0.0487 × 10−3

(4.147)u3 = −0.0189 × 10−3 v3 = 0.0703 × 10−3; θ3 = −0.0108 × 10−3 (4.148)

u4 = −0.0189 × 10−3 v4 = 0 θ4 = −0.1596 × 10−3 (4.149)

O cálculo das reações externas seria feito da mesma forma dos exemplosapresentados anteriormente. Já a interpolação para estimativa de u(s), v(s) edos esforços F (s), M (s) e V (s) é calculada elemento a elemento realizando asrespectivas transformações dos dofs globais para locais e usando as funções deforma e suas derivadas em cada elemento. O programa a seguir em Matlab Rmostra os passos necessários para solução computacional deste problema.

%--------------------------------------------------------------------------% This program computes the mass and stiffness matrix% for a frame%% Variable descriptions%% x and y = global x and y coordiates of each node% m = element mass matrix% k = element stiffness matrix

114



% M = system mass matrix

% K = system stiffness matrix% ff = system force vector% index = a vector containing system dofs associated with each element% bcdof = a vector containing dofs associated with boundary conditions% bcval = a vector containing boundary condition values associated with% the dofs in ’bcdof’%--------------------------------------------------------------------------

clc; clear

%--------------------------------------------------------------------------nel=3; % number of elementsnnel=2; % number of nodes per elementndof=3; % number of dofs per nodennode=(nnel-1)*nel+1; % total number of nodes in systemsdof=nnode*ndof; % total system dofs

el=1; % length of element [m]

%--------------------------------------------------------------------------

% nodal coordinates%--------------------------------------------------------------------------

gcoord(1,1)=0; gcoord(1,2)=0; % x, y coord. values of node 1 in terms of tgcoord(2,1)=el; gcoord(2,2)=0; % x, y coord. values of node 2 in terms of tgcoord(3,1)=el; gcoord(3,2)=-el; % x, y coord. values of node 3 in terms of tgcoord(4,1)=2*el; gcoord(4,2)=-el; % x, y coord. values of node 4 in terms of t

%--------------------------------------------------------------------------% nodal connectivity%--------------------------------------------------------------------------

nodes(1,1)=1; nodes(1,2)=2; % element (1) between nodes 1 and 2nodes(2,1)=2; nodes(2,2)=3; % element (2) between nodes 2 and 3nodes(3,1)=3; nodes(3,2)=4; % element (3) between nodes 3 and 4

%--------------------------------------------------------------------------% Geometrical and physical properties%--------------------------------------------------------------------------

115



Es=210e9; % elastic modulus [N/m^2]

area=4e-2; % cross-sectional areaI=4e-4; % moment of inertia of cross-section [m^4]rho=2700; % mass density per volume (dummy value for static analysis) [k

%--------------------------------------------------------------------------

mm=zeros(sdof,sdof); % initialization of system mass matrixkk=zeros(sdof,sdof); % initialization of system stiffness matrixindex=zeros(nel*ndof,1); % initialization of index vector

%--------------------------------------------------------------------------

for iel=1:nel % loop for the total number of elements

nd(1)=nodes(iel,1); % 1st connected node for the (iel)-th elementnd(2)=nodes(iel,2); % 2nd connected node for the (iel)-th element

x1(iel)=gcoord(nd(1),1); y1(iel)=gcoord(nd(1),2); % coordinate of 1st nodex2(iel)=gcoord(nd(2),1); y2(iel)=gcoord(nd(2),2); % coordinate of 2nd node

leng(iel)=sqrt((x2(iel)-x1(iel))^2+(y2(iel)-y1(iel))^2); % element length

if (x2(iel)-x1(iel))==0;beta=2*atan(1); % angle between local and global axes

elsebeta=atan((y2(iel)-y1(iel))/(x2(iel)-x1(iel)));

end

% ls = (x2(iel)-x1(iel))/leng(iel);% beta = acos(ls);

index=feeldof(nd,nnel,ndof); % extract system dofs for the element

[k,m]=feframe2(Es,I,leng(iel),area,rho,beta,1);% compute element stiffness and mass matrix

% Ms=feasmbl1(Ms,m,index); % assemble each element matrix into system mass makk=feasmbl1(kk,k,index); % assemble each element matrix into system stiffnes

end

116



%--------------------------------------------------------------------------% Plot the frame

figure(1)plot(gcoord(:,1),gcoord(:,2),’ro’,’linewidth’,1); hold on;

for iel=1:nelplot([x1(iel) x2(iel)],[y1(iel) y2(iel)],’linewidth’,2)

end

set(gca,’fontsize’,12,’fontweight’,’bold’);title(’2-D Frame’)xlabel(’Width [m]’,’fontsize’,12,’fontweight’,’bold’);ylabel(’Height [m]’,’fontsize’,12,’fontweight’,’bold’);axis([-0.5 2.5 -2 0.5])saveas(1,’2Dframe.eps’,’psc2’);

%--------------------------------------------------------------------------% applied constraints%--------------------------------------------------------------------------

% 1, 2 and 11

K = kk([3 4 5 6 7 8 9 10 12],[3 4 5 6 7 8 9 10 12]);

%--------------------------------------------------------------------------% Esforços aplicados% dof 5 - P = 10kN% dof 12 - M = -20 kN.mP = -10e3;M = -20e3; % POR FAVOR, TOMEM CUIDADO COM OS SINAIS!!

F = [0 0 P 0 0 0 0 0 M]’; % vetor com esforços aplicados

%--------------------------------------------------------------------------% Respostas Nodais

xd = K\F;

% vetor de deslocamentos coordenadas globaisd = [ 0; 0 ;xd(1:8); 0; xd(9)];

117



% calculo das cargas externas aplicadasFR = kk*d;

%--------------------------------------------------------------------------% Calculo das forças axiais, momento fletor e força cortante no frame%--------------------------------------------------------------------------

aux = [0:0.1:1;0:-.1:-1;1:.1:2]; % coordenada aux

xnode1 = [0 1 1]

disp = []; % deslocamento transversal no frameMf = []; % momento fletor no frame;Cort = []; % cortante no frame;

for iel = 1:nel

nd(1)=nodes(iel,1); % 1st connected node for the (iel)-th element

nd(2)=nodes(iel,2); % 2nd connected node for the (iel)-th element





end

% Matriz de Transformação de coordenadas

T=[ cos(beta) sin(beta) 0 0 0 0;...-sin(beta) cos(beta) 0 0 0 0;...

0 0 1 0 0 0;...0 0 0 cos(beta) sin(beta) 0;...0 0 0 -sin(beta) cos(beta) 0;...

118



0 0 0 0 0 1];


% Transformação de global para local

dl = T*d(index);

L = leng(iel); % comprimento do elemento

for j = 1:length(aux)

% Força axial% revise teoria de treliça% Deslocamento axialu(iel,j) = [(L-aux(iel,j))/L aux(iel,j)/L]*dl([1 4]);

% Eps = du/ds = [Bu]*dl;Eps(iel,j) = [-1/L 1/L]*dl([1 4]);; % deformação axial

Sigma(iel,j) = Es*Eps(iel,j); % Tensão Axial

Fax(iel,j) = Sigma(iel,j)*area; % Foça axial na barra

% Deslocamento transversals = aux(iel,j)-xnode1(iel);

% chama os vetores funções de forma e suas derivadas[N,H,D]=funforma(L,s);

% deslocamento transversal no elementov(iel,j) = N*dl([2 3 5 6]);

% Momento fletorM(iel,j) = Es*I*H*dl([2 3 5 6]);

% Cortante

V(iel,j) = Es*I*D*dl([2 3 5 6]);

end

119



disp = [disp v(iel,:)];Mf = [Mf M(iel,:)];Cort = [Cort V(iel,:)];

end

%--------------------------------------------------------------------------

A função que calcula a matriz de rigidez para o elemento é chamada de feframe2.m e apresentada a seguir.

function [k,m]=feframe2(el,xi,leng,area,rho,beta,ipt)

%--------------------------------------------------------------% Purpose:% Stiffness and mass matrices for the 2-d frame element% nodal dof u_1 v_1 theta_1 u_2 v_2 theta_2%% Synopsis:% [k,m]=feframe2(el,xi,leng,area,rho,beta,ipt)%% Variable Description:% k - element stiffness matrix (size of 6x6)% m - element mass matrix (size of 6x6)% el - elastic modulus% xi - second moment of inertia of cross-section% leng - element length% area - area of beam cross-section% rho - mass density (mass per unit volume)% beta - angle between the local and global axes% ipt = 1: consistent mass matrix

% is positive if the local axis is in the ccw direction from% the global axis% ipt = 1 - consistent mass matrix% = 2 - lumped mass matrix% = 3 - diagonal mass matrix%--------------------------------------------------------------------------

% stiffness matrix at the local axis

120



a=el*area/leng;

c=el*xi/(leng^3);kl=[a 0 0 -a 0 0;...

0 12*c 6*leng*c 0 -12* c 6*leng*c;...0 6*leng*c 4*leng^2*c 0 -6*leng*c 2*leng^2*c;...

-a 0 0 a 0 0;...0 -12*c -6*leng*c 0 12*c -6*leng*c;...0 6*leng*c 2*leng^2*c 0 -6*leng*c 4*leng^2*c];

% rotation matrix

r=[ cos(beta) sin(beta) 0 0 0 0;...-sin(beta) cos(beta) 0 0 0 0;...

0 0 1 0 0 0;...0 0 0 cos(beta) sin(beta) 0;...0 0 0 -sin(beta) cos(beta) 0;...0 0 0 0 0 1];

% stiffness matrix at the global axis

k=r’*kl*r;


if ipt==1

mm=rho*area*leng/420;ma=rho*area*leng/6;ml=[2*ma 0 0 ma 0 0;...

0 156*mm 22*leng*mm 0 54*mm -13*leng*mm;...0 22*leng*mm 4*leng^2*mm 0 13*leng*mm -3*leng^2*mm;...

ma 0 0 2*ma 0 0;...0 54*mm 13*leng*mm 0 156*mm -22*leng*mm;...0 -13*leng*mm -3*leng^2*mm 0 -22*leng*mm 4*leng^2*mm];


elseif ipt==2

ml=zeros(6,6);

121



mass=rho*area*leng;

ml=mass*diag([0.5 0.5 0 0.5 0.5 0]);

% diagonal mass matrix

else

ml=zeros(6,6);mass=rho*area*leng;ml=mass*diag([0.5 0.5 leng^2/78 0.5 0.5 leng^2/78]);

end

% mass in the global system

m=r’*ml*r;

Os gráficos de momento e cortante em cada trecho (elemento) poderiamser levantados. Porém, optou-se por mostrar os resultados com uma tabelailustrando os valores dos esforços nos nós de cada elemento, tabela (4.6).

Tab. 4.6: Esforços atuantes nos nós do frame.

Elemento x(m) y(m) F (kN) M (kN.m) V (kN)1 0 0 0 0 -51 1 0 0 -5 -52 1 0 15 -5 02 1 -1 15 -5 03 1 -1 0 -5 -153 2 -1 0 -20 -15

4.3.2 Exemplo de aplicação com carregamento distri-buído

O segundo exemplo considera um frame com carregamento distribuídolinearmente ao longo de um elemento conforme ilustrado na fig. (4.25). Asunidades para este problema estão no sistema inglês: q = 1 kip/ft, L = 15ft, E = 30 × 303 kip/in2, A = 100 in2 e I = 1000 in2.

122



Fig. 4.25: Frame no plano com carregamento distribuído.

Para este frame são necessários apenas dois elementos, um para o membroinclinado e outro para o membro horizontal. Note que no primeiro elementodeverá ser somado ao deslocamento transversal o termo vf (s) para obter asolução exata com apenas um elemento neste trecho, assim como utilizar suas

derivadas segunda e terceira para cálculo exato do momento e da cortanteneste trecho.A solução computacional para este problema é feita com o uso do pro-

grama seguinte.

%--------------------------------------------------------------------------% This program computes the mass and stiffness matrix% for a frame%% Variable descriptions%

% x and y = global x and y coordiates of each node% m = element mass matrix% k = element stiffness matrix% M = system mass matrix% K = system stiffness matrix% ff = system force vector% index = a vector containing system dofs associated with each element% bcdof = a vector containing dofs associated with boundary conditions% bcval = a vector containing boundary condition values associated with

123



% the dofs in ’bcdof’

%--------------------------------------------------------------------------

clc; clear

%--------------------------------------------------------------------------nel=2; % number of elementsnnel=2; % number of nodes per elementndof=3; % number of dofs per nodennode=(nnel-1)*nel+1; % total number of nodes in systemsdof=nnode*ndof; % total system dofs

el=15*12; % length of element [in]

%--------------------------------------------------------------------------% nodal coordinates%--------------------------------------------------------------------------

gcoord(1,1)=0; gcoord(1,2)=0; % x, y coord. values of node 1 in terms of tgcoord(2,1)=el/sqrt(2); gcoord(2,2)=el/sqrt(2); % x, y coord. values of nodgcoord(3,1)=el+el/sqrt(2); gcoord(3,2)=el/sqrt(2); % x, y coord. values of no

%--------------------------------------------------------------------------% nodal connectivity%--------------------------------------------------------------------------

nodes(1,1)=1; nodes(1,2)=2; % element (1) between nodes 1 and 2nodes(2,1)=2; nodes(2,2)=3; % element (2) between nodes 2 and 3

%--------------------------------------------------------------------------% Geometrical and physical properties%--------------------------------------------------------------------------

Es=30000; % elastic modulus [N/m^2]area=100; % cross-sectional areaI=1000; % moment of inertia of cross-section [m^4]rho=2700; % mass density per volume (dummy value for static analysis) [k

q = -1/(12); % kip-in%--------------------------------------------------------------------------

124



mm=zeros(sdof,sdof); % initialization of system mass matrix

kk=zeros(sdof,sdof); % initialization of system stiffness matrixindex=zeros(nel*ndof,1); % initialization of index vector

%--------------------------------------------------------------------------

for iel=1:nel % loop for the total number of elements






end

% ls = (x2(iel)-x1(iel))/leng(iel);% beta = acos(ls);


[k,m]=feframe2(Es,I,leng(iel),area,rho,beta,1); % compute element stiffness

% Ms=feasmbl1(Ms,m,index); % assemble each element matrix into system mass makk=feasmbl1(kk,k,index); % assemble each element matrix into system stiffnes


0 0 1 0 0 0;...0 0 0 cos(beta) sin(beta) 0;...0 0 0 -sin(beta) cos(beta) 0;...0 0 0 0 0 1];

if iel == 1

125



rl = [0 q*leng(iel)/2 q*leng(iel)^2/12 0 ...

q*leng(iel)/2 -q*leng(iel)^2/12 ]’; % esforços locaisF(:,iel) = T’*rl; % esforçoos nas coordenadas globais - soment

else

F(:,iel) = zeros(6,1);

end

end

Forca = [F(:,1)’ 0 0 0]’;

%--------------------------------------------------------------------------% Plot the frame

figure(1)plot(gcoord(:,1),gcoord(:,2),’ro’,’linewidth’,1); hold on;

for iel=1:nelplot([x1(iel) x2(iel)],[y1(iel) y2(iel)],’linewidth’,2)end

set(gca,’fontsize’,12,’fontweight’,’bold’);title(’2-D Frame’)xlabel(’Width [in]’,’fontsize’,12,’fontweight’,’bold’);ylabel(’Height [in]’,’fontsize’,12,’fontweight’,’bold’);saveas(1,’2Dframe1.eps’,’psc2’);

%--------------------------------------------------------------------------

% applied constraints%--------------------------------------------------------------------------% 1, 2, 3, 7,8 e 9

K = kk([4 5 6],[4 5 6]); % Matriz de rigidezFap = F([4 5 6])’; % Forças aplicadas

%--------------------------------------------------------------------------

126



%--------------------------------------------------------------------------

% Respostas Nodais

xd = K\Fap;

% vetor de deslocamentos coordenadas globaisd = [ 0 0 0 xd(1:3)’ 0 0 0]’;

%--------------------------------------------------------------------------% Calculo das forças axiais, momento fletor e força cortante no frame%--------------------------------------------------------------------------

aux = [0:(el/sqrt(2))/100:el/sqrt(2);el/sqrt(2):180/100:(el+el/sqrt(2))]; % coordenada aux

xnode1 = [0 el/sqrt(2)];

disp = []; % deslocamento transversal no frameMf = []; % momento fletor no frame;Cort = []; % cortante no frame;

for iel = 1:nel




if (x2(iel)-x1(iel))==0;

beta=2*atan(1); % angle between local and global axeselse

beta=atan((y2(iel)-y1(iel))/(x2(iel)-x1(iel)));end

% Matriz de Transformação de coordenadas


127



0 0 1 0 0 0;...

0 0 0 cos(beta) sin(beta) 0;...0 0 0 -sin(beta) cos(beta) 0;...0 0 0 0 0 1];


% Transformação de global para local

dl = T*d(index);

L = leng(iel); % comprimento do elemento

for j = 1:length(aux)

% Força axial% Deslocamento axialu(iel,j) = [(L-aux(iel,j))/L aux(iel,j)/L]*dl([1 4]);

% Eps = du/ds = [Bu]*dl;Eps(iel,j) = [-1/L 1/L]*dl([1 4]);; % deformação axial

Sigma(iel,j) = Es*Eps(iel,j); % Tensão Axial

Fax(iel,j) = Sigma(iel,j)*area; % Força axial na barra

% Deslocamento transversals = aux(iel,j)-xnode1(iel);s[N,H,D]=funforma(L,s);

if iel ==1 ; % elemento 1 há correção (carregamento)

% deslocamento transversal no elemento + corrigindo com vf(s)vf(iel,j) = q*((L-s)^2*s^2)/(24*Es*I);v(iel,j) = N*dl([2 3 5 6])+vf(iel,j);

% Momento fletord2vfds2 = (2*L^2-12*s*L+12*s^2)*q/(24*Es*I);M(iel,j) = Es*I*H*dl([2 3 5 6])+Es*I*d2vfds2;

% Cortante

128



d3vfds3 = (-12*L+24*s)*q/(24*Es*I);

V(iel,j) = Es*I*D*dl([2 3 5 6])+Es*I*d3vfds3;

else

v(iel,j) = N*dl([2 3 5 6]);M(iel,j) = Es*I*H*dl([2 3 5 6]);V(iel,j) = Es*I*D*dl([2 3 5 6]);

end

end

disp = [disp v(iel,:)];Mf = [Mf M(iel,:)];Cort = [Cort V(iel,:)];

end

%--------------------------------------------------------------------------

Como no exemplo anterior, os gráficos de momento e cortante em cadatrecho (elemento) poderiam ser levantados. Porém, optou-se por mostrar osresultados com uma tabela ilustrando os valores dos esforços nos nós de cadaelemento, tabela (4.7).

Tab. 4.7: Esforços atuantes nos nós do frame com carregamento distribuído.

Elemento x(in) y(in) F (kip) M (kip.in) V (kip)1 0 0 -7.69 -288.46 -8.511 127.279 127.279 -7.69 -120.89 -2.04

2 127.279 127.279 -10.0268 -105.4 0.8587072 307.279 127.279 -10.0268 49.1988 0.858707

4.4 Considerações finais

Este capítulo apresentou três elementos de extrema importância em aná-lise estrutural: elemento de treliça, viga e frame. Deve ficar claro para os

129



alunos que apesar de termos trabalhado nas aulas com sistemas mecânicos

simples, com estes elementos é possível modelar inúmeros problemas comple-xos de interesse em engenharia mecânica. Alguns exemplos:

• Otimização de formas, geometria e material para atingir requisitos demenor peso, menor preço e melhor desempenho em um sistema de en-genharia.

• Possibilidade de modelar sistemas mecânicos complexos com baixocusto computacional usando elementos de barra, treliça, viga e frame,por exemplo eixos de rotores, turbinas, robôs, estruturas metálicas, etc.

• Análise de falha e proposição de modificação estrutural para atingirnovo desempenho, suportar novas cargas, etc.

• Análise de problemas dinâmicos.

Também esperasse que após este capítulo os alunos tenham adquiridomaior conhecimento em práticas de programação computacional, visto quetodas as soluções apresentadas foram escritas em termos de rotinas em com-putador. Estes programas usam rotinas e funções que são usadas para realizartodo o tipo de trabalho repetitivo e de cálculo. Um engenheiro faz apenaso trabalho de análise e conferência dos valores. Com isto estas rotinas sãogerais. Além disto, as estruturas dos algoritmos numéricos usadas ao longodeste capítulo podem ser bem adaptadas para outros tipos de elementosmais complexos, por exemplo, elementos de viga de Timoshenko, placa deKirchhoff, ou até mesmo elementos sólidos. Uma extensão imediata se referea modelagem de treliças e frames 3D, onde a mudança é feita na matriz detransformação.

Por fim, as soluções computacionais apresentadas neste capítulo foramfeitas usando o Matlab R, porém a solução do problema foi feita usando ométodo dos elementos finitos (FEM). Assim, não obtemos uma solução Ma-tlab, como alguém pode vir a pensar, e sim uma solução via FEM, que seria

a mesma caso se implementasse estas rotinas em qualquer linguagem compu-tacional, como por exemplo, PASCAL, BASIC, SCILAB, C++, FORTRAN,DELPHI, etc. ou até mesmo no Excel, em uma calculadora HP ou ainda emum programa dedicado a FEM, como o ANSYS! Um bom engenheiro quequeira trabalhar com projetos deve dominar bem ao menos uma linguagemde programação para utilização como usuário.

130



4.5 Exercícios

Ex. 4.1 Considere que as barras das treliças das figs. (4.30) a (4.33) sãotodas feitas de aço com módulo de elasticidade de 210 GPa e área da seçãodas barras igual a A = 0.001 m 2. Pede-se:

1. Construa e plote para cada caso uma malha.

2. Calcule as matrizes de rigidez para cada cada elemento de barra e monte as matrizes globais.

3. Aplicando as restrições para cada caso calcule os deslocamento nodais,

forças externas aplicadas (reações de apoio), deformações axiais, ten-sões normais e forças internas axiais nas barras.

4. Implemente sua solução usando as rotinas Matlab R fornecidas.

5. Repita os itens anteriores, porém usando algum software comercial como o Ansys R. Compare e discuta através de tabelas os resultados obtidos com o Matlab R, Ansys R e com os cálculos feitos de forma manual usando, por exemplo, o método dos nós estudado em mecânica geral.

Fig. 4.26: Treliça 1.

Ex. 4.2 Considere as vigas mostradas das figs. (4.34) até (4.41). Todas as vigas são de aço com módulo de elasticidade de E = 210 GPa e momentode inércia I = 4 × 10−4 m 4. Para cada uma destas vigas calcule o deslo-camento transversal, momento fletor e cortante usando FEM. Implemente a solução usando as rotinas fornecidas pelo docente. Compare a solução encon-trada com o método das seções e com a solução obtida com algum programa comercial como Ansys R ou Nastran R.

131




Fig. 4.28: Treliça 3. Considere d = 2.0 m e L = 1000 N.

Fig. 4.29: Treliça 4. Considere a = 3.0 m

Ex. 4.3 Considere os frames das figuras de (4.42) até (4.48). Para cada uma destas estruturas calcule os deslocamentos axiais, transversais, forças

132




axiais, momentos fletores e cortante. Faça os gráficos para cada uma destas variáveis em função do número de elementos usados. Compare a soluçãoFEM obtida com o Matlab R e Ansys R.

133




Fig. 4.32: Treliça 7. α = 30o e P = 500 N.

Ex. 4.4 A seção transversal de uma barragem de concreto é mostrada na fig. (4.49). Usando um modelo simples de frame no plano estime as tensões normais devido a flexão, assim como o diagrama de momento fletor, cortante,

134




Fig. 4.34: Viga 1 com w0 = 10 kN/m, L = 5 m, a = 2 m, b = 4 m e M = 20kNm.

Fig. 4.35: Viga 2 com L = 4 m e w0 = 10 kN/m.

135



Fig. 4.36: Viga 3 com a = 1 m, L = 4 m e w0 = 10 kN/m.

Fig. 4.37: Viga 4 com a = 2 m, b = 3 m, L = 4 m, P = 10 kN e w0 = 10kN/m.

carga axial e os deslocamentos axiais e transversais na barragem devido aopeso próprio da barragem e a pressão d’água agindo na barragem 19. Utilize um número razoável de elementos e calcule as propriedades médias em cada elemento (área e momento de inércia). Assuma que a largura da barragem

é 18 m. A densidade do concreto é 2400 kg/m

3

e da água é 1000 kg/m

3

. O módulo de elasticidade do concreto é 30 GPa.

19Consulte um livro de estática sobre hidroestática.

136



Fig. 4.38: Viga 5 com a = 1 m e w0 = 10 kN/m.

Fig. 4.39: Viga 6.

137



Fig. 4.40: Viga 7.

Fig. 4.41: Viga 8.

138



Fig. 4.42: Frame 1, q = 10 kN/m, L = 2 m, E = 210 GPa, A = 4 × 10

−2

m

2

e I = 4 × 10−4 m4.

Fig. 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.

Fig. 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



Fig. 4.45: Frame 4 com barras quadradas de 200 mm × 200 mm, E = 10GPa e P = 20 kN.

Fig. 4.46: Frame 5.

140



Fig. 4.47: Frame 6.

Fig. 4.48: Frame 7.

Fig. 4.49: Barragem de concreto.

141



Capítulo 5

Análise Dinâmica de Estruturas

via FEM

Uma grande variedade de problemas transientes são de interesse em enge-nharia, como por exemplo nas áreas de mecânica dos fluídos e transferênciade calor. Estes problemas são descritos por equações diferenciais parciais eenvolvem variáveis do tipo u(x , y, z, t) no caso mais geral.

FEM pode ser uma ferramenta bastante versátil para resolver este tipode problema. De maneira geral a ideia é transformar um problema de valorde contorno e com condição inicial em um sistema de equações diferenciaisordinárias lineares que podem ser resolvidas numericamente usando váriosmétodos, por exemplo, algoritmos da família Runge-Kutta ou métodos deintegração implicíta, como o Newmark.

Neste capítulo, em particular, o enfoque será em mecânica estrutural,onde a meta é formular as matrizes de massa para os elementos estudados nocapítulo anterior. No fim, do capítulo uma pequena revisão de conceitos deanálise modal analítica é feita tendo como base o conhecimento das matrizesestruturais de massa e rigidez.

5.1 Matriz de massa para elementos estrutu-rais

Um problema transiente 3D geral em mecânica pode ser descrito pelaseguinte equação diferencial parcial geral:

∂

∂x

kx

∂u

∂x

+

∂

∂y

ky

∂u

∂y

+

∂

∂z

kz

∂u

∂z

+ pu + q = m

∂ 2u

∂t2(5.1)

142



sendo kx, ky, kz, p, q e m funções conhecidas de (x,y,z ). A incógnita a

ser descrita é definida como u sendo dependente de (x,y,z ) e do tempo t.Para solucionar este problema de forma exata é necessário conhecer as con-dições de contorno essenciais e naturais em toda a fronteira e as condiçõesiniciais. Destaca-se que apenas para problemas 3D simples é possível descre-ver uma solução exata para a eq. (5.1). Para os problemas unidimensionais,normalmente a solução exata é facilmente obtida.

O método de Galerkin pode ser aplicado nesta equação a partir da propo-sição de uma função para interpolar u em função de valores conhecidos nosnós. No próximo capítulo isto será detalhado para problemas tridimensionaisa forma geral considerando apenas a parte estática. Podemos adiantar que

a aplicação do método de aproximação de Galerkin na equação diferencialparcial 3D irá envolver o resultado do teorema de Green-Gauss, visto nocapítulo 2 deste texto.

Aqui iremos omitir bastante passagens e concentraremos apenas no resul-tado final. Para um problema dinâmico 3D as equações do elemento finitoformam um sistema de equações diferenciais ordinárias lineares de segundaordem1:

md + (kk + k p + kα) d = rq + rβ (5.2)

sendo m a matriz de massa do elemento, d os graus de liberdade nodais,

d a derivada segunda no tempo dos graus de liberdade nodais, kk, k p e kα

os termos de rigidez e rq e rβ os carregamentos externos e internos aplicadosao sistema.

No caso de um problema de mecânica estrutural a equação do elementoé descrita por:

md + kd = R (5.3)

sendo R o vetor de esforços aplicados, que pode depender dos estadosiniciais de tensão, forças de corpo e forças externas aplicadas; m é a matrizde massa para o elemento descrita por2:

m =

V

ρNNT dV (5.4)

sendo ρ a densidade no elemento, N o vetor da função forma usada e V o volume. A matriz de rigidez para o elemento no caso geral é descrita por:

1Dependendo da natureza física do problema esta equação pode ser de primeira ordem.2Para o caso 3D

143



k = V

BCBT dV (5.5)

sendo C a matriz de propriedades do material3, B a derivada primeira dovetor função de forma.

As matrizes de rigidez k para os elementos mais comuns já foram estuda-das nos capítulos anteriores. A seguir são apresentadas descrições sucintasdas matrizes de massa explicitas para os elementos mais comuns.

5.2 Matriz de massa para barra

Lembrando dos capítulos anteriores que o vetor de funções de forma parao elemento linear de barra envolve N 1(x) e N 2(x) e é dado por:

NT = [N 1(x) N 2(x)] (5.6)

sendo as funções de forma:

N 1(x) = −x − x2

L(5.7)

N 2(x) =

x

−x1

L (5.8)

x1 e x2 são as coordenadas do primeiro e do segundo nó, respectivamente,L é o comprimento do elemento e x a coordenada generalizada.

Usando estas funções de forma a matriz de massa para o elemento serádada pela eq. (5.4):

m =

V

ρNNT dV = ρA

x2x1

N 1(x)N 2(x)

N 1(x) N 2(x)

dx (5.9)

sendo A a área da seção transversal do elemento. Após efetuar a integra-ção obtém-se a matriz de massa para o elemento de barra:

m =ρAL

6

2 11 2

(5.10)

Note que nesta formulação a matriz de massa é distribuída, ou seja, hátermos cruzados entre as coordenadas no dof 1 e dof 2. Uma forma alternativa

3Para o caso de um estado tridimensional ou plano de tensões, a matriz C dá a relaçãoconstitutiva entre o tensor tensão σ e o tensor deformação .

144



de obter a matriz de massa é considerar que no elemento a massa é contínua

e única:

m = ρAL (5.11)

Considerando que a massa seja concentrada e dividida apenas nos nós,pode-se propor que:

m =ρAL

2

1 00 1

(5.12)

Na formulação assumindo massa concentrada o sistema é desacoplado, eq.(5.12), sendo a massa diagonal, já para o caso de assumir massa distribuídaisto não ocorre, eq. (5.10).

A montagem da matriz de massa global M é feita da mesma maneira quea montagem da matriz de rigidez.

5.3 Matriz de massa para elemento de treliça

plana

A matriz de massa para o elemento treliça é obtida a partir das funçõesde interpolação na coordenada s. Assim:

N 1(s) = −s − L

L(5.13)

N 2(s) =s

L(5.14)

Interpolando os dofs:

u

v = N 1 0 N 2 0

0 N 1 0 N 2

u1

v1

u2v2

= NT d (5.15)

Usando estas funções de interpolação, a matriz de massa é calculada por:

m =

V

ρNNT dV = ρA

L0

−s−L

L0

0 − s−LL

sL

00 s

L

ds (5.16)

145



Integrando a equação anterior obtém-se a matriz de massa para o ele-

mento de treliça plana:

m =ρAL

6

2 0 1 00 2 0 11 0 2 00 1 0 2

(5.17)

Não será apresentado aqui, mas a partir da transformação de coordenadaglobal para local, pode-se obter a matriz de massa concentrada para esteelemento:

m = TT mlT (5.18)

5.4 Matriz de massa para elemento de viga

Para a viga de Euler-Bernoulli as funções de forma usada são:

N 1(s) =2s3

L3− 3s2

L2+ 1 (5.19)

N 2(s) =s3

L2 −2s2

L+ s (5.20)

N 3(s) =3s2

L2− 2s3

L3(5.21)

N 4(s) =s3

L2− s2

L(5.22)

Usando estas funções de forma para obtenção da matriz de massa, tem-se:

m = ρA L

0

N 1(s)N 2(s)

N 3(s)N 4(s)

N 1(s) N 2(s) N 3(s) N 4(s) ds (5.23)

Integrando cada termo, obtém-se a matriz de massa distribuída do ele-mento de viga:

m =ρAL

420

156 22L 54 −13L22L 4L2 13L −3L2

54 13L 156 −22L−13L −3L2 −22L 4L2

(5.24)

146



5.5 Matriz de massa para elemento de frame

A matriz de massa para elemento de frame no plano em coordenadaslocais é uma combinação da matriz de massa de viga e barra. Assim:

ml =ρAL

420

140 0 0 70 0 00 156 22L 0 54 −13L0 22L 4L2 0 13L −3L2

70 0 0 140 0 00 54 13L 0 156 −22L0

−13L

−3L 0

−22L 4L2

(5.25)

Lembrando do capítulo anterior que a matriz de transformação de coor-denada global para local é dada por:

T =

ls ms 0 0 0 0−ms ls 0 0 0 0

0 0 1 0 0 00 0 0 ls ms 00 0 0 −ms ls 00 0 0 0 0 1

(5.26)

Usando esta transformação a matriz de massa m para o frame em coor-denadas globais é dada por:

m = TT mlT (5.27)

5.6 Análise modal analítica - resposta livre

Agora que vimos uma forma de obter a equação do movimento de um sis-tema usando o método de elementos finitos devemos resolver este conjuntode EDOs. Uma questão que deve ficar clara de antemão é que o sistema re-

presentado pelas matrizes de massa M e rigidez K corresponde a um sistemaacoplado de equações, o que pode dificultar determinadas análises, além denão podermos fazer uma generalização direta com sistemas mais simples deum grau de liberdade. Nestes casos a transformação do sistema para umaoutra base de coordenadas pode ser muito útil. Dentre as várias coordenadasa do tipo modal ocupa um lugar de destaque em dinâmica estrutural. Estetópico é muitas vezes referenciado na literatura como análise modal analítica.

147



Será apresentad na sequencia, como tratar um problema de vibrações livres

e forçadas como base nos parâmetros modais4.Considere que um sistema global descrito pela Eq. (5.3) com amorteci-

mento nulo C = 0 e não possui excitação externa alguma R = 0. Assim aequação do movimento para o sistema todo se reduz ao seguinte

Mx + Kx = 0. (5.28)

Uma forma de solucionar este problema é propor uma solução do tipo

x = Φe jωt (5.29)

Sendo o vetor Φ formado por amplitudes que indicam quais as formasmodais do problema. Substituindo a Eq. (5.29) dentro da Eq. (5.28) obtém-se

−ω2MΦe jωt + KΦe jωt = 0 (5.30)

que após uma simples manipulação torna-se

e jωt

K − ω2M

Φ = 0 (5.31)

uma vez que e jωt = 0 temos o seguinte problema a solucionar

K − ω2MΦ = 0 (5.32)

A Eq. (5.32) representa um problema clássico de auto-valor e auto-vetor.Este problema também pode ser descrito como

M−1K − λI

Φ = 0 (5.33)

sendo I matriz idêntidade de ordem n × n, sendo n o número de graus deliberdade empregados e λ = ω2. Agora o problema de auto-valor e auto-vetorpode ser escrito em uma forma padrão

M−1KΦ = λΦ (5.34)

AΦ = λΦ (5.35)

sendo A = M−1K. Os auto-valores são dados por λ = ω2 e neste casosão relacionados diretamente as freqüências naturais dos sistemas. Os auto-vetores são dados por Φ e representam os modos de vibrar (ou formas mo-dais), e nada mais são do que uma razão de amplitudes. A solução do pro-

4Freqüência naturais, fatores de amortecimento e modos de vibrar.

148



blema de auto-valor e auto-vetor pode ser feita via inúmeros métodos nu-

méricos, por exemplo o método de Choleski5, porém aqui será resolvido deforma clássica a partir do cálculo do determinante. Assim

det

M−1K − λI

= 0 (5.36)

ou ainda a partir de

det

K − ω2M

= 0 (5.37)

O problema de auto-valor leva à uma equação algébrica em ω2. Comoos coeficientes M e K são normalmente reais e simétricos teremos n raízes

reais, o que implicará em n freqüências naturais. Se o sistema for estável,K é definida positiva e as raízes são positivas. Um sistema não-restringidoapresentará modos de corpo rígido correspondendo a freqüências naturaisnulas.

Importante constatar que os modos de vibrar representam uma base orto-gonal no espaço. Assim a matriz modal Φ apresenta as seguintes propriedadespara i = j

ΦT i MΦ j = 0 (5.38)

ΦT

iKΦ

j =0

(5.39)(5.40)

sendo Φi o i-ésimo modo associado com a i-ésima freqüência natural ωni

e Φ j o j-ésimo modo associado com a j-ésima freqüência natural ωnj . Assim

ΦT i MΦi = 1 (5.41)

ΦT i KΦi = ω2

i (5.42)

Neste caso os modos Φ são normalizados em relação a matriz massa, o

que implica que a matriz modal é ortonormal6. A matriz modal Φ contémas formas de vibrar Φ1 quando o sistema é excitado na primeira freqüêncianatural ωn1, Φ2 quando o sistema é excitado na segunda freqüência naturalωn2 e assim por diante, assim esta matriz é dada por

Φ = [Φ1 Φ2 · · · Φn] (5.43)

5Consulte o help do comando Matlab R eig para maiores informações.6Na maioria das vezes é necessário normalizar a matriz Φ para se ter esta propriedade.

149



Um vez calculados os modos de vibrar e as freqüência naturais pode-se

substituir estes valores na solução proposta, Eq. (5.29) e obter a solução daresposta de vibração do sistema (se conhecendo as condições iniciais x(0)).

O sistema mecânico de MDOF de coordenadas físicas também pode serconvertido em coordenadas modais através da transformação da base físicapara a base modal (representada pela matriz modal Φ) tal que

x = Φq (5.44)

sendo q o vetor deslocamento em coordenadas modais. Substituindo aEq. (5.44) dentro da Eq. (5.28) e pré-multiplicando por ΦT tem-se

ΦT MΦq + ΦT KΦq = 0 (5.45)

Assumindo que a matriz modal Φ é normalizada em relação a matriz demassa M e com a propriedade de ortonormalidade

ΦT MΦ = Mm = I (5.46)

ΦT KΦ = KmΩ (5.47)

sendo Ω = diag(ω21 , ω2

2, · · · , ω2n). Aplicando o resultado acima dentro da

Eq. (5.45) chega-se a equação para o sistema MDOF livre e sem amorteci-mento escrita em uma base modal que é dada por

q + Ωq = 0 (5.48)

Note que a equação acima significa que o sistema de equações diferenciaisordinárias representado pela Eq. (5.28) é totalmente desacoplado em váriossistemas de 1 dof, caso se resolva escreve-lo em coordenadas modais7.

O amortecimento do sistema, no geral, pode ser descrito como sendo dotipo proporcional a massa e rigidez:

C = αM + β K (5.49)

sendo α e β dois valores escalares que são constantes de proporcionali-dade8. A matriz de amortecimento nesta forma também pode ser diagonali-zada via um procedimento modal:

Cm = ΦT CΦ (5.50)

7Ou seja, caso represente-o em outra base ortogonal ou ortonormal (no caso específicoda matriz modal estar normalizada com relação a matriz de massa do sistemas).

8Normalmente ajustadas em procedimentos experimentais.

150



Este modelo de amortecimento normalmente é assumido por conveniência

analítica. A justificativa para isto é devido ao fato da natureza do amorte-cimento não ser exatamente conhecida, sendo que seus valores são muitoaproximados na prática.

A matriz de coeficientes modais pode ser descrita por:

Z = 0.5M−1m CmΩ−1 (5.51)

Assim, o modelo modal para o sistema de 2.o ordem pode ser obtido nasua forma final:

q + 2ZΩq + Ω2q = 0 (5.52)

É importante destacar aos alunos algumas aplicações práticas de análisemodal analítica:

• Para uma estrutura ou sistema existente pode ser difícil se ter a exis-tência de um modelo FEM, nestes casos, deve-se proceder um análisemodal experimental com base nos dados de resposta da estrutura emcondições operacionais ou laboratório9.

• O conhecimento dos modos de vibrar analíticos e experimentais10 po-dem ser usados como métricas (funções objetivo) para indicar a existên-

cia de alguma variação estrutural correspondente a um possível dano(monitorar a saúde da estrutura) ou mesmo para auxiliar algum proce-dimento de ajuste de modelos visando modificação estrutural para atin-gir novo desempenho, minimizar vibrações estruturais, verificar a va-lidade de algum modelo FEM, realizar prognóstico do comportamentodinâmico ou estático de um sistema com variação de carga, etc.

Portanto, análise modal analítica via FEM e análise modal experimen-tal são duas ferramentas complementares amplamente usadas em projetosmodernos de automóveis, aviões, indústria aeroespacial, etc.

5.7 Algoritmo de Newmark

O sistema de equações diferenciais de segunda ordem em mecânica estru-tural pode ser resolvido por qualquer método considerando a existência dealguma excitação F externa sendo aplicado no sistema ou mesmo condiçãoinicial de deslocamento e velocidade em algum nó. Entre estes, o método

9Este tópico é estudado de forma introdutória no curso de vibrações.10E até certo ponto das frequencias naturais e fatores de amortecimento.

151



de Newmark é o mais versátil e popular para solução de grandes sistemas de

equações diferenciais de segunda ordem. Aqui não será dada nenhuma prova.Apenas apresentado sucintamente o método e mostrado um algoritmo efetivopara solução do sistema de EDOs.

Considerando a equação do movimento do sistema descrita pelas matrizesde massa e rigidez obtidas via FEM e com o amortecimento sendo do tipoproporcional a massa e/ou rigidez:

Mx + Cx + Kx = F. (5.53)

sendo x, x e x os vetores aceleração, velocidade e deslocamento, respec-tivamente.

A equação acima pode ser integrada usando algum método numérico.Em essência, a integração numérica direta é baseada em duas ideias. Naprimeira, ao invés de tentar satisfazer a equação acima em todo tempo t,busca-se satisfaze-lá apenas em intervalos discretos de tempo ∆t. A segundaideia consiste em variar os deslocamentos, velocidades e acelerações dentrodo intervalo de tempo ∆t assumido.

Em seguida, considera-se que os vetores deslocamento, velocidade e ace-leração no tempo inicial t0, denotados por x(0), x(0) e x(0) respectivamente,são conhecidos e implementa-se a solução das equações de equilíbrio paraum tempo de t0 até tN . Na solução, todo o tempo considerado é dividido

em N intervalos iguais ∆t(∆t = tN /N ) e o esquema de integração empre-gado estabelece uma solução aproximada para os tempos ∆t, 2∆t, 3∆t, · · ·,t,t + ∆t,· · ·, T N .

O esquema geral de integração de Newmark assume que:

x (t + ∆t) = x(t) + ∆t [(1 − γ )x(t) + γ x(t + ∆t)] (5.54)

x(t + ∆t) = x(t) + ∆tx(t) +

1

2− β

x(t) + β x(t + ∆t)

∆t2 (5.55)

As constantes γ e β são conhecidas como parâmetros de Newmark e sãodeterminados visando obter exatidão e estabilidade numérica. Na literaturaexistem muitas variações deste algoritmo. Newmark originalmente propôs oesquema conhecido como aceleração média constante, conhecida como regratrapezoidal, neste caso γ = 1/2 e β = 1/6. A fig. (5.1) mostra o esquemade integração. Porém outros esquemas podem ser usados, como por exemploγ = 1/2 e β = 1/4, que será empregado na rotina computacional do finaldesta seção.

A ideia é fazer com que a equação do movimento, eq. (5.53), seja válidanos intervalos de tempo de 0 até tN :

152



Fig. 5.1: Esquema de aceleração média constante de Newmark.

Mx(0) + Cx(0) + Kx(0) = F(0)....

Mx(t) + Cx(t) + Kx(t) = F(t).

M¨x

(t + ∆t) +C

˙x

(t + ∆t) +Kx

(t + ∆t) =F

(t + ∆t).... (5.56)

Mx(tN ) + Cx(tN ) + Kx(tN ) = F(tN ).

Com base nesta ideia e no esquema de integração de Newmark pode-se es-crever um algoritmo computacional para integração de equações diferenciaisde segunda ordem de sistemas lineares descrito por quatro passos básicos:

• Inicialização.

•Predição.

• Equação de equilíbrio

• Correção.

Escrevendo explicitamente cada passo temos:

1. Dados do problema: M, C, K

153



2. Inicialização:

x(0) = M−1

F(0) − Cx(0) − Kx(0)

(5.57)

3. Incremento temporal:

tk+1 = tk + ∆t (5.58)

4. Predição:

xtk+1 = xtk + (1 − γ )∆txtk (5.59)

xtk+1 = xtk + ∆txtk + (1

2− β )∆t2xtk (5.60)

5. Equação de equilíbrio:

S = M + γ ∆tC + β ∆t2K (5.61)

xtk+1 = S−1 (Ftk

−Cxtk

−Kxtk) (5.62)

6. Correção:

xtk+1 = xtk + ∆tγ xtk (5.63)

xtk+1 = xtk + ∆t2β xtk (5.64)

7. Critério de parada: atingir tN .

Abaixo é apresentado um pequeno script em Matlab R

que pode ser usado

para integrar a equação do movimento de um sistema mecânico qualquer deforma geral a partir do conhecimento de suas matrizes estruturais, eventualforça de excitação e condições iniciais de velocidade e deslocamento. Estealgoritmo é uma adaptação simples do procedimento apresentado anterior-mente, considerando esquema de integração retangular. Usando pequenos∆t esta rotina é efetiva para aplicações práticas em grandes esquemas deintegração de equações do movimento.

154



%--------------------------------------------------------------------------

%% Classical Newmark Method for time integration considering linear% system%% Date: June 2007, Besançon, France% Author: DA SILVA, Samuel%% function [disp,vel,acc] = newmark_linear(M,K,Damp,Bt,F,t,vel,disp)%% Input:

% M: mass matrix% K: stiffness matrix% Damp: damping matrix% Bt: dof to excitation% F: excitation force% t: time vector% disp: initial displacement vector% vel: initial velocity vector%% Output

% disp: displacement vector% vel: velocity vector% acc: acceleration vector%%--------------------------------------------------------------------------

function [disp,vel,acc] = newmark_linear(M,K,Damp,Bt,F,t,vel,disp)

dt=t(2)-t(1);

% Newmark parameters

gamma = 1/2;Beta = 1/4;

% Initializationacc(:,1)=inv(M)*(Bt*F(1)-Damp*vel(:,1)-K*disp(:,1));

for i =2:length(t);

% Prediction

155



vel(:,i) = vel(:,i-1)+(1-gamma)*dt*acc(:,i-1);

disp(:,i) = disp(:,i-1)+dt*vel(:,i-1)+(.5-Beta)*(dt^2)*acc(:,i-1);

% Equilibrium equationS = M + gamma*dt*Damp + Beta*(dt^2)*K;acc(:,i) = inv(S)*(Bt*F(i)-Damp*vel(:,i)-K*disp(:,i));

% Correctionvel(:,i) = vel(:,i) + dt*gamma*acc(:,i);disp(:,i) = disp(:,i) + (dt^2)*Beta*acc(:,i);

end

%--------------------------------------------------------------------------

5.8 Exemplo de aplicação

Para ilustrar um problema de análise dinâmica irá ser considerada umaviga de alumínio com 0.5 m de comprimento, 25 mm de largura e 5 mm deespessura. Assume-se módulo de elasticidade de 70 GPa e densidade de 2710

kg/m2. A condição de contorno usada será engastada-livre. Discretizandoa viga em 6 elementos, com dois dofs por nó (deslocamento transversal erotação), o sistema restringido terá 12 graus de liberdade. Considere queuma força F atua na extremidade livre na direção do dof x11.

A rotina a seguir apresenta a programação visando extrair analiticamenteos parâmetros modais da estrutura e integrar a equação do movimento dosistema em função de dois tipos diferentes de forças aplicadas. No esquemade integração foi usado uma taxa de amostragem de 3 kHz com 2048 amostrastemporais.

%-------------------------------------------------------------------------%% Análise Dinâmica de uma Viga Engastada-livre% Autor: Samuel da Silva, maio de 2009%-------------------------------------------------------------------------%

clc; clear; close all

%-------------------------------------------------------------------------%% Dados FEM para a viga %

156



%-------------------------------------------------------------------------%

nel = 6; % numero de elementosE = 70e9; % modulo de elasticidade (aluminio: 70 GPa)rho = 2710; % densidade (aluminio 2710)b = 0.025; % largurah = 0.005; % espessuraL = 0.5; % comprimento total da viga

%-------------------------------------------------------------------------%% Matrizes de massa e rigidez linear

%-------------------------------------------------------------------------%

% A rotina viga retorna a matriz de massa e rigidez global% para a viga[M,K] = viga(nel,E,rho,b,h,L);

% Aplicando condições de contorno engastada-livreK=K(3:length(K),3:length(K));M=M(3:length(M),3:length(M));

% Amortecimento proporcional a rigidez e a massaDamp =3e-5*K;

%-------------------------------------------------------------------------%% Parametros para aquisição de dados%-------------------------------------------------------------------------%

N = 2048; % number of samplesFs = 3000; % sampling rate [Hz] - cuidado na escolha deste valordt = inv(Fs); % sampling time [s]t = 0:dt:(N-1)*dt; % time vector [s]

fN=Fs/2; % Nyquist frequencydf=Fs/(N-1); % intervalo entre as frequenciasfreq = 0:df:fN; % frequency vector [Hz]

%-------------------------------------------------------------------------%% Integração da equação do movimento%-------------------------------------------------------------------------%

157



ndof = length(M);

Bt = zeros(ndof,1);Bt(11,1) = 1; %dof que é excitado

% Initial conditionvel(:,1)=zeros(ndof,1); % initial zero velocitydisp(:,1)=zeros(ndof,1); % initial zero displacement

% Dados para integração usando método de Newmark% Força de excitação consideradaFf1=.05*chirp(t,1,1,5000,’q’); % Start 1 Hz, cross 5000Hz at t=1sec

[x1,xdot1,x2dot1]= newmark_linear(M,K,Damp,Bt,Ff1,t,vel,disp);

figure(1)plot(t,Ff1,’linewidth’,1)set(gca,’fontsize’,12,’fontweight’,’bold’);xlabel(’Tempo [s]’,’fontsize’,12,’fontweight’,’bold’);ylabel(’Força [N]’,’fontsize’,12,’fontweight’,’bold’);saveas(1,’fchirp.eps’,’psc2’)

figure(2)plot(t,1e6*x1(11,:),’linewidth’,1)set(gca,’fontsize’,12,’fontweight’,’bold’);xlabel(’Tempo [s]’,’fontsize’,12,’fontweight’,’bold’);ylabel(’Deslocamento [\mum]’,’fontsize’,12,’fontweight’,’bold’);saveas(2,’disp.eps’,’psc2’)

%-------------------------------------------------------------------------%% Força aleatóriaFf2 = 1*randn(N,1);[x1,xdot1,x2dot1]= newmark_linear(M,K,Damp,Bt,Ff2,t,vel,disp);

figure(3)plot(t,Ff2,’linewidth’,1)set(gca,’fontsize’,12,’fontweight’,’bold’);xlabel(’Tempo [s]’,’fontsize’,12,’fontweight’,’bold’);ylabel(’Força [N]’,’fontsize’,12,’fontweight’,’bold’);saveas(3,’frand.eps’,’psc2’)

figure(4)

158



plot(t,1e6*x1(11,:),’linewidth’,1)

set(gca,’fontsize’,12,’fontweight’,’bold’);xlabel(’Tempo [s]’,’fontsize’,12,’fontweight’,’bold’);ylabel(’Deslocamento [\mum]’,’fontsize’,12,’fontweight’,’bold’);saveas(4,’disprand.eps’,’psc2’)

%-------------------------------------------------------------------------%% Análise modal analitica%-------------------------------------------------------------------------%

% Problema de auto-valor e autovetor usando o Método de Choleski

[Phi,Omega]=eig(K,M,’chol’);

Wn = sqrt(abs(Omega)); % frequencia natural em rad/sFn = diag(Wn)/(2*pi); % frequencia natural em Hz

Mm = Phi’*M*Phi; % Matriz de massa modalKm = Phi’*K*Phi; % Matriz de ridiez modalDm = Phi’*Damp*Phi; % Matriz de amortecimento no domÌnio modal

Z = diag(0.5*inv(Mm)*Dm*inv(Wn)); % fatores de amortecimento modal

%-------------------------------------------------------------------------%% Modos de Vibrar

figure(5)subplot(211)plot([0 ;Phi(1:2:11,1)/max(abs(Phi(1:2:11,1)))],’-o’,’linewidth’,2); grid onset(gca,’fontsize’,12,’fontweight’,’bold’);xlabel(’Nós’,’fontsize’,12,’fontweight’,’bold’);ylabel(’1.oModo’,’fontsize’,12,’fontweight’,’bold’);

subplot(212)plot([0 ;Phi(1:2:11,2)/max(abs(Phi(1:2:11,2)))],’-o’,’linewidth’,2); grid onset(gca,’fontsize’,12,’fontweight’,’bold’);xlabel(’Nós’,’fontsize’,12,’fontweight’,’bold’);ylabel(’2.oModo’,’fontsize’,12,’fontweight’,’bold’);saveas(5,’modo12.eps’,’psc2’)

figure(6)subplot(211)

159



plot([0 ;Phi(1:2:11,3)/max(abs(Phi(1:2:11,3)))],’-o’,’linewidth’,2); grid on

set(gca,’fontsize’,12,’fontweight’,’bold’);xlabel(’Nós’,’fontsize’,12,’fontweight’,’bold’);ylabel(’3.oModo’,’fontsize’,12,’fontweight’,’bold’);subplot(212)plot([0 ;Phi(1:2:11,4)/max(abs(Phi(1:2:11,4)))],’-o’,’linewidth’,2); grid onset(gca,’fontsize’,12,’fontweight’,’bold’);xlabel(’Nós’,’fontsize’,12,’fontweight’,’bold’);ylabel(’1.o Modo’,’fontsize’,12,’fontweight’,’bold’);saveas(6,’modo34.eps’,’psc2’)

%-------------------------------------------------------------------------%

A resposta de deslocamento no dof x11 assim como as excitações aplica-das para os dois casos são mostradas nas figuras (5.2) e (5.3). Ressalta-seque no processo de integração da equação do movimento se tem acesso aodeslocamento, velocidade e rotação em todos os dofs do sistema.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7−0.05

−0.04

−0.03

−0.02

−0.01

0

0.01

0.02

0.03

0.04

0.05

Tempo [s]

F o r ç a [ N ]

(a) Força aplicada na extremidade livre.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7−300

−200

−100

0

100

200

300

400

Tempo [s]

D e s l o c a m e n t o [ µ m ]

(b) Deslocamento transversal na extremi-dade livre

Fig. 5.2: Resposta para o sistema.A aplicação do problema de auto-valor e auto-vetor fornece os modos de

vibrar da viga. Deve ficar claro aos alunos que em função da discretizaçãorealizada apenas é possivel obter com acurácia os primeiros modos de vibrar.Na prática, mesmo com malhas bem refinidas, analise dinâmica com FEMsomente fornece bons resultados para faixas de frequencia inferiores a 2 kHz.Esta é uma das grandes desvantagens de FEM, sobretudo para aplicaçõesvisando detectar modificações estruturais minimas que não afetam a estru-tura globalmente (região de baixa frequencia). Para estas situações existe na

160



0 0.1 0.2 0.3 0.4 0.5 0.6 0.7−5

−4

−3

−2

−1

0

1

2

3

4

Tempo [s]

F o r ç a

[ N ]

(a) Força aplicada na extremidade livre.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7−4000

−3000

−2000

−1000

0

1000

2000

3000

4000

Tempo [s]

D e s l o c a m e n t o [ µ m ]

(b) Deslocamento transversal na extremi-dade livre

Fig. 5.3: Resposta para o sistema.

literatura um método relativamente novo de análise, chamado de Método dosElementos Espectrais (SEM), onde a ideia é discretizar com apenas um ele-mento por descontinuidade e interpolar a solução pela própria solução exatano domínio da frequencia. Este método ainda não está popularizado, sobre-tudo por que não existem softwares dedicados. Porém, a sua utilização futurapoderá revolucionar aplicações estruturais por permitir se atingir análise na

faixa de médias e altas frequencias.Para nosso simples exemplo, extraiu-se a partir das matrizes M e K os

quatro primeiros modos de vibrar. As figuras (5.4) e (5.5) ilustram estasformas modais.

A tabela (5.1) apresenta os valores das quatro primeiras frequencias na-turais e fatores de amortecimento para esta viga. Estes valores poderiam sercorrelacionados com parâmetros extraídos experimentalmente visando ajus-tar propriedades da viga ou correlacionar com a saúde estrutural do elemento.

Tab. 5.1: Algumas frequencias naturais e fatores de amortecimento na viga.

Modo Frequencia Natural (Hz) Fator de amortecimento1 16.42 0.00152 102.92 0.00973 288.65 0.02724 568.25 0.0536

161



1 2 3 4 5 6 70

0.2

0.4

0.6

0.8

1

Nós

1 . º

M o d o

1 2 3 4 5 6 7−1

−0.5

0

0.5

1

Nós

2 . º

M o

d o

Fig. 5.4: 1.o e 2.o Modo de vibrar da viga calculado analiticamente.

1 2 3 4 5 6 7−1

−0.5

0

0.5

1

Nós

3 . º

M o d o

1 2 3 4 5 6 7−1

−0.5

0

0.5

1

Nós

1 . º

M o d o

Fig. 5.5: 3.o e 4.o Modo de vibrar da viga calculado analiticamente.

162



Capítulo 6

Elementos Finitos Bidimensionais

Como visto nos capítulos anteriores, a grande parte dos problemas de en-genharia são descritos por equações diferenciais parciais envolvendo um pro-blema de valor de contorno. A solução exata de problemas unidimensionaispode ser obtida em muitos casos, porém geometria e condições de contornocomplexos exigiram a aplicação de um procedimento numérico. Dentre asdiversas abordagens possíveis, o método dos elementos finitos acabou sendoum dos mais populares para aplicações práticas.

Os conceitos fundamentais de FEM válidos para problemas unidimen-sionais também são aplicados em problemas de maior dimensão, como oscasos 2D e 3D. O maior problema é que os procedimentos de integração ediferenciação se tornam um pouco mais complexos, exigindo a aplicação deconceitos de cálculo vetorial mais sofisticados, como por exemplo, a aplicaçãodo teorema de Green-Gauss.

Um problema 2D pode ser escrito pela seguinte equação diferencial par-cial:

∂

∂x

kx

∂u

∂x

+

∂

∂y

ky

∂u

∂y

+ pu + q = 0 (6.1)

que é uma generalização do caso geral em uma dimensão. Nesta equação

kx(x, y), ky(x, y), p(x, y) e q (x, y) são funções conhecidas definidas em umaárea. A solução desconhecida desta equação é u(x, y). Alguns exemplos desistemas que podem ser descritos por uma equação deste tipo:

• Fluxo de potência 2D.

• Fluxo de calor em regime permanente.

• Difusão.

163



•Barras sujeitas a torção.

• Deformação em placas e membranas.

• Inúmeras outras aplicações.

Somente em poucas situações práticas a eq. (6.1) terá uma solução exataobtida de forma analítica. Algumas observações devem ser feitas sobre a eq.(6.1):

1. É uma equação diferencial parcial de segunda ordem.

2. As condições de contorno que envolvem u diretamente são essenciais.3. As condições de contorno que envolvem as derivadas de u são naturais.

4. A área do domínio de solução é denotada por A.

5. O contorno do problema é bidimensional e denotado por C .

6. O contorno C é descrito em termos de uma coordenada c e um vetorunitário normal n.

A fig. (6.1) mostra o domínio de solução de um problema 2D. As com-

ponente x e y do vetor normal ao contorno n são dadas por nx e ny:

n =

nx

ny

(6.2)

|n| =

n2x + n2

y = 1 (6.3)

As condições essenciais são especificadas em uma porção do contorno C indicada por C e e dada por u(c) especificada em C e, onde c é a coordenadaao longo do contorno. Já as condições de contorno naturais são descritas em

uma porção do contorno C indicada por C n e dada pela expressão:

kx

∂u

∂xnx + ky

∂u

∂yny = α(c)u(c) + β (c) em C n (6.4)

sendo α(c) e β (c) parâmetros específicos ao longo do contorno. Quandokx = ky = k a condição natural pode ser escrita como:

k

∂u

∂xnx +

∂u

∂yny

≡ k

∂u

∂n= αu + β (6.5)

164



Fig. 6.1: Domínio de solução de um problema 2D.

6.1 Integração por partes em duas dimensões -

Teorema de Green-Gauss

O estudante deve se recordar do capítulo 3 da forma fraca usando aaproximação proposta por Galerkin, assumindo que u(x, y) seja interpoladapor funções de forma. A solução do problema envolvia integração por partes.Integração por partes em problemas bidimensionais não é uma tarefa dasmais simples e vai exigir a aplicação do teorema de Green-Gauss, que foiapresentado de forma rápida no capítulo 2 desta apostila. A seguir iremosrever os teoremas da divergência e de Green-Gauss, para depois aplicar esteresultado para a forma fraca visando obter as equações FEM usando o métodode Galerkin.

6.1.1 Teorema da divergência de Gauss

Assumindo um vetor de funções F = (F 1(x, y), F 2(x, y))T . O divergentede F é definido como:

div F =

∇ ·F =

∂F 1

∂x

+∂F 2

∂y

(6.6)

de acordo com o teorema da divergência de Gauss, a área sobre a integraldo ∇ · F é igual a integral de linha do produto escalar do vetor F pelo vetorunitário normal n. Assim:

A

∇ · FdA =

c

FT ndc (6.7)

sendo dA = dxdy a área diferencial. Todos os termos na integral decontorno devem ser expressados em termos da condição de contorno c. dc é

165



um elemento diferencial do contorno dado por:

dc =

dx2 + dy2 (6.8)

Escrevendo explicitamente o teorema da divergência de Gauss: A

∂F 1∂x

+∂F 2∂y

dA =

c

(F 1nx + F 2ny) dc (6.9)

sendo nx e ny as componentes do vetor unitário normal ao contorno.

6.1.2 Teorema de Green-Gauss

Caso seja conhecida uma outra função escalar g(x, y), a integral do diver-gente do produto gF pode ser escrita como:

A

div (gF) dA =

c

gFT ndc (6.10)

Usando a regra do produto, o divergente do produto gF pode ser escritocomo:

div (gF) =∂ (gF 1)

∂x+

∂ (gF 2)

∂y=

= g∂F 1∂x

+ F 1∂g

∂x+ g

∂F 2∂y

+ F 2∂g

∂y

Assim a equação anterior pode ser rearranjada da forma:

div(gF) = g

∂F 1∂x

+∂F 2∂y

+

∂g

∂x

∂g

∂y

F 1F 2

(6.11)

Notando que o vetor linha no segundo termo é o transposto do vetorgradiente do campo escalar g:

∇g ≡ ∂g

∂x

∂g

∂yT

(6.12)

Assim:

∇ · (gF) = g(∇ · F) + (∇g)T F (6.13)

Assim a integral do divergente ∇ · (gF) é escrita como:

A

∇· (gF)dA =

A

g(∇·F)dA +

A

(∇g)T FdA =

c

gFT ndc (6.14)

166



Rearranjando estes termos temos uma forma conhecida comumente na

literatura como teorema de Green-Gauss: A

g(∇ · F)dA =

c

gFT ndc −

A

(∇g)T FdA (6.15)

O teorema de Green-Gauss pode ser escrito em uma forma muito similar aconhecida integração por partes no caso unidimensional assumindo o vetor defunções F 1 = f e F 2 = 0, ou seja, um caso particular onde F = (f (x, y), 0)T .O divergente desta função é simplesmente:

∇ ·F =

∂f

∂x

(6.16)

e assim aplicando o teorema de Green-Gauss, tem-se: A

g

∂f

∂x

dA =

c

gf nxdc −

A

∂g

∂xf dA (6.17)

Similarmente fazendo F 1 = 0 e F 2 = f , tem-se: A

g

∂f

∂y

dA =

c

gf nydc −

A

∂g

∂yf dA (6.18)

Note que se as derivadas da função f (x, y) foram feitas com respeito a x,

então o contorno integral envolve nx e vice-versa. De forma muito parecidacom o caso 1D, esta foram do teorema de Green-Gauss fornece uma áreaintegral em que as derivadas da função são trocadas.

Assim, o teorema de Green-Gauss é uma das chaves para obter a formafraca no momento de formular as equações de elementos finitos usando aaproximação de Galerkin. Este procedimento pode ser estendido para casos3D, porém este texto será restrito a aplicações envolvendo domínio 2D.

6.2 Equações de elementos finitos usando o mé-

todo de GalerkinEm um problema 2D ou 3D o domínio de solução é discretizado usando

elementos. Para simplificar a notação não serão introduzidos os subscritose inscritos e assumiremos A e C para ser referir a um elemento arbitráriode área e seu contorno, respectivamente. De início não iremos considerarnenhuma função de forma específica relacionada a nenhum elemento. Nasoutras seções irá se detalhar as funções de forma para elementos retangularese triangulares.

167



O mesmo procedimento adotado para o caso 1D será válido aqui, ou seja,

propor uma solução polinomial envolvendo os graus de liberdade desconhe-cidos dos nós e funções de forma. Obviamente que para os casos 2D ou 3Destas funções não vão ser válidas em todos os pontos do domínio pelo fato dese tratar de uma aproximação. O objetivo é que o erro de usar esta funçãoseja mínimo, daí a necessidade de ponderar por funções de peso N i, que sãoponderações de Galerkin e conhecidas como funções de forma. Todo esteprocesso dá origem a chamada forma fraca de Galerkin.

Assim, o procedimento é similar ao caso 1D, sendo que o resíduo deGalerkin para um problema de valor de contorno 2D é:

A

∂

∂x

kx

∂u

∂x

+

∂

∂y

ky

∂u

∂y

+ pu + q

N idA = 0; i = 1, 2, · · ·

(6.19)No caso 1D usou-se o processo de integração por partes, porém aqui isto

é inviável. Devemos aplicar o teorema de Green-Gauss nos dois primeirostermos do resíduo ponderado de Galerkin, assumindo g = N i e f = kx

∂u∂x

.Assim aplicando as eqs. (6.17) e (6.18):

A∂

∂x kx

∂u

∂xN idA = C kx

∂u

∂x

N inxdc

− A kx

∂u

∂x

∂N i

∂x

dA (6.20) A

∂

∂y

ky

∂u

∂y

N idA =

C

ky

∂u

∂yN inydc −

A

ky

∂u

∂y

∂N i∂y

dA (6.21)

Substituindo os resultados anteriores no resíduo ponderado de Galerkin,eq. (6.19) tem-se:

C

kx

∂u

∂xnx + ky

∂u

∂yny

dc +

A−kx

∂u

∂x

∂N i∂x − ky

∂u

∂y

∂N i∂y + puN i + qN i dA = 0 (6.22)

Neste ponto é muito interessante que o aluno retorne ao capítulo 3 ecompare as eqs. (6.22) e (3.19). As expressão são bastante similares. Aintegral do contorno todo C ele deve ser válida com as condições assumidasessenciais que agem em uma parte do contorno C e e as condições naturaisque atuam em um parte do contorno C n. Assim, é interessante dividir aintegral no contorno neste dois trechos:

168



C e

kx

∂u

∂xnx + ky

∂u

∂yny

N idc +

C n

(αu + β ) N idc

+

A

−kx

∂u

∂x

∂N i∂x

− ky

∂u

∂y

∂N i∂y

+ puN i + qN i

dA = 0 (6.23)

Esta separação em duas integrais no contorno é justificada, pois o alunodeve se atentar ao fato que nas condições essenciais se assumem os valoresconhecidos dos graus de liberdade no contorno. Assim, é de se esperar que oresíduo ponderado de Galerkin seja nulo para este termo. Portanto, aplicando

este resultado na eq. (6.23) a solução assumida da forma fraca equivalenteao problema de valor de contorno bidimensional é dada por:

C n

(αu + β ) N idc

+

A

−kx

∂u

∂x

∂N i∂x

− ky

∂u

∂y

∂N i∂y

+ puN i + qN i

dA = 0, i = 1, 2, · · ·(6.24)

A solução assumida em cada elemento é interpolada por funções de forma2D N i(x, y). Considerando n funções de forma para interpolar u(x, y):

u(x, y) =

N 1(x, y) N 2(x, y) · · · N n(x, y)

u1

u2...

un

= NT d (6.25)

As derivadas parciais de u(x, y) podem ser escritas em função das deriva-das parciais das funções de forma N i(x, y):

∂u

∂x=

∂N 1∂x

∂N 2∂x

· · · ∂N n∂x

u1

u2...

un

= BT x d (6.26)

∂u

∂y=

∂N 1∂y

∂N 2∂y

· · · ∂N n∂y

u1

u2...

un

= BT y d (6.27)

169



sendo u1, u2,

· · ·, un soluções desconhecidas nos nós do elemento. Na

condição de contorno a solução assumida deve ser escrita em termos da co-ordenada de contorno c na forma:

u(c) =

N 1(c) N 2(c) · · · N n(c)

u1

u2...

un

= NT c d (6.28)

Substituindo todos estes resultados na forma fraca eq. (6.24) tem-se:

C n

αN i(c)NT c d + βN i dc

+

A

−kxBT

x d∂N i∂x

− kyBT y d

∂N i∂y

+ pN iNT d + qN i

dA = 0 (6.29)

A eq. (6.29) representa um conjunto de n equações com N i funções deforma 2D, i = 1, 2, · · · , n.

Escrevendo todas a n equações de (6.29) em uma forma matricial:

C n

α

N 1(c)N 2(c)

...N n(c)

NT c d + β

N 1(c)N 2(c)

...N n(c)

dc +

A

−kx

∂N 1/∂x∂N 1/∂x

...∂N n/∂x

BT x d−

dA =

00...0

(6.30)

que pode ser escrita de forma compacta como:

C n

αNcNc

T d + β Nc

dc +

A

−kxBxBT x d − kyByBT

y d + pNNT d + q N

dA = 0 (6.31)

Rearranjando os termos da eq. (6.31) de forma a deixar do lado esquerdotodos os termos envolvendo os dofs do sistema d, obtém-se as equações doelemento:

170



− C n

αNcNT c dc +

A

(kxBxBT x + kyByBT

y − pNNT )dA

d = C n

β Ncdc +

A

q NdA (6.32)

Os dois primeiros termos da integral de área na equação anterior podemser arranjados de forma mais compacta usando matrizes, assim:

kxBxBT x + kyByBT

y = Bx By kx 00 ky

BT x

BT y ≡ BCBT (6.33)

sendo B dada por:

B =

∂N 1∂x

∂N 2∂x

· · · ∂N n∂x

∂N 1∂y

∂N 2∂y

· · · ∂N n∂y

(6.34)

e C uma matriz com as propriedades constitutivas do sistema:

C =

kx 00 ky

(6.35)

Com estes resultados a eq. (6.32) é escrita como:

− C

nαNcNT c dc +

A

(BCBT − pNNT )dA

d =

C n

β Ncdc+

A

q NdA

(6.36)Neste ponto pode-se definir as matrizes ×n:

kk =

A

BCBT dA (6.37)

k p = − A pNNT dA (6.38)

kα = − C n

αNcNT c dc (6.39)

e os vetores n × 1:

rq =

A

q NdA (6.40)

rβ =

C n

β Ncdc (6.41)

171



Assim a equação de elementos finitos fica:

(kk + k p + kα) d = rq + rβ (6.42)

As matrizes kα e rβ são resultantes da aplicação das condições de contornonaturais e irão afetar somente os elementos que tem contorno com C n. Paraos elementos no interior do domínio estas matrizes são nulas. A matriz kk

é a matriz de rigidez do elemento envolvendo as relações constitutivas daspropriedades do sistema, a matriz k p o termo de rigidez devido a aplicaçãodo esforço p e o vetor rq é devido ao carregamento externo aplicado q .

Para escrever explicitamente estas matrizes e vetores é necessário§ seconhecer as funções de forma de cada elemento para efetuar as integraçõesduplas acima. Qualquer forma de elementos finitos pode ser usada e é possíveldescrever qualquer geometria 2D com elementos apropriados. Na sequênciaserá apresentado a forma explicita destas equações envolvendo dois elementos2D típicos: elemento retangular com 4 nós e elemento triangular com 3 nós.

6.3 Elementos finitos retangulares

O elemento retangular é o mais simples que pode ser imaginado paraproblemas bidimensionais. A solução assumida é escrita em função de po-

linômios de forma em x e y. Porém, do ponto de vista de aplicação práticaum elemento retangular não é tão versátil uma vez que pode se tornar muitodificil modelar geometria irregular.

Um típico elemento retangular com quatro nós1 com dimensão 2a por2b é ilustrado na fig. (6.2). O contorno é definido pelos quatro lados doelemento. O vetores unitários normais ao contorno e tangentes também sãoapresentados na fig. (6.2). A análise torna-se um pouco mais fácil em umsistema de coordenadas local descrito por s e t com origem no centro doelemento. Para cada elemento estas coordenadas locais são relacionadas ascoordenadas globais por uma simples translação. Denotando as coordenadasglobais do centroide do elemento por (xc, yc), tem-se as seguintes relaçõesentre as coordenadas locais e globais:

s = x − xc; t = y − yc (6.43)

Notando que com a mudança de variáveis:

ds = dx dt = dy (6.44)

1Também é possível definir elementos retangulares com 6 ou 8 nós.

172



Fig. 6.2: Elemento finito retangular com quatro nós.

ou seja, as derivadas com relação a s são as mesmas com relação a x, assimcomo as derivadas em relação a t serem iguais as derivadas com relação a y.No sistema de coordenadas locais as coordenadas dos nós são:

Nó 1: coordenada (−a,−b).

Nó 2: coordenada (a,−b).

Nó 3: coordenada (a,b).

Nó 4: coordenada (−a,b).

Assumindo a solução desconhecida nos quatro nós como graus de liber-dade do sistema pode-se considerar a solução FEM assumida na área dodomínio A da seguinte forma:

u(s, t) =

N 1(s, t) N 2(s, t) N 3(s, t) N 4(s, t)

u1

u2

u3

u4

= NT d (6.45)

e a solução no contorno como:

u(c) =

N 1(c) N 2(c) N 3(c) N 4(c)

u1

u2

u3

u4

= NT c d (6.46)

173



Em uma área A a solução assumida pode ser escrita com um polinômio em

duas dimensões com quatro coeficientes. Um polinômio linear bidimensionalsatisfatório para um elemento retangular pode ser descrito como:

u(s, t) = c0 + c1s + c2t + c3st =

1 s t st

c0c1c2c3

(6.47)

os coeficientes c0, c1, c2 e c3 pode ser expressos em termos de graus deliberdade nodais pela avaliação dos polinômios nos nós (contorno dos ele-

mentos): u1

u2

u3

u4

=

1 −a −b ab1 a −b −ab1 a b ab1 −a b −ab

c0c1c2c3

(6.48)

invertendo a matriz 4 × 4 tem-se:

c0c1c2c3

=

14

14

14

14− 1

4a14a

14a

− 14a

−1

4b −1

4b

1

4b

1

4b14ab − 1

4ab14ab − 1

4ab

u1

u2

u3

u4

(6.49)

Substituindo estes resultados no polinômios tem-se:

u(s, t) =

1 s t st

14

14

14

14− 1

4a14a

14a

− 14a− 1

4b− 1

4b14b

14b

14ab

− 14ab

14ab

− 14ab

u1

u2

u3

u4

(6.50)

Realizando as multiplicações na equação anterior, obtém-se finalmente asfunções de forma para o elemento retangular de 4 nós:

N 1(s, t) =(a − s)(b − t)

4ab(6.51)

N 2(s, t) =(a + b)(b − t)

4ab(6.52)

N 3(s, t) =(a + s)(b + t)

4ab(6.53)

N 4(s, t) =(a − s)(b + t)

4ab(6.54)

174



Observe que assim como nos casos 1D, a função de forma N i é 1 no i-ésimo

nó e zero nos demais.A solução ao longo dos lados do elemento são escritas substituindo os

valores apropriados de s e t em cada lado. A funções de interpolação nocontorno c para o lado 1 é:

s = c, t = −b =⇒ NT c,1−2 =

a−c2a

a+c2a

0 0

; −a ≤ c ≤ a (6.55)

Lado 2:

s = a, t = c =⇒ N

T

c,2−3 = 0b−c

2b

b+c

2b 0 ; −b ≤ c ≤ b (6.56)Lado 3:

s = −c, t = b =⇒ NT c,3−4 =

0 0 a−c

2aa+c2a

; −a ≤ c ≤ a (6.57)

Lado 4:

s = −a, t = −c =⇒ NT c,4−1 =

b+c2b

0 0 b−c2b

; −b ≤ c ≤ b (6.58)

Após conhecidas as funções de forma do elemento retangular com quatronós deve-se integrar as equações das matrizes que definem este elemento. Porsimplicidade irá se considerar que kx, ky, p, q , α e β são constantes em cadaelemento.

Assim o vetor funções de forma é:

NT =

(a−s)(b−t)4ab

(a+b)(b−t)4ab

(a+s)(b+t)4ab

(a−s)(b+t)4ab

(6.59)

já a matriz derivada do vetor função de forma é:

BT = ∂N i∂s∂N i

∂t = t−b4ab

(b−t)4ab

(b+t)4ab

− (b+t)4ab

s−a

4ab −(a+s)

4ab

(a+s)

4ab

(a−s)

4ab (6.60)

Uma vez que será necessário integrar BCBT :

BCBT =

ky(a−s)2+kx(b−t)2

16a2b2ky(a2−s2)−kx(b−t)2

16a2b2ky(s2−a2)+kx(t2−b2)

16a2b2kx(b2−t2)−ky(a−s)2

16a2b2ky(a+s)2+kx(b−t)2

4a2b2kx(b2−t2)−ky(a+s)2

4a2b2ky(s2−a2)+kx(t2−b2)

4a2b2ky(a+s)2+kx(b+t)2

16a2b2ky(a−s2)−kx(b+t)2

4a2b2

sim. ky(a−s)2+kx(b+t)2

16a2b2

(6.61)

175



Com isto a rigidez kk é:

kk =

a−a

b−b

BCBT dsdt =kya

3b+ kxb

3akya

6b− kxb

3a−kya

6b− kxb

6akxb6a

− kya

3bkya

3b+ kxb

3akxb6a

− kya

3b−kya

6b− kxb

6akya

3b+ kxb

3akya

6b− kxb

3a

sim. kya

3b+ kxb

3a

(6.62)

já o termo k p é:

k p = − a−a

b−b

pNNT dsdt =−4

9abp −2

9abp −1

9abp −2

9abp

−49

abp −29

abp −19

abp−4

9abp −2

9abp

sim. −49

abp

(6.63)

O vetor de carga interna é rq:

rq =

a−a

b−b

q Ndsdt =

abq abq abq abq

(6.64)

Já o vetor rβ e a matriz kα dependem do vetor função de forma Nc

no contorno c. Para o elemento retangular com quatro nós existem quatrocondições de contorno natural (uma em cada face do elemento). Assim, paraa condição de contorno C n no lado 1:

N T c = 1

2 −c2a

c2a +

12 0 (6.65)

kα = − a−a

αNcNT c dc =

−2aα

3−aα

30 0

−2aα3

0 00 0

sim. 0

(6.66)

rβ =

a−a

β NT c dc =

aβ aβ 00

(6.67)

176



No contorno do lado 2 é:

N T c =

0 −12− c

2bc2b

+ 12

0

(6.68)

kα = − b−b

αNcNT c dc =

0 0 0 0

−2bα3

− bα3

0−2bα

30

sim. 0

(6.69)

rβ =

b−b

β NT c dc =

0

bβ bβ

0

(6.70)

No contorno do lado 3 é:

N T c =

0 0 12− c

2ac2a

+ 12

(6.71)

kα = − a−a

αNcNT c dc =

0 0 0 0

0 0 0−2aα

3−aα

3

sim. −2aα3

(6.72)

rβ = a−a

β NT c dc =

0

0aβ aβ

(6.73)

e por fim no contorno do lado 4 é:

N T c =

c2b

+ 12

0 0 12− c

2b

(6.74)

kα = − b−b

αNcNT c dc =

−2bα

30 0 − bα

3

0 0 00 0

sim. −2bα

3

(6.75)

rβ =

b−b

β NT c dc =

bβ 00

bβ

(6.76)

Usando este elemento e seguindo o procedimento padrão de discretização emontagem de equações, uma solução aproximada pode ser obtida para umagrande variedade de problemas práticos que são governados por equaçõesdiferenciais parciais envolvendo domínio 2D.

177



6.4 Exemplo: equação de Laplace em um do-

mínio quadrado

Os problemas básicos de engenharia podem ser descritos por três grandesclasses de equações diferenciais parciais:

Equações Hiperbólicas:∂ 2u

∂t2= a2

∂ 2u

∂x2(6.77)

Este tipo de equação é usado no estudo de vibrações transversais decorda, vibrações longitudinais, oscilações elétricas, vibração torcional

de eixos, etc. Este tipo de equação é muitas vezes denominado deEquação da Onda. Este tipo de equação pode ser bi ou tridimensio-nal.

Equações Parabólicas:∂u

∂t= a2

∂ 2u

∂x2(6.78)

Este tipo de equação é útil para o estudo de propagação de calor, difu-são, problemas de probabilidade, etc. Esta equação é conhecida comoEquação de Fourier e usada em problemas transientes de transmissão

de energia.Equações Elípticas:

∂ 2u

∂x2+

∂ 2u

∂y2= ∇2 · u = 0 (6.79)

Este tipo de equação é usada no estudo de problemas de condução deenergia, campos magnéticos e hidrodinâmica em regime estacionário,ou seja problemas de valor de contorno envolvendo equilíbrio. Estaequação é conhecida como Equação de Laplace.

Este exemplo terá como meta solucionar a equação de Laplace usando

FEM com elementos retangulares. Assuma a seguinte EDP:

∂ 2u

∂ 2x+

∂ 2u

∂ 2y= 0 (6.80)

com domínio:

0 < x < 1; 0 < y < 1 (6.81)

e condições de contorno:

178



Fig. 6.3: Elemento triangular com três nós.

u(0, y) = 0 u(1, y) = 0u(x, 0) = x(1 − x) u(x, 1) = 0

(6.82)

Comparando este problema com a equação geral vista na eq. (6.1) tem-seque kx = ky = 1 e p = q = 0. As condições de contorno nos quatro lados dodomínio são todas essenciais. Uma vez que a equação de Laplace é clássicaconhece-se a solução exata:

u(x, y) =∞n=1

−4sin(nπx) [(−1)n − 1] sinh(nπ(1 − y))

sinh(nπ)n3π3(6.83)

Aqui será obtida uma solução FEM com quatro elementos finitos quadra-dos totalizando ao todo nove graus de liberdade em todo o sistema.

6.5 Elementos finitos triangulares

O elemento finito triangular com 3 nós é apresentado na fig. (6.3). As

coordenadas dos três nós é definida pelos pontos (x1, y1), (x2, y2) e (x3, y3).A origem da coordenada c do contorno é o primeiro nó do elemento.

Assumindo que as soluções desconhecidas nos três nós são dadas por u1,u2 e u3, a solução aproximada FEM tem a seguinte forma no domínio da áreaA da solução:

u(x, y) =

N 1(x, y) N 2(x, y) N 3(x, y) u1

u2

u3

= NT d (6.84)

179



Já no contorno c esta solução é descrita por:

u(c) =

N 1(c) N 2(c) N 3(c) u1

u2

u3

= NT c d (6.85)

Uma solução polinomial linear completa para este elemento 2D têm trêstermos e pode ser descrita por:

u(x, y) = c0 + c1x + c2y =

1 x y

c0c1c2

(6.86)

Os coeficientes desconhecidos c0, c1 e c2 podem ser escritos em termos decoordenadas nodais avaliando os dofs nos nós. Assim: u1

u2

u3

=

1 x1 y11 x2 y21 x3 y3

c0c1c2

(6.87)

Escrevendo a expressão anterior de forma inversa tem-se:

c0c1

c2

=

1

2A

f 1 f 2 f 3b1 b2 b3

c1 c2 c3

u1

u2

u3

(6.88)

sendo A a área do triângulo e as outras constantes escritas em função dascoordenadas dos nós deste elemento:

f 1 = x2y3 − x3y2; f 2 = x3y1 − x1y3; f 3 = x1y2 − x2y1;

b1 = y2 − y3; b2 = y3 − y1; b3 = y1 − y2;

c1 = x3 − x2; c2 = x1 − x3; c3 = x2 − x1;

A =1

2(f 1 + f 2 + f 3)

Assim substituindo estes valores no polinômio que interpola a soluçãou(x, y) obtém-se:

u(x, y)Ψleft( 1 x y1

2A

f 1 f 2 f 3b1 b2 b3c1 c2 c3

u1

u2

u3

(6.89)

Realizando as multiplicações obtém-se as três funções de forma do ele-mento triângular:

180



N 1(x, y) =1

2A(xb1 + yc1 + f 1) (6.90)

N 2(x, y) =1

2A(xb2 + yc2 + f 2) (6.91)

N 3(x, y) =1

2A(xb3 + yc3 + f 3) (6.92)

Notando que a função de interpolação N i é 1 no i-ésimo nó e 0 nos outrosnós. As soluções nos três lados que definem o contorno são descritas porestas funções. Assim, para o lado 1:

u(c) = − c−L12

L12

cL12

0 u1

u2

u3

; 0 ≤ c ≤ L12 (6.93)

sendo L12 =

(x2 − x1)2 + (y2 − y1)2 o comprimento do lado que vai donó 1 até o nó 2.

Para o lado 2 é:

u(c) =

0 − c−L23

L23

cL23

u1

u2

u3

; 0 ≤ c ≤ L23 (6.94)

sendo L23 =

(x3 − x2)2 + (y3 − y2)2 o comprimento do lado que vai donó 2 até o nó 3.

Por fim, o lado 3 é dado por:

u(c) =

cL31

0 − c−L31L31

u1

u2

u3

; 0 ≤ c ≤ L31 (6.95)

sendo L31 =

(x3 − x1)2 + (y3 − y1)2 o comprimento do lado que vai do

nó 3 até o nó 1.Uma vez conhecida as funções de forma para o elemento triângular pode-

se realizar as derivações em relação a x e y:

BT =

∂N 1∂x

∂N 2∂x

∂N 3∂x

∂N 1∂y

∂N 2∂y

∂N 3∂y

(6.96)

com isto a matriz de rigidez kk é escrita como:

181



kk =

A

BCBT dA = ABCBT = kxb21 + kyc21 kxb1b2 + kyc1c2 kxb1b3 + kyc1c3kxb22 + kyc22 kxb2b3 + kyc2c3

sim. kxb23 + kyc23

(6.97)

é interessante notar que uma vez que a matriz B é constante a integraçãoanterior é trivial. Porém, a integração para obter a matriz k p envolve oproduto das funções de forma N que são funções de x e y e devem serintegrados em uma região triangular. Esta integração não é trivial e deve serrealizada com cuidado. Inúmeras referência básicas de FEM apresentam osprocedimentos para obtenção desta matriz. Aqui irá ser apenas fornecido oresultado final deste procedimento:

k p = −

A

pNNT dA =

− pA

12

2 1 1

2 1sim. 2

(6.98)

Já o vetor rq é dado por:

rq =

A

q NdA =qA

3

111

(6.99)

As integrais no contorno são avaliadas separadamente em cada um dostrês lados do elemento triangular.

Assim o lado 1 é:

kα = − C n

αNcNT c dc = −αL12

6

2 1 02 0

sim. 0

(6.100)

rβ =

C n

β NT c dc =

βL12

2

110

(6.101)

O lado 2 é:

182



kα = − C n


6

0 0 02 1 sim.2

(6.102)

rβ =

C n

β NT c dc =

βL23

2

011

(6.103)

E por fim o lado 3:

kα = − C n


6

2 0 10 0

sim. 2

(6.104)

rβ =

C n

β NT c dc =

βL31

2

101

(6.105)

Uma vez definidas todas estas matrizes e vetores pode-se montar a equa-ção do elemento:

(kk + k p + kα) d = rq + rβ (6.106)A montagem e a solução do sistema completo é realizado de forma con-

vencional e igual aos casos 1D estudados até aqui.

183



Capítulo 7

Introdução à Mecânica dos

Sólidos Computacional

Grande parte do projeto de estruturas aeronáuticas, veículos terrestres emarítimos, máquinas e infra-estrutura se baseia na determinação dos cam-pos de tensão e deformação envolvendo sólidos elásticos. As exigências nosprojetos modernos envolvem escolha de materiais com características espe-ciais, geometria complexa, carregamento variável, etc. que necessitam seranalisados no momento do projeto a partir de inúmeras simulações do com-portamento estático e dinâmico do sistema frente a estas cargas.

O uso de computadores revolucionou o conceito de projeto mecâ-nico/estrutural. Até a década de 60/70 grande parte do projeto se base-ava em seguir um roteiro de cálculo do tipo inline . Isto fica evidente casoo leitor busque ler um livro didático desta época, onde se encontrava emanexo inúmeras tabelas, equações prontas, métodos gráficos, etc. Todo esteprocedimento era exaustivo, sem grande precisão e mudanças e ajustes eramcomplicados de serem realizados. Com a popularização dos computadorese das estações de trabalho foi possível também ficar factível a análise deaspectos pouco explorados até então. Com o uso do computador ficou fácilanalisar problemas assumindo elementos sólidos, inúmeros graus de liberdade

e dimensão e empregando técnicas de otimização visando obter formas, dis-tribuição e geometria ótima de estruturas. Também é importante destacaras facilidades gráficas para visualizar simulações envolvendo técnicas avan-çadas de realidade virtual. Os resultados disto foram que o tempo entre umanteprojeto e a construção de um protótipo foi muito encurtado, o que repre-sentou economia de recursos materiais, financeiros e humanos. Um excelenteexemplo é o projeto nacional dos jatos comerciais da família 170/190 da Em-braer. O projeto do EMB 170 foi anunciado em 1999 sendo que o protótipovoôu em 2000. Seria impossível imaginar um projeto desta dimensão sem o

184



Fig. 7.1: Simulação em realidade virtual de ingestão de água nos motores deum avião.

uso de ferramentas modernas de análise. Com esta visão, muitos cursos deengenharia mecânica no mundo e no Brasil se adaptaram bem as mudanças.Toda esta área ficou conhecida como Mecânica Computacional, no caso maisgeral onde o interesse é análise de sólido e fluído, ou mecânica dos sólidoscomputacional no caso mais específico.

Neste sentido este capítulo tem como meta fornecer uma pequena intro-dução neste tipo de concepção de análise, onde a utilização de FEM é indis-pensável e sine qua non . O foco aqui será primeiro rever alguns conceitosbásicos de elasticidade e das relações constitutivas de materiais metálicos.Na sequencia serão apresentadas as equações FEM básicas envolvidas emproblemas no estado plano de tensão. Em especial será detalhada a soluçãousando o elemento triangular 2D.

7.1 Revisão básica de elasticidade

7.1.1 O conceito de tensãoQuando um corpo é submetido a cargas externas forças internas normais

e tangencias são produzidas. A melhor forma para caracterizar estas forças éusar o conceito de tensão. O vetor tensão t em um ponto P qualquer em umadireção n normal ao plano de ação desta tensão é definido matematicamentecom o conceito de limite:

tn(P ) = lim∆A→0∆F

∆A(7.1)

185



Fig. 7.2: Teste em vôo em condição real.

Decompondo a força ∆F em componentes normais ∆Fn e tangenciais Ft

pode-se descrever a tensão normal e tangencial, respectivamente:

σn(P ) = lim∆A→0∆Fn

∆A(7.2)

τ n(P ) = lim∆A→

0

∆Ft

∆A (7.3)

Tensão é definida essencialmente em função do plano que atravessa umponto na estrutura. Como infinitos planos passam por um ponto no es-paço é comum definir as tensões nos planos formados pelo sistema cartesiano(x,y,z ).

Nesta concepção o estado completo de tensões em um ponto P é descritocompletamente em função de nove componentes.

Vetores tensão no plano normal com a direção x:

tx = (σx τ xy τ xz)T (7.4)

Vetores tensão no plano normal com a direção y:

ty = (τ yx σy τ yz)T (7.5)

Vetores tensão no plano normal com a direção z :

tz = (τ zx τ zy σz)T (7.6)

186



Em termos gráficos as componentes de tensão formam um bloco diferen-

cial. Os vetores que descrevem completamente a tensão em ponto podem serescritos como um tensor tensão S:

S =

σx τ xy τ xzτ yx σy τ yzτ zx τ zy σz

(7.7)

Se for considerado uma condição de equilíbrio nas coordenadas retangu-lares as componentes de cisalhamento são iguais:

τ xy = τ yx τ yz = τ zy τ xz = τ zx (7.8)

o que significa que o tensor S é simétrico e pode ser descrito completa-mente por seis componentes:

S =

σx σy σz τ xy τ yz τ xzT

(7.9)

O tensor tensãoS pode ser projetado em qualquer plano normal a qualquervetor unitário normal n = (nx ny nz)T no espaço. Esta projeção é feitapor:

tn = Sn (7.10)

sendo que tn é o tensor S projetado no plano normal ao vetor unitário n.

7.1.2 Tensão principal e direção principal

Deve ficar bem claro aos estudantes que o tensor tensão em um pontoé descrito em função dos planos que atravessam este ponto. Uma vez queinfinitos planos atravessam um ponto é interessante analisar o plano ondeonde as tensão normais são máximas e não há cisalhamento, ou seja, encon-trar um vetor unitário normal n p tal que a projeção do tensor S no planoortogonal a este plano seja diagonal. Esta projeção ortogonal são as tensõesprincipais e o as componentes do vetor unitário n p as direções principais. Isto

pode ser escrito como um problema de auto-valor e auto-vetor:

(S − σiI) np = 0 (7.11)

sendo I a matriz identidade 3 × 3, σi a i-esima tensão principal (auto-valor) e n p um vetor unitário normal ao plano onde a projeção de S forneceas direções principais. O vetor n pi é um auto-vetor, assim para cada tensãoprincipal (auto-valor) existe um vetor unitário n p, ou seja, um auto-vetorcorrespondente. As três tensões principais podem ser ordenadas do maiorpara o menor.

187



Máxima: σ1

Intermediária: σ2

Mínima: σ3

Pode ser mostrado que a máxima tensão cisalhante em qualquer plano édada por:

τ max =σ1 − σ3

2(7.12)

o que significa que τ max age em um plano inclinado 45o com as direções

principais de σ1 e σ2.

7.1.3 Critérios de falhas

Os projetos mecânicos são feitos a partir da escolha de geometria e ma-terial tal que o sistema estrutural/máquina não alcance os valores limitantesmáximos permitidos para que o sistema opere com segurança. Três critériossão comumente usados: Critério de Tresca, Von Mises e critério de Mohr-Coulomb. Abaixo uma revisão rápida de cada um deles. Os alunos devemrelambrar que estes critérios são conservativos.

Critério de Tresca

Neste critério a falha de um sistema ocorre quando a máxima tensãoem qualquer plano em um material alcança o valor da tensão máxima decisalhamento em um teste uniaxial, ou seja, quando σ3 = 0. Assim, conformea eq. (7.12):

τ max =σ1 − σ3

2=

σf

2(7.13)

Com este critério o fator de segurança seg é dado por:

seg =σf /2

τ max

=σf

σ1 − σ3

(7.14)

Note que σf pode ser obtido para um material em um teste experimentaluniaxial. Em alguns textos σf é chamado de tensão admissível.

188



Critério de Von Mises

Este é o critério mais usado em estruturas metálicas. Ele assume que omaterial falha quando a tensão de cisalhamento octaetral é igual a tensãooctaedral em um teste uniaxial. A primeira coisa é relembrar o que é umplano octaedral? Basicamente é um plano ortogonal a um vetor unitárionormal n que tem todos os ângulos iguais com as direções principais. Astensões de cisalhamento projetadas neste plano são conhecidas como tensõesde cisalhamento octaedrais.

A tensão octaedral normal é dada por:

σoct =

σ1 + σ2 + σ3

3 (7.15)Já a tensão de cisalhamento octaedral é:

τ oct =1

3

(σ1 − σ2)2 + (σ2 − σ3)2 + (σ3 − σ1)2 (7.16)

sendo σ1 > σ2 > σ3 as tensões principais. As tensão de cisamelhamentooctaedral pode ser projetada no sistema retangular (x,y,z ) em termos dotensor σ assim:

τ oct =

1

3 (σx − σy)

2

+ (σy − σz)

2

+ (σz − σx)

2

+ 6 τ 2

xy + τ 2

yz + τ 2

zx (7.17)Com esta equação a tensão de cisalhamento octaedral em um teste unixial

é:

τ oct1 =1

3

2σ2

f =

√ 2

3σf (7.18)

Assim, de acordo com a definição do critério de Von Misses a falha domaterial ocorre se:

1√ 2 (σx − σy)2 + (σy − σz)2 + (σz − σx)2 + 6 τ 2xy + τ 2yz + τ 2zx = σf

(7.19)O lado esquerdo da equação anterior é chamado de tensão efetiva σe.

Portanto o coeficiente de segurança para este critério é:

seg =σf

σe

(7.20)

Por fim, salienta-se que o critério de Von Misses é menos conservativo doque o critério de Tresca.

189



Critério de Mohr-Coulomb

O critério de Mohr-Coulomb se aplica a materiais frágeis que tem ten-dência a falhar por tração e compressão com valores diferentes. Com estecritério a falha de um material poder ser prevista quando:

σ1

σtf

− σ3

σcf

= 1 (7.21)

sendo σ1 a tensão principal máxima, σ3 a tensão principal mínima, σtf éa tensão que o material falha por tensão uniaxial e σcf a tensão compressivaque o material falha em um teste uniaxial. Neste critério, o fator de segurançaé definido como:

seg =1

σ1/σtf − σ3/σcf

(7.22)

7.1.4 Relações constitutivas

Quando um corpo sólido elástico sofre um esforço externo ou um carrega-mento aplicado ele gera um campo de deslocamentos espaciais definidos poru(x,y,z ), v(x,y,z ) e w(x,y,z ). No caso mais geral este campo de desloca-mento é variável no tempo, mas aqui iremos nos atentar ao caso de equilíbrioestático. O interesse é saber qual a taxa de variação destes deslocamentos nasdireções cartesianas, o que leva ao conceito de deformação que deve ter sidoestudado em detalhes em um curso formal de resistência dos materiais. O ve-tor deformação completo no espaço considerando pequenos deslocamentos1

é descrito por:

=

xyz

γ xyγ yz

γ zx

=

∂u∂x∂v∂y∂w∂z

∂u∂y

+ ∂v∂x

∂v∂z

+ ∂w∂y

∂u∂z + ∂w

∂x

(7.23)

As tensões e deformações são relacionadas entre si através das equaçõesconstitutivas do material em questão. A forma mais simples de fazer estarelação é generalizar a lei de Hooke do caso unidimensional, assim:

1Neste caso derivadas de ordem superior podem ser desprezadas sem prejuízo da análise.

190



σx

σy

σz

τ xyτ yzτ zx

=

C 11 C 12 C 13 C 14 C 15 C 16C 21 C 22 C 23 C 24 C 25 C 26C 31 C 32 C 33 C 34 C 35 C 36C 41 C 42 C 43 C 44 C 45 C 46C 51 C 52 C 53 C 54 C 55 C 56C 61 C 62 C 63 C 64 C 65 C 66

xyz

γ xyγ yzγ zx

(7.24)

A equação anterior é a base para toda a teoria de elasticidade linear,sendo C ij são constantes elásticas estabelecidas para cada tipo de material.No caso mais geral são 36 constantes. Escrevendo a lei de Hooke generalizada

matricialmente:

σ = C (7.25)

Várias simplificações podem ser assumidas para reduzir o número de cons-tantes elásticas, grande parte assumindo simetria das propriedades em rela-ção aos planos envolvidos. Uma destas simplificações assume que as proprie-dades do material são simétricas em relação aos três eixos de planos coorde-nados, assim o número de constantes pode ser reduzido para 9. Tais materiaissão chamados de ortotrópicos e são descritos pela relação constitutiva:

σx

σy

σz

τ xyτ yzτ zx

=

C 11 C 12 C 13 0 0 0C 21 C 22 C 23 0 0 0C 31 C 32 C 33 0 0 0

0 0 0 C 44 0 00 0 0 0 C 55 00 0 0 0 0 C 66

xyz

γ xyγ yzγ zx

(7.26)

Notando que neste tipo de material as tensão de cisalhamento são com-pletamente definidas apenas pelas deformações transversais e de forma de-sacoplada. Por fim, a forma mais simples ainda é considerar que todas aspropriedades são simétricas e iguais. Este é o caso de materiais chamados deisotrópicos. Para definir completamente estas relações são necessárias apenasduas constantes: E o módulo elasticidade e o coeficiente de Poisson ν :

191



σx

σy

σz

τ xyτ yzτ zx

=E

(1 + ν )(1 − 2ν )

1 − ν ν ν 0 0 0ν 1 − ν ν 0 0 0ν ν 1 − ν 0 0 00 0 0 1−2ν

20 0

0 0 0 0 1−2ν 2

00 0 0 0 0 1−2ν

2

xyz

γ xyγ yzγ zx

(7.27)

O efeito inverso pode ser obtido invertendo a matriz C:

= Dσ (7.28)

sendo a matriz D para os materiais isotrópicos definida por:

D =

1E

− ν E

− ν E

0 0 0− ν

E 1E

−E

0 0 0− ν

E − ν

E 1E

0 0 00 0 0 1

G0 0

0 0 0 0 1G

00 0 0 0 0 1

G

(7.29)

sendo que G é conhecido com módulo de cisalhamento transversal:

G =E

2(1 + ν )(7.30)

Outra consideração importante muitas vezes feita para simplificar é assu-mir que o material é perfeitamente homogêneo, ou seja, E e ν são constantesem todo o corpo sólido.

7.2 Equação de elementos finitos em sólidos

elásticos

Aplicando as condições de equilíbrio de um corpo frente a forças de corpoe campo pode-se obter um sistema acoplado de equações diferenciais parciaisque descreve o campo de tensões σij de um corpo sólido no espaço. Nãoserá fornecido os detalhes para obtenção deste sistema. Qualquer bom livrode elasticidade descreve em detalhes toda a formulação para se chegar aoseguinte problema de equilíbrio:

192



∂σx

∂x+

∂τ xy∂y

+∂τ xz∂z

+ bx = 0 (7.31)

∂τ yx∂x

+∂σy

∂y+

∂τ yz∂z

+ by = 0 (7.32)

∂τ zx∂x

+∂τ xy∂y

+∂σz

∂z + bz = 0 (7.33)

(7.34)

sendo bx, by e bz as forças de corpos atuantes que formam o vetor b. Estasequações podem ser escritas de forma compacta usando notação tensorial:

∂σij

∂i+ bi = 0, i = x,y,z ; j = x,y,z (7.35)

Uma vez que o sistema está acoplado são necessárias três funções deponderação para cálcular o resíduo ponderado de Galerkin, denotados por:u,v e w:

V

∂σx

∂x+

∂τ xy∂y

+∂τ xz∂z

+ bx

u +

∂τ yx∂x

+∂σy

∂y+

∂τ yz∂z

+ by

v +

∂τ zx∂x

+ ∂τ xy∂y

+ ∂σz

∂z + bz wdV = 0

to

(7.37)

Com o conhecimento das relações constitutivas é possível descrevero campo de tensões em função dos deslocamentos u(x,y,z ), v(x,y,z ) ew(x,y,z ), porém é mais fácil manipular algebricamente as equações escri-tas em função do tensor.

Para obter a forma fraca deve-se aplicar o teorema de Green-Gauss na

eq. (7.36). O resultado final da forma fraca é: A

(σx()x + σy()y + σz()z + τ xy(γ )xy + τ yz(γ )yz + τ xz(γ )xz) dV = S

(q xu + q yv + q zw) dS +

V

(bxu + byv + bzw) dV (7.38)

sendo S a parte da face onde as forças são aplicadas e q x = σxnx +τ xyny +τ xznz,etc. Além disto, tem-se que:

193



∂ u∂x

= ()x; ∂ v∂y

= ()y∂ w∂z

= ()z∂ u∂y

+ ∂ v∂x

= (γ )xy∂ v∂z

+ ∂ w∂y

= (γ )yz∂ u∂z

+ ∂ w∂x

= (γ )zx

As equações de elementos finitos podem ser derivadas usando funções deinterpolação nos três deslocamentos:

u(x,y,z ) = N 1u1 + N 2u2 + · · · (7.39)

v(x,y,z ) = N 1v1 + N 2v2 + · · · (7.40)

w(x,y,z ) = N 1w1 + N 2w2 + · · · (7.41)sendo u1, v1, w1, u2,· · · os os graus de liberdade nos nós do elemento

sólido e N i(x,y,z ) as funções de formas do elemento usado. Estes termospodem ser descritos de forma matricial:

u(x,y,z ) =

uvw

=

N 1 0 0 N 2 0 · · ·0 N 1 0 0 N 2 · · ·0 0 N 1 0 0 · · ·

u1

v1w1

u2

...

= NT d

(7.42)Da solução assumida o vetor deformação para o elemento pode ser calcu-

lado por:

=

xyz

γ xy

γ yzγ zx

=

∂N 1∂x

0 0 ∂N 2∂x

0 0 · · ·0 ∂N 1

∂y0 0 ∂N 2

∂y0 · · ·

0 0 ∂N 1∂z

0 0 ∂N 2∂z

· · ·∂N 1∂y

∂N 1∂x

0 ∂N 2∂y

∂N 2∂x

0 · · ·0

∂N 1

∂z

∂N 1

∂y 0∂N 2

∂z

∂N 2

∂y · · ·∂N 1∂z

0 ∂N 1∂x

∂N 2∂z

0 ∂N 2∂x

· · ·

u1

v1w1

u2

...

= BT d

(7.43)Com as relações constitutivas a relação tensão-deformação é escrita como:

σ = C = CBT d (7.44)

A equação do elemento pode ser escrita por um problema algébrico como:

kd = rq + rb (7.45)

194



As forças concentradas nos nós podem ser aplicadas diretamente nos nós

e escritas na forma global no final. A matriz de rigidez k do elemento éescrita como:

k =

V

BCBT dV (7.46)

já o vetor de carga rq é:

rq = S Nsq x

q yq z ds =

S

N 1sq xdS

S N 1sq ydS

S N 1sq zdS S

N 2sq xdS ...

(7.47)

sendo Ns a função de forma específica na face S onde q atua.

7.2.1 Problema no estado plano de tensões

7.2.2 Elemento triangular com três nós

195



Referências Bibliográficas

[1] K. J. Bathe. Finite Element Procedures in Engineering Analysis . PrenticeHall, 1982.

[2] M. A. Bhatti. Fundamental Finite Element Analysis and Applications with Mathematica and Matlab Computations . John Wiley & Sons, 2005.

Documents

Analise-Estruturas1 (1)