UNIVERSIDADE FEDERAL DE UBERLÂNDIA
FACULDADE DE ENGENHARIA CIVIL
PROGRAMA DE PÓS-GRADUAÇÃO EM
ENGENHARIA CIVIL
MÉTODO DOS ELEMENTOS FINITOS
Prof. Francisco Antonio Romero Gesualdo
" [...] embora o método dos elementos finitos possa fazer um
bom engenheiro melhor, ele pode fazer um engenheiro fraco
mais perigoso." Cook (1989)
- 2010 –
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
SUMÁRIO
1 Histórico sobre o uso do MEF................................................................................................. 3 2 Considerações gerais ............................................................................................................... 3 3 GENERALIZAÇÃO DO MÉTODO DOS ELEMENTOS FINITOS .................................... 4 4 FUNDAMENTOS DA TEORIA DA ELASTICIDADE........................................................ 6 4.1 Relações tensões × deformações para elementos em estado plano de tensões ................ 8 4.1.1 Tensões principais nos elementos infinitesimais ...................................................... 9 4.1.2 Procedimento para determinação do módulo de elasticidade Ex, Ey e coeficiente de
Poisson para material ortotrópico ..................................................................................... 10 4.1.3 Determinação do módulo de elasticidade transversal (G) ....................................... 11 5 Manipulação de matrizes ....................................................................................................... 14 5.1 Autovalores e Autovetores de uma Matriz ..................................................................... 14 5.2 Matriz Definida Positiva ................................................................................................ 14 6 Matriz de rigidez ................................................................................................................... 15 6.1 Generalidades ................................................................................................................. 15 6.2 Exemplo de obtenção da matriz de rigidez .................................................................... 17 6.3 Agrupamentos da matriz de rigidez ............................................................................... 21 6.3.1 Forma geral ............................................................................................................. 21 6.3.2 Recalques de apoio .................................................................................................. 21 6.3.3 Reações de apoio ..................................................................................................... 22 7 Função aproximadora ............................................................................................................ 22 8 Funções de interpolação e funções de forma ........................................................................ 24 9 Formulação geral usando o Princípio dos Trabalhos Virtuais .............................................. 25 10 Formulação usando o Princípio dos Trabalhos Virtuais para o caso de barras de treliça
plana ......................................................................................................................................... 26 11 Exemplo numérico .............................................................................................................. 28 12 Formulação usando o Princípio dos Trabalhos Virtuais para o caso de barras de pórtico
plano ......................................................................................................................................... 28 13 Formulação pelo Método dos Elementos Finitos ................................................................ 33 13.1 Generalidades ............................................................................................................... 33 13.2 Método direto ............................................................................................................... 33 13.3 Método da energia ........................................................................................................ 36
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
APÊNDICE A – Componentes de tensões para o caso geral ............................................... 39 APÊNDICE B – Relações entre tensões e deformações ...................................................... 40 APÊNDICE C – Deformações na flexão ............................................................................. 41 APÊNDICE D – Elástica para viga bi-apoiada com força uniformemente distribuída ....... 46 APÊNDICE E – Regras para diferenciais ............................................................................ 48 APÊNDICE F – Um problema específico............................................................................ 49 APÊNDICE G – Quadratura de Gauss ................................................................................. 50 APÊNDICE H – Galerkin .................................................................................................... 51 APÊNDICE I – Uso do ANSYS .......................................................................................... 52 APÊNDICE J – Matriz de Rigidez do elemento no Sistema de Coordenadas Global [Ke] . 53 Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
3
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
1 HISTÓRICO SOBRE O USO DO MEF
1943: Courant (matemático) apresentou solução polinomial de problema de torção, seguindo
procedimentos considerados como MEF
1950: indústria aeronáutica inicia o uso do MEF para avaliar asas de aviões (Boeing)
1960: recebe o nome de MEF atribuído por Clough
1970: começam a surgir os primeiros softwares: ANSYS, NASTRAN, ASKA etc.
1980: disseminação do método
Elementos mais simples têm formulação mais simples, caso de triângulos com 6 graus de
liberdade (gl) – 3 nós -, ou retangulares com 8 gl – 4 nós. Contudo, modelar sistemas de
contorno complexo pode exigir a combinação de diferentes formas de elementos.
Usar elementos mais simples exige melhor discretização, maior refinamento da malha, maior
capacidade e tempo computacional. Pode-se então utilizar elementos mais precisos (de alta
ordem) que são mais complexos do ponto de vista de sua formulação matemática.
2 CONSIDERAÇÕES GERAIS
Para o caso de modelos bidimensionais, basicamente, empregam-se elementos do tipo
triangular e quadrangular. Estes devem ter no mínimo 3 e 4 nós, respectivamente, definidos
para cada um dos seus vértices. Também é possível o uso de elementos com nós adicionais ao
longo de suas arestas, tendo-se, então, os chamados elementos de alta ordem, mais complexos
e que permitem representações polinomiais quadráticas.
nó
nó
a) comportamento linear
b) alta ordem – comportamento não linear
Figura 1 – Elementos e quantidade de nós
Elementos com forma mais regular são sempre preferíveis. Para o caso bidimensional, a
forma quadrangular é a recomendada. Para modelos tridimensionais o hexágono deve ser o
tipo desejado. No entanto, isto nem sempre é possível, especialmente nas regiões de
transições de malhas, onde pode ser exigido o uso de elementos com formas diferenciadas
ajustáveis aos contornos, como ilustrado na Figura 2.
elemento
triangular
Y
Z
X
Figura 2 – Malha com elementos de diferentes formatos
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
4
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
Como o método é numérico, ou seja, o domínio é subdividido em partes não contínuas
(discretizadas), então, quanto maior for o número de nós e, conseqüentemente, o número de
elementos, haverá uma tendência de se obter maior precisão dos resultados. O usuário não
deverá se preocupar somente com este aspecto, pois isto pode gerar um custo computacional
relevante em termos de tempo de processamento. Outros pontos devem ser observados na
geração do modelo. Por exemplo, a escolha da malha. É uma importante etapa para análise
pelo MEF. Devem ser procurados elementos com proporcionalidades de dimensões
(regulares) evitando elementos distorcidos dos tipos mostrados na Figura 3.
b
a
β
α
a >>> b
α >>> β
quadrilátero se aproxima de um
triângulo
nós não centralizados nos lados
lados curvos
h2
h1
h1 >>> h2
Figura 3 – Tipos de elementos que devem ser evitados
O modelo completo corresponde à associação dos diversos elementos conectados pelos nós.
No processo de cálculo, apenas são garantidos que os deslocamentos associados aos graus de
liberdade dos nós de conectividade são iguais. Assim, ao longo dos lados dos elementos as
deformações são específicas para cada elemento e dependem do seu comportamento interno.
Figura 4 – Conexão entre os elementos – compatibilidade de deslocamentos apenas nos nós
3 GENERALIZAÇÃO DO MÉTODO DOS ELEMENTOS FINITOS
Fenômeno é representado por uma função matemática (equação diferencial) válida para um
certo domínio. É tratada de forma aproximada e transformada para o conjunto de subdomínios
chamados de elementos finitos.
•
subdomínios são chamados de elementos
•
pontos de contato entre elementos são chamados de nós
Se "u" representa a incógnita do problema, pode-se escrever:
m
u(x) = ∑ u (e)
e =1
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
5
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
onde u(e) são valores de "u" dentro de cada subdomínio (elemento). Como a solução é
numérica, então, u e u(e) são aproximados.
Através de procedimentos diretos, ou por resíduos ponderados, ou por método variacional
(balanço de energia) encontra-se uma expressão do tipo:
[K] {u} = {F}
Aplicando-se as condições de contorno, encontram-se valores de "u" (aproximados) para
todos os nós que pertencem aos subdomínios. Com os valores "u" conhecidos para cada nó
dos elementos determinam-se as informações internas ao elemento.
A incógnita "u" pode ser temperatura, velocidade, potencial magnético, deslocamento, fluxo
etc.
Vantagem: o problema global (mais complexo) é transformado em um problema local. Assim,
a matriz K pode ser obtida mais facilmente por acoplamento das matrizes individuais de cada
elemento.
Precisão: depende da forma, quantidade dos elementos finitos e tipo da função de
interpolação. A experiência do usuário é importante para obter bons resultados. Regiões com
variações bruscas do fenômeno exigem maior refinamento da malha a ser usada. A eficiência
da distribuição da malha (discretização) produzirá resposta eficiente.
Figura 5 – Problemas de divisão em malhas
Procedimentos:
- adota-se uma função aproximadora "u" (global) do tipo:
u(x) = u1 ϕ1 + u2 ϕ2 + u3 ϕ3 + ....... + um ϕn
onde ϕi são funções de forma global, linearmente independentes, satisfazendo as condições de
contorno.
Para um domínio considerado com 3 elementos e 4 nós, pode-se representar graficamente esta
função na forma indicada na Figura 6.
Observar a superposição de efeitos.
As funções φi, chamadas de função de forma global, são normalizadas. Podem ser lineares,
quadráticas, cúbicas etc.
Problema real é aproximado e subdividido em quantidade equivalente ao número de nós.
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
6
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
(a) = (b) : (global) equivalente a: u1 (c) + u2 (d) + u3 (g) + u4 (f) : (local)
Figura 6 – Representação gráfica das funções de forma
4 FUNDAMENTOS DA TEORIA DA ELASTICIDADE
Considerando que a maioria das formulações envolve os fundamentos da teoria da
elasticidade, seja pelas relações entre tensões e deformações ou pelas distribuições de tensões,
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
7
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
então será feita uma rápida abordagem sobre as informações básicas mencionadas
anteriormente.
Inicialmente, é importante entender como os elementos se deformam e as relações
diferenciais envolvendo estas deformações. As tensões atuantes sobre um elemento
infinitesimal tridimensional são dadas na Figura 7. Este é o caso mais geral de solicitação,
onde em cada face do elemento atuam duas tensões paralelas à face (tensões de cisalhamento)
e um tensão normal (perpendicular) ao plano.
Y
σyy
τyx
τyz
τzy
τxy
τzx
σzz
σxx
τxz
X
Z
Figura 7 – Tensões sobre um elemento infinitesimal tridimensional
Em cada direção considerada ocorrerão deformações produzidas pelas respectivas tensões,
dadas por εxx, εyy, εzz, γxy, γyz e γxz. Considerando o caso de material isótropo (mesma
propriedade em todas as direções, ou seja, E (módulo de elasticidade) e ν (coeficiente de
Poisson) iguais (constantes) em qualquer direção, tem-se a Equação 1, que relaciona tensões
com deformações, para o caso tridimensional.
υ
υ
⎡1 − ν
⎢ υ 1 −ν
⎧σ xx ⎫
υ
⎢
⎪σ ⎪
⎢ υ
υ 1 −ν
⎪ yy ⎪
⎢
⎪⎪σ zz ⎪⎪
E
0
0
⎢ 0
⎨ ⎬=
τ
(
)(
)
1
1
2
ν
ν
+
−
⎢
xy
⎪ ⎪
⎢ 0
0
0
⎪τ yz ⎪
⎢
⎪ ⎪
⎪⎩τ xz ⎪⎭
⎢
0
0
⎢⎣ 0
0
0
0
1 - 2ν
2
0
0
0
0
0
1 - 2ν
2
0
0
⎤
⎥ ⎧ε xx ⎫
⎥⎪ ⎪
⎥ ⎪ε yy ⎪
⎥
0 ⎥ ⎪⎪⎨ε zz ⎪⎪⎬
⎥ ⎪γ xy ⎪
0 ⎥ ⎪γ yz ⎪
⎥⎪ ⎪
1 - 2ν ⎥ ⎪⎩γ xz ⎪⎭
2 ⎥⎦
0
0
0
1
A seguir serão feitas algumas considerações para casos mais específicos de distribuição de
tensões. Serão tratados os casos planos (não tridimensionais), que podem ser considerados em
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
8
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
dois casos: estado plano de tensões e estado plano de deformações, conforme será descrito a
seguir.
4.1 Relações tensões × deformações para elementos em estado plano de tensões
A Figura 8 ilustra um elemento sujeito ao estado plano de tensões. Este elemento é
considerado de espessura pequena em relação às demais dimensões e, portanto, é um
elemento bidimensional com tensões (σz) nulas na direção perpendicular ao seu plano. É o
caso típico de elementos associados a vigas. Nestes elementos existem deformações
transversais, porém não existem tensões uma vez que não há nenhum impedimento nesta
direção, permitindo a livre deformação.
Para o desenvolvimento das equações relacionadas com tensões, foi considerado que as
tensões de tração são positivas e as de compressão são negativas.
σyy
τyx
τxy
σxx
σxx
τxy
τyx
σyy
Figura 8 – Elemento submetido ao estado plano de tensões
Se o material é considerado com propriedades diferentes nas duas direções, então isto
significa que para cada nível de carregamento é preciso caracterizar cinco diferentes
constantes: E x , E y ,ν xy ,ν yx e G .
Destas constantes, Ex e Ey são imediatas, pois são obtidas diretamente pelas equações que
relacionam tensão e deformação do concreto, ou seja, para cada instante do carregamento
têm-se as deformações calculadas, através das quais se obtêm as tensões e, conseqüentemente,
os módulos de elasticidade secantes Ex e Ey pela relação σ/ε. Para os coeficientes de Poisson é
necessário estabelecer relações que permitam seu cálculo a partir das solicitações, como será
visto na Seção 4.1.2. Para a determinação do módulo de deformação transversal será feito um
desenvolvimento específico para relacioná-lo com os parâmetros associados, conforme Seção
4.1.3.
Por definição, os coeficientes de Poisson aqui empregados correspondem à relação entre
deformações na direção passiva (perpendicular à direção da ação da força) e a deformação
ativa (direção da força), com um sinal negativo, pois tração na direção ativa implica em
encurtamento na direção passiva – Equação 2. Importante é notar que as deformações
indicadas nestas equações correspondem exclusivamente àquelas associadas à aplicação
isolada da tensão na respectiva direção x ou y.
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
9
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
ν xy = −
εy
εx
e
ν yx = −
εx
εy
2
4.1.1 Tensões principais nos elementos infinitesimais
Nos elementos infinitesimais de uma estrutura sob ação de solicitações surgem tensões
normais de tração e/ou compressão (σx e σy), bem como tensões tangenciais (τxy), conforme
ilustra a Figura 9. Estas tensões podem ser representadas matricialmente como indicado na
Equação 3.
⎡σ x
[σ ] = ⎢
⎣τ xy
τ xy ⎤
σ y ⎥⎦
3
σy
τxy
σx
dy
σx
dx
τxy
σy
Figura 9 – Elemento infinitesimal sujeito a tensões normais e cisalhantes
Podem ser decompostas em direções preferenciais chamadas de principais, onde as tensões
normais ortogonais têm um valor máximo e mínimo, como mostrado na Figura 10. Estas
tensões ocorrem para o elemento infinitesimal girado de um ângulo α em relação à posição
inicial, tal que a tensão tangencial seja nula.
σ1
σ2
α
σ2
σ1
Figura 10 – Tensões principais em elemento infinitesimal
Para obter as tensões equivalentes sobre o elemento infinitesimal girado pelo ângulo α devese aplicar a Equação 4, onde αij é uma matriz transformação dada pela Equação 5. Isto
corresponde à rotação do sistema de coordenadas retangulares.
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
10
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
[σ ] = [α ij ]
T
⎡σ 1α ⎤
⋅ [σ ]⋅ α ij = ⎢⎢σ 2α ⎥⎥
⎢⎣ τ α ⎥⎦
α
[α ] = ⎡⎢−cos
senα
ij
⎣
[ ]
4
senα ⎤
cos α ⎥⎦
5
Desenvolvendo a Equação 4 e aplicando as relações trigonométricas apropriadas, obtêm-se as
Equações 6, 7 e 8 para as tensões [σ ].
⎛σ x +σ y ⎞ ⎛σ x +σ y ⎞
⎟⎟ + ⎜⎜
⎟⎟ ⋅ cos 2α + τ xy ⋅ sen2α
⎝ 2 ⎠ ⎝ 2 ⎠
6
⎛σ x +σ y ⎞ ⎛σ x −σ y ⎞
⎟⎟ − ⎜⎜
⎟⎟ ⋅ cos 2α − τ xy ⋅ sen2α
⎝ 2 ⎠ ⎝ 2 ⎠
7
⎛σ x −σ y ⎞
⎟⎟ ⋅ sen2α + τ xy ⋅ cos 2α
2
⎝
⎠
8
σ 1α = ⎜⎜
σ 2α = ⎜⎜
τ α = −⎜⎜
Derivando as Equações 6, 7 e 8 parcialmente em relação a α e igualando-as a zero, obtêm-se
os valores máximos e mínimos referentes às tensões principais, respectivamente, σ1 e σ2.
2
⎛σ +σ y ⎞
⎟⎟ +
σ 1 = ⎜⎜ x
⎝ 2 ⎠
⎛σ x −σ y ⎞
⎜⎜
⎟⎟ + τ xy 2
2
⎝
⎠
⎛σ +σ y ⎞
⎟⎟ −
σ 2 = ⎜⎜ x
⎝ 2 ⎠
⎛σ x −σ y ⎞
⎜⎜
⎟⎟ + τ xy 2
⎝ 2 ⎠
9
2
10
Estas tensões ortogonais ocorrem para a rotação α que vale:
1⎡
2 ⎢⎣
⎞⎤
⎟⎥
⎟⎥
σ
σ
−
x
y
⎝
⎠⎦
⎛ 2τ xy
α = ⎢ tan −1 ⎜⎜
11
4.1.2 Procedimento para determinação do módulo de elasticidade Ex, Ey e coeficiente
de Poisson para material ortotrópico
Serão tratados separadamente os efeitos das deformações longitudinais (εx e εy) e da
deformação transversal (γxy). Assim, no desenvolvimento a seguir serão considerados apenas
os efeitos de εx e εy. Com esta simplificação, a Equação 12 representa as relações constitutivas
para um material ortotrópico para o estado plano de tensões, aqui adotado como referência
para o material em estudo. Estas equações são facilmente deduzidas, como é apresentado em
detalhes em Carroll (1999, p.164), Lekhnitskii (1981), dentre outros.
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
11
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
εx =
ν yx
1
σx −
σy
Ex
Ey
ν xy
1
εy = −
σx +
σy
Ex
Ey
ou
⎡ 1
⎧ε x ⎫⎢
⎪ ⎪⎢ E x
⎨ ⎬⎢
⎪ ⎪⎢ ν xy
⎩ε y ⎭⎢− E
x
⎣
ν yx ⎤
⎥ ⎧σ x ⎫
E y ⎥⎪ ⎪
⎥⎨ ⎬
1 ⎥⎪ ⎪
σy
E y ⎥⎦ ⎩ ⎭
−
ou
{ε } = [D]−1 {σ }
12
Como a matriz das constantes elásticas [D]-1 é simétrica, os termos da diagonal ascendente
podem ser trocados. Assim a Equação 12 pode ser reescrita na forma indicada na Equação 13.
εx =
ν
1
1
(σ x −ν xyσ y )
σ x − xy σ y =
Ex
Ex
Ex
ν yx
1
1
(−ν yxσ x + σ y )
εy = − σx + σy =
Ey
Ey
Ey
13
Desta forma, o módulo de elasticidade na direção x pode ser escrito usando a primeira
equação de 13, resultando na Equação 14.
Ex =
σ x − ν xyσ y
εx
14
O módulo de elasticidade Ey é então obtido usando a Equação 12 vinculada a εy. O uso desta
equação mantém a dependência de νxy, resultando na Equação 15.
Ey =
σy
ν xyσ x
εy +
15
Ex
Tendo-se os valores de Ex, Ey e νxy, torna-se possível determinar νyx tendo em vista a simetria
E
ν
ν
da matriz [D]-1, ou seja, yx = xy . Portanto: ν yx = y ν xy .
Ex
E y Ex
Deve ser lembrado que a matriz [D] completa para materiais ortotrópicos (caso em questão) é
dada pela Equação 16.
⎡ Ex
⎢1 −ν ν
xy yx
⎢
⎢
ν E
[D ] = ⎢⎢ xy y
1 −ν xyν yx
⎢
⎢
0
⎢
⎢⎣
ν yx E x
1 −ν xyν yx
Ey
1 −ν xyν yx
0
⎤
0⎥
⎥
⎥
0⎥
⎥
⎥
⎥
G⎥
⎥⎦
16
4.1.3 Determinação do módulo de elasticidade transversal (G)
Como hipótese básica, será considerada a validade da superposição de efeitos, ou seja, a ação
conjunta das tensões σx, σy e τxy equivale à soma dos efeitos de cada uma agindo
isoladamente. Portanto, nesta seção será considerado o caso de um elemento infinitesimal
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
12
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
sujeito somente ao efeito das tensões de cisalhamento τxy, admitindo que esta tensão é a única
responsável pela geração da distorção angular (γ). A Figura 11 esquematiza a situação
considerada.
dx
τ
τ
y’
γ’
d
dy
γ = 2γ’
γ’
τ
Δd1
γ’
τ
Δd2
x’
Figura 11-Elemento finito sujeito ao cisalhamento puro
As deformações em um instante de carregamento podem ser obtidas pelas Equações 17 e 18,
onde Δd1 e Δd2 são os valores absolutos dos deslocamentos nas diagonais caracterizadas pelos
eixos x’ e y’, respectivamente, na forma de alongamento (positivo) e encurtamento
(negativo).
Δd1 = ε x ' ⋅ d
17
Δd 2 = ε y ' ⋅ d
18
Geometricamente, a distorção angular γ vale 2 γ’, portanto pode ser obtida pela Equação 19.
γ =
Δd1 Δd 2
+
d
d
19
Substituindo a Equação 17 e 18 em 19 e lembrando que ε x ' é positivo e ε y ' é negativo,
obtêm-se 20.
γ = ε x' − ε y'
20
As equações de ε x ' e ε y ' são estabelecidas conforme a Equação 13.
ε x' =
(
1
σ x ' −ν xy' ⋅ σ y '
Ex '
(
)
1
ε y' =
σ y ' −ν yx' ⋅ σ x '
Ey'
)
21
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
13
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
Sabendo-se que σx’ = -σy’ = τ (por equilíbrio estático de forças do elemento da Figura 11), e
substituindo as Equações 21 em 20, obtêm-se através da 22 a equação que relaciona distorção
angular com a tensão cisalhante.
'
'
γ 1 ⎛⎜ 1 +ν xy 1 +ν yx ⎞⎟
= =
+
E y ' ⎟⎠
τ G ⎜⎝ E x '
22
Portanto, com o inverso do valor dado na Equação 22, tem-se o módulo de deformação
'
, determinados para cada instante do
transversal (G) em função dos valores de E x ' , E y ' , ν xy' e ν yx
carregamento. Este mesmo valor é mostrado por Carroll (1999, p.168), cuja dedução é feita
usando conceitos da energia de deformação.
Assim, a equação completa que transforma tensão em deformação para tensões planas, em
função das constantes básicas, é dada pela Equação 23.
⎡ 1
⎧εx ⎫ ⎢
⎪ ⎪ ⎢ Ex
⎪ ⎪ ⎢
⎪ ⎪ ⎢ ν xy
⎨ ε y ⎬ = ⎢−
⎪ ⎪ ⎢ Ex
⎪ ⎪ ⎢
⎪γ ⎪ ⎢ 0
⎩ xy ⎭ ⎢
⎣
−
ν yx
Ey
1
Ey
0
⎤
⎥ ⎧σ x ⎫
⎥⎪ ⎪
⎥⎪ ⎪
⎥⎪ ⎪
0
⎥ ⎨σ y ⎬
⎥⎪ ⎪
⎪ ⎪
'
⎞⎥ ⎪ ⎪
1 + ν yx
⎟⎥ τ xy
+
E y ' ⎟⎠⎥⎦ ⎩ ⎭
0
⎛ 1 + ν xy'
⎜
⎜ E
⎝ x'
23
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
14
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
5 MANIPULAÇÃO DE MATRIZES
As formulações computacionais, em geral, estão baseadas em operações matriciais, tendo em
vista a facilidade de manipulação das informações. Assim, é indispensável conhecer
informações relacionadas às características de matrizes e suas aplicações.
5.1 Autovalores e Autovetores de uma Matriz
Os autovalores de um matriz [K] de ordem n×n são definidos como os valores de λ que
satisfazem a relação [K]{u}=λ{u} para {u}≠0. Esta expressão pode ser rearranjada pelo
posicionamento de todos os termos no primeiro membro e, também, introduzindo a matriz
identidade [I], ou seja, ([K] - λ[I]){u}=0.
A solução existe somente se a matriz ([K] - λ[I]) for singular, o que significa que det([K] λ[I])=0. O polinômio resultante desta condição é chamado de polinômio característico, cujas
raízes são os autovalores da matriz [K].
⎡7 3⎤
Exemplo: encontrar os autovalores para a matriz [K ] = ⎢
⎥.
⎣3 5⎦
3 ⎤
⎡7 − λ
Det([K]- λ[I]) = det ([K ] − λ [I ]) = ⎢
= λ2 − 12λ + 26 = 0 . O polinômio terá duas
⎥
5 − λ⎦
⎣ 3
soluções: λ1 = 9,162 e λ2 = 2,838, que são os autovalores de [K].
A todo autovalor λi tem-se um autovetor associado {ui}, tal que, [K]{ui}=λi{ui}.
5.2 Matriz Definida Positiva
Uma matriz [K] de ordem n×n é definida positiva se para qualquer vetor {u}≠0 a relação
{u}T[K]{u} > 0 é atendida. A matriz [K] é definida positiva se e somente se todos seus
autovalores são maiores que zero. Outra condição é que todos os determinantes das matrizes
da série indicada a seguir sejam positivos, incluindo o determinante da matriz completa.
K
[K11 ], ⎡⎢ 11
⎣ K 21
⎡ K11
K12 ⎤ ⎢
, K 21
K 22 ⎥⎦ ⎢
⎢⎣ K 31
K12
K 22
K 32
K13 ⎤
K 23 ⎥⎥ etc.
K 33 ⎥⎦
Duas outras condições são necessárias, mas não suficientes: Kii > 0 e K ij < K ii K jj (i≠j)
A matriz é semi-definida positiva se pelo menos um autovalor é zero e os demais são
positivos.
Exemplos:
⎡7 3⎤
No exemplo do caso anterior, a matriz [K ] = ⎢
⎥ tem os autovalores λ1 = 9,162 e λ2 =
⎣3 5⎦
2,838, ambos positivos, portanto, a matriz [K] é definida positiva.
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
15
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
⎡8 1 2 ⎤
Outro exemplo: seja a matriz [K ] = ⎢⎢1 7 5⎥⎥ . Esta apresenta os autovalores λ1 = 11,795, λ2
⎢⎣2 5 4⎥⎦
= 7,061 e λ3 = 0,144, todos positivos, ou seja, é uma matriz definida positiva.
6 MATRIZ DE RIGIDEZ
6.1 Generalidades
Como já comentado a equação básica que rege os problemas estruturais via método dos
elementos finitos obedece à equação [K] {u} = {F} do tipo indicado na Equação 24.
⎡k 11
⎢k
⎢ 21
⎢k 31
⎢
⎢ ...
⎢k n1
⎣
k 12
k 13
k 22
k 23
k 32
k 33
...
...
k n2
k n3
... k 1n ⎤ ⎧u1 ⎫ ⎧F1 ⎫
... k 2n ⎥⎥ ⎪⎪u2 ⎪⎪ ⎪⎪F2 ⎪⎪
⎪ ⎪ ⎪ ⎪
... k 3n ⎥ ⎨u3 ⎬ = ⎨F3 ⎬
⎥⎪ ⎪ ⎪ ⎪
... ... ⎥ ...
...
⎪ ⎪ ⎪ ⎪
... k 9n ⎥⎦ ⎪⎩un ⎪⎭ ⎪⎩Fn ⎪⎭
Equação 24
Isto significa que a relação entre força e deslocamento é diretamente interligada por um
coeficiente chamado de rigidez (kij). A matriz de rigidez é, então, formada por este conjunto
de coeficientes de rigidez que interconecta qualquer coordenada definida para o problema. Na
prática, é necessário impor as condições de contorno − deslocamentos impostos − para tornar
o sistema de equações invertível, ou seja, com solução. A partir deste sistema de equações é
possível determinar quaisquer deslocamentos e reações de apoio da estrutura. Na prática
torna-se interessante agrupar o sistema de equações em partes que tenham características
afins, como será visto na Seção 6.3.
Para um caso simples de uma mola de constante de rigidez k, a força aplicada corresponde ao
deslocamento da mola multiplicado pela sua constante, ou seja, F = k u (Figura 12). Isto
significa que para qualquer u pode-se obter o valor de F e vice-versa.
c
d
Matricialmente, o sistema fica representado por:
F
k = cte da mola
⎡ k11 k12 ⎤ ⎧ u1 ⎫ ⎧ F1 ⎫
⎢
⎥⎨ ⎬ = ⎨ ⎬
⎣k 21 k 22 ⎦ ⎩u 2 ⎭ ⎩F2 ⎭
Figura 12 – Mola com força aplicada
Como u1=0, então, a única incógnita restante será u2. Equivale a eliminar a 1ª. linha e 1ª.
coluna, ou seja, [k22]{u2}={F2}. No exemplo, k22=k e F2=F. Então: u2 = F/k (u2 é o
deslocamento do nó 2).
Seja agora um arranjo estrutural formado por duas molas, em seqüência, de constantes ka e kb,
solicitadas pelas forças Pa e Pb que atuam sobre os nós 2 e 3, respectivamente, como mostrado
na Figura 13.
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
16
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
c
d
e
Pa
ka
Pb
kb
Figura 13 – Sistema formado por duas molas
Tabela 1 – Matriz de rigidez para sistema formado por duas molas
Deslocamentos impostos {u}
Equilíbrio estático
ka
1
1
c
nó
d
e
2
F1 = k11 = ka
F1
0
ka
F2 = k12 = -ka
F2
⎧u1 ⎫ ⎧1⎫
⎪ ⎪ ⎪ ⎪
⎨u2 ⎬ = ⎨0⎬
⎪u ⎪ ⎪0⎪
⎩ 3⎭ ⎩ ⎭
0
3
F3 = k13 = 0
F3
ka
1
F1 = k21 = -ka
F1
1
kb
ka
2
F2=k22=ka + kb
F2
⎧u1 ⎫ ⎧0⎫
⎪ ⎪ ⎪ ⎪
⎨u2 ⎬ = ⎨1⎬
⎪u ⎪ ⎪0⎪
⎩ 3⎭ ⎩ ⎭
kb
3
F3 = k32=-kb
F3
0
1
1
kb
0
2
⎧u1 ⎫ ⎧0⎫
⎪ ⎪ ⎪ ⎪
⎨u2 ⎬ = ⎨0⎬
⎪u ⎪ ⎪1⎪
⎩ 3⎭ ⎩ ⎭
F1=k13 = 0
F1
F2 = k23=-kb
F2
kb
3
F3 = k33=kb
F3
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
17
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
O sistema de equações que irá representar este arranjo estrutural é formado por uma matriz de
rigidez 3×3, um vetor de deslocamentos e outro de forças com três elementos cada um, pois o
arranjo é formado por 3 nós com um grau de liberdade para cada nó, dado pela Equação 25.
⎡k 11 k 12
⎢k
⎢ 21 k 22
⎢⎣k 31 k 32
k 13 ⎤ ⎧u1 ⎫ ⎧F1 ⎫
⎪ ⎪ ⎪ ⎪
k 23 ⎥⎥ ⎨u2 ⎬ = ⎨F2 ⎬
k 33 ⎥⎦ ⎪⎩u3 ⎪⎭ ⎪⎩F3 ⎪⎭
Equação 25
Para obtenção da matriz de rigidez utiliza-se o artifício algébrico de impor deslocamento
unitário na coordenada correspondente à coluna da matriz de rigidez a ser determinada. A
Tabela 1 apresenta os deslocamentos impostos, as forças envolvidas no problema e o
equilíbrio dos nós que resultam na obtenção da matriz de rigidez dos arranjo estrutural em
questão. Resulta, portanto, na matriz de rigidez dada pela Equação 26.
0 ⎤ ⎧u1 ⎫ ⎧F1 ⎫
− ka
⎡ ka
⎢− ka ka + kb − kb⎥ ⎪u ⎪ = ⎪F ⎪
⎢
⎥⎨ 2 ⎬ ⎨ 2 ⎬
⎢⎣ 0
kb ⎥⎦ ⎪⎩u3 ⎪⎭ ⎪⎩F3 ⎪⎭
− kb
Equação 26
Ao impor u1 = 0, a multiplicação da 1ª. linha de [K] pelo vetor u resulta em F1 = -ka×u2. Notar
que F1 é a reação de apoio no nó 1 (ponto fixo), conhecida após a determinação de u2.
Uma vez que se conhece u1, resta apenas conhecer u2 e u3. Assim, a solução do sistema
resume-se à um sistema de equações envolvendo estes deslocamentos incógnitos. Isto
equivale a eliminar a 1ª. linha e a 1ª. coluna de [K] gerando, enfim, a Equação 27.
⎡ka + kb − kb⎤ ⎧u2 ⎫ ⎧F2 ⎫
⎨ ⎬=⎨ ⎬
⎢ − kb
kb ⎥⎦ ⎩u3 ⎭ ⎩F3 ⎭
⎣
Equação 27
Fazendo F2 = Pa e F3 = Pb, que são as forças atuantes na estrutura (Figura 13 – pág. 16) e
resolvendo este sistema de equações, tem-se:
u2 =
Pa + Pb
ka
e
u3 =
1
(Pa + Pb ) + Pb
ka
kb
Substituindo u2 na equação F1 = -ka×u2 (deduzida anteriormente), tem-se F1 = -(Pa + Pb),
como esperado. O sinal negativo significa que a reação F1 é uma força com sentido da direita
para a esquerda, pois F1 é positiva da esquerda para a direita, sentido da coordenada positiva
imposto no início da solução do problema.
Sendo a estrutura analisada bastante simples, como exercício, confirme os valores de u2 e u3
avaliando a estrutura de forma direta, sem o procedimento empregado.
6.2 Exemplo de obtenção da matriz de rigidez
Para analisar a obtenção da matriz de rigidez será usado o exemplo dado na Figura 14, sobre a
qual atuam as forças FA, FB, FC. Os nós foram numerados de 1 a 3, de forma aleatória. Mais
adiante será visto como otimizar a numeração de nós.
O desenvolvimento do exemplo em questão tem como objetivo mostrar a manipulação
numérica de informações estruturais sem, contudo, representar um procedimento padrão. É
apenas uma preparação do ambiente usando ferramentas não convencionais, mas que traz a
idéia matricial e mostra recursos importantes no trato numérico do cálculo de estruturas.
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
18
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
FA
FB
nó 1
FC
nó 2
nó 3
Figura 14 – Exemplo numérico
O tipo de estrutura (pórtico plano) apresenta três tipos de incógnitas por nó: translação
horizontal, translação vertical e rotação. Isto significa que cada nó terá três graus de liberdade
e, portanto, para a estrutura da Figura 14 serão nove graus de liberdade, conforme indicado na
Figura 15. A numeração destes nós deve seguir uma seqüência iniciada pela numeração da
coordenada horizontal, a vertical e finalmente a de rotação. Assim, haverá uma regra de
numeração seqüencial para favorecer o desenvolvimento do problema, conforme Figura 15b.
2
1
3
5
3i-1
8
4
3i-2
7
3i
nó i
9
6
(a)
(b)
Figura 15 – Graus de liberdade
A estrutura analisada, que tem nove graus de liberdade, terá uma matriz de rigidez 9×9, nove
deslocamentos incógnitos (ui) e nove possíveis solicitações externas (Fi) aplicadas nas
coordenadas indicadas (Equação 35).
⎡ k 11
⎢k
⎢ 21
⎢ k 31
⎢
⎢ ...
⎢⎣ k 91
k 12
k 13
k 22
k 23
k 32
...
k 33
...
k 92
k 93
... k 19 ⎤ ⎧ u 1 ⎫ ⎧ F1 ⎫
... k 29 ⎥⎥ ⎪⎪u 2 ⎪⎪ ⎪⎪F2 ⎪⎪
⎪ ⎪ ⎪ ⎪
... k 39 ⎥ ⎨u 3 ⎬ = ⎨F3 ⎬
⎥
... ... ⎥ ⎪ ... ⎪ ⎪ ... ⎪
⎪ ⎪ ⎪ ⎪
... k 99 ⎥⎦ ⎪⎩u 9 ⎪⎭ ⎪⎩F9 ⎪⎭
Equação 28
Um procedimento para obter os valores de kij corresponde a impor deslocamentos (vetor ui)
unitários da seguinte forma:
•
para obter a primeira coluna da matriz [K]: impor {u}T = {1 0 0 ... 0}*
•
para obter a segunda coluna da matriz [K]: impor {u}T = {0 1 0 ... 0}
•
para obter a terceira coluna da matriz [K]: impor {u}T = {0 0 1 ... 0}
*
a representação do vetor {u} na forma transposta {u}T é usada apenas para facilitar a diagramação do texto.
Lembrar que {u} é um vetor coluna vertical e {u}T é horizontal
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
19
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
•
assim sucessivamente até a nona coluna.
Observe que ao impor um deslocamento unitário e os demais nulos significa obter um vetor
de forças que corresponde aos valores da correspondente coluna da matriz [K]. Seja por
exemplo impor {u}T = {1 0 0 ... 0} (cada coluna multiplicada pelo vetor {u} resulta em ki1 =
Fi):
⎡ k 11
⎢k
⎢ 21
⎢ k 31
⎢
⎢ ...
⎢⎣ k 91
k 12
k 13
k 22
k 23
k 32
...
k 33
...
k 92
k 93
... k 19 ⎤ ⎧ 1 ⎫ ⎧ F1 ⎫
... k 29 ⎥⎥ ⎪ 0 ⎪ ⎪⎪F2 ⎪⎪
⎪⎪ ⎪⎪ ⎪ ⎪
... k 39 ⎥ ⎨ 0 ⎬ = ⎨F3 ⎬
⎥
... ... ⎥ ⎪...⎪ ⎪ ... ⎪
⎪ ⎪ ⎪ ⎪
... k 99 ⎥⎦ ⎩⎪ 0 ⎭⎪ ⎩⎪F9 ⎭⎪
→
⎧ F1 ⎫
⎧ k 11 ⎫
⎪F ⎪
⎪k ⎪
⎪⎪ 2 ⎪⎪
⎪⎪ 21 ⎪⎪
⎨k 31 ⎬ = ⎨F3 ⎬
⎪ ... ⎪
⎪ ... ⎪
⎪ ⎪
⎪ ⎪
⎪⎩F9 ⎪⎭
⎩⎪k 91 ⎭⎪
O vetor {F} representa o conjunto de forças aplicadas em cada coordenada definida para os
nós.
Para a determinação das forças Fi deve-se usar, por exemplo, o equilíbrio dos nós
correspondentes às forças Fi a serem determinadas. Como ilustração, ao impor {u}T = {1 0 0
... 0}, determinam-se as forças F1, F2 e F3 pelo equilíbrio do nó 1, como indicado na Figura
16. Nesta figura foi utilizada a representação típica de engaste perfeito como forma de
representar a ausência das translações e da rotação dos nós. Veja que o nó 1 é deslocado
horizontalmente (coordenada 1), mas são mantidos os deslocamentos verticais (translações) e
as rotações nulos. Também foi considerado que o elemento 1 (barra 1) faz a ligação entre o nó
2 e o nó 1, o elemento 2 conecta o nó 2 ao 3 e o elemento 3 interliga o nó 1 ao 3.
Na Figura 16b foi empregada a notação Pi( j ) para representar as solicitações atuantes nas
extremidades das barras que, conseqüentemente atuam sobre o nó. O índice i indica o tipo de
esforço (força axial: 1, força cortante: 2 e momento fletor: 3). O índice j indica o número da
barra à qual os esforços estão associados.
A determinação das forças atuantes nas coordenadas 1, 2 e 3 é feita por simples equilíbrio de
forças nas direções x, y e z. Esta última direção (z) significa equilíbrio de momentos. A
determinação dos esforços Pi( j ) é feita de forma simples ao considerar cada barra com os
“recalques de apoio” calculados a partir do deslocamento unitário imposto, bastando saber
que estes esforços são obtidos pela expressão indicada na Figura 17. Os deslocamentos
(recalques diferenciais) entre as extremidades das barras são obtidos pela projeção do
deslocamento unitário imposto sobre o nó 1.
Para se determinar os demais valores de Fi (i de 4 a 9) será necessário fazer o equilíbrio dos
nós 2 e 3 de forma análoga ao feito para o nó 1.
É interessante notar que a matriz de rigidez será o resultado de uma série de equilíbrio de
forças, ou seja, apesar de todas as técnicas empregadas o fundamento básico é sempre o
equilíbrio estático para a solução de estruturas.
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
20
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
1
1
2
3
a) estrutura com deslocamento u1 = 1 imposto
F2
y
F1
x
F3
P2(5)
P1(3)
P1(4)
P3(6)
P2(3)
P3(3)
b) equilíbrio de forças do nó 1
Figura 16 – Estrutura com deslocamentos impostos {u}T={1 0 0 ... 0}.
y
A
EI
x
B
u1=δ
L
MBA
MAB
VA
VB
Esforços nas extremidades da barra AB
com sentidos positivos de acordo com o
sistema de coordenadas (x, y):
Momentos fletores:
M AB = M BA =
Forças cortantes:
VA = −VB =
6EI
δ
L2
12EI
δ
L3
Figura 17 – Momentos fletores nas extremidades de barra com recalque
Ao final da montagem da matriz de rigidez será observado que para o processo ser
automatizado é mais adequado adotar uma ordem em que as contribuições à matriz de rigidez
sejam dadas pela seqüência de barras. Isto significa que as parcelas de forças envolvidas no
equilíbrio dos nós sejam referidas às contribuições das barras, como será abordado mais
adiante de forma mais detalhada, especialmente porque representa um procedimento dentro
do método dos elementos finitos.
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
21
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
6.3 Agrupamentos da matriz de rigidez
6.3.1 Forma geral
Uma matriz de rigidez é originalmente composta por todos os coeficientes de rigidez que
relacionam as coordenadas definidas para a estrutura. A solução do sistema de equações
somente é possível depois de se impor as condições de contorno relativas às restrições de
apoio. Para facilitar o desenvolvimento matemático, agrupam-se os coeficientes não
vinculados às coordenadas sem restrições − sem condições de contorno de vinculações −, e
aqueles com restrições de apoio. Isto é feito por renumeração das coordenadas. A Equação 9
representa a ação geral de uma estrutura ([K]{u}={F}), com os valores agrupados.
⎡ k1,1
⎢ ...
⎢
⎢ ...
⎢
⎢ ...
⎢ k m,1
⎢
⎢k m+1,1
⎢ ...
⎢
⎣⎢ k n,1
Ka
... ... ...
k1,m
k1,m+1
..
... ... ...
...
...
..
... ... ...
...
...
..
... ... ...
...
...
..
... ... ...
k m ,m
k m ,m+1
..
... ... ... k m+1,m
k m+1,m+1 ..
... ... ...
...
...
..
... ... ...
k n ,m
k n ,m+1
..
Ka
ua
Fa
k1,n ⎤ ⎧ u1 ⎫ ⎧ F1 ⎫
... ⎥⎥ ⎪ ... ⎪ ⎪ ... ⎪
⎪
⎪ ⎪
⎪
... ⎥ ⎪ ... ⎪ ⎪ ... ⎪
⎥⎪
⎪ ⎪
⎪
... ⎥ ⎪ ... ⎪ ⎪ ... ⎪
⎨
⎬=⎨
⎬
k m ,n ⎥ ⎪ u m ⎪ ⎪ Fm ⎪
⎥
k m+1,n ⎥ ⎪u m+1 ⎪ ⎪ Fm+1 ⎪
⎪
⎪ ⎪
⎪
... ⎥ ⎪ ... ⎪ ⎪ ... ⎪
⎥
k n ,n ⎦⎥ ⎪⎩ u n ⎪⎭ ⎪⎩ Fn ⎪⎭
Forças efetivamente aplicadas
sobres os nós
≡
Kb
Kb
ub
Fb
Reações
de apoio
29
⎡ K aa K ab ⎤ ⎧ua ⎫ ⎧ Fa ⎫
⎢ K K ⎥ ⎨u ⎬ = ⎨ F ⎬
bb ⎦ ⎩ b ⎭
⎣ ba
⎩ b⎭
valores nulos (quando existem apoios fixos)
ou valores de recalques de apoio (impostos)
6.3.2 Recalques de apoio
Deve se lembrado que o sistema útil a ser resolvido corresponde a [Ka]{ua}={Fa}, as demais
equações são redundantes. Contudo, carregamento na forma de recalque de apoio implica na
imposição de deslocamentos fora das coordenadas com graus de liberdade. Assim, esta
situação exige um tratamento especial, facilmente entendido manipulando as equações
indicadas na Equação 29. Os deslocamentos de recalques de apoio serão armazenados na
região ub, ou seja, deslocamentos armazenados como um+1 a un. Isto implica em rearranjar o
sistema de equações pela alteração do vetor {Fa}. Este passa a ser composto pelas forças
externas aplicadas diretamente sobre o nó correspondente mais o efeito do recalque
encontrado pela multiplicação da linha da matriz de rigidez (coordenadas m+1 a n) pelos
respectivos recalques impostos. Isto significa passar para o segundo membro o efeito do
recalque. Para a montagem do primeiro valor do vetor {F}, resulta em F1 – k1,m+1×um+1 – k1,
m+2×um+2 – ... ... ... – k1×un. Para cada valor de F haverá essa reestruturação do sistema de
equações resultando nas mesmas m equações, apenas alteradas pelos termos independentes
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
22
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
representados por F1, F2, ... Fm, provenientes da multiplicação de [kab]×{ub}, respectivamente.
A partir disto, basta resolver o sistema m × m.
6.3.3 Reações de apoio
A solução de uma estrutura tendo qualquer tipo de carregamento implica na solução do
sistema [Ka]{ua}={Fa}, onde {Fa} pode corresponder às forças aplicadas sobre os nós ou
também estar modificado pelas ações de recalques de apoio. De qualquer forma, a solução
indicará todos os valores de deslocamentos desconhecidos, além daqueles impostos. No
entanto, o vetor {Fb} continuará desconhecido, mas corresponde às forças atuantes nas
coordenadas dos vínculos de apoio e, portanto, são as reações de apoio. Desta forma, estas
incógnitas serão facilmente conhecidas resolvendo a equação [Kba Kbb]{u} = {Fb}. Este
sistema envolve a multiplicação da matriz de ordem m-n × n por um vetor {u} de ordem n,
resultando no vetor {Fb}. de ordem m-n, que são as reações de apoio.
7 FUNÇÃO APROXIMADORA
Os valores de momentos fletores nas extremidades das barras − muitos chamam de
engastamento perfeito − indicados na Figura 17 são provenientes de cálculos convencionais e
aparecem na literatura da área. Contudo, estes resultados podem ser facilmente encontrados se
utilizarmos um polinômio para representar os deslocamentos “u” associados ao problema. A
princípio isto significa uma aproximação, porém será notado que neste caso coincidirá com a
resposta exata.
Então, seja adotado um polinômio do tipo w = c1 + c2 x + c3 x2 + c4 x3 para representar a
elástica da barra (foi adotada a letra w para representar os deslocamentos verticais (para
baixo) de uma barra isolada do tipo ilustrado na Figura 18). Notar que o grau máximo deste
polinômio deve ser 3, pois nele existirão 4 incógnitas (c1, c2, c3 e c4), tendo em vista que à
barra somente poderão ser definidas quatro condições de contorno (translação vertical e
rotação das duas extremidades da barra), associadas às coordenadas definidas na Figura 18.
1
2
3
4
Figura 18 – Sistema de coordenadas
Lembrando-se que a rotação da barra é dada pela primeira derivada da equação da elástica,
tem-se:
w = c1 + c 2 x + c 3 x 2 + c 4 x 3
w′ =
dw
= θ = c2 + 2 ⋅ c3 x + 3 ⋅ c4 x 2
dx
Impondo as condições de contorno associadas ao caso, é possível criar um sistema de
equações envolvendo os deslocamentos (w) e as incógnitas (ci). Tem-se então:
⎧w = c1
x=0⇒⎨ 1
⎩w2 = c 2
⎧⎪w3 = c1 + c 2L + c 3L2 + c 4L3
x =L ⇒ ⎨
⎪⎩w 4 = c 2 + 2 ⋅ c 3L + 3 ⋅ c 4L2
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
23
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
As equações anteriores podem ser reescritas matricialmente na forma da Equação 30:
⎧w1 ⎫ ⎡1
⎪ ⎪ ⎢
⎪w2 ⎪ ⎢0
⎨ ⎬=
⎪w3 ⎪ ⎢⎢1
⎪⎩w 4 ⎪⎭ ⎣0
⎤ ⎧ c1 ⎫
⎥⎪ ⎪
⎥ ⎪⎨c 2 ⎪⎬
⎥ ⎪c 3 ⎪
2 ⎥⎪
1 2L 3L ⎦ ⎩c 4 ⎪⎭
0 0
1 0
L L2
0
0
L3
Equação 30
A Equação 30 é genérica, ou seja, vale para qualquer {w} imposto. Especificamente para o
caso em estudo {w}T={1 0 0 0}, ver convenção de sinais observando a Figura 33b na pág. 43
(Apêndice A), o sistema resulta em:
⎧1⎫ ⎡1
⎪0⎪ ⎢0
⎪ ⎪ ⎢
⎨ ⎬=
⎪0⎪ ⎢⎢1
⎪⎩0⎪⎭ ⎣0
⎤ ⎧c1 ⎫
⎥ ⎪c ⎪
⎥ ⎪⎨ 2 ⎪⎬
L L2 L3 ⎥ ⎪c3 ⎪
⎥
1 2L 3L2 ⎦ ⎪⎩c 4 ⎪⎭
0
1
0
0
0
0
Resolvendo o sistema de equações, obtem-se:
c1 = 1
c2 = 0
c3 = -3/L2
c4 = 2/L3
Substituindo estes valores em w = c1 + c2 x + c3 x2 + c4 x3 tem-se a equação da elástica para o
caso avaliado:
w =1 −
3
2
L
x2 +
2
3
L
x3
Equação 31
A partir da equação da elástica pode-se obter a equação do momento fletor ao longo do
comprimento da barra empregando-se a expressão clássica da flexão (ver Apêndice A):
EI
d2 w
dx 2
= −M
Equação 32
Esta expressão com o sinal negativo associado a M, significa que o momento fletor traciona a
borda superior, uma vez que os valores de w são positivos para cima. O Apêndice A traz as
justificativas detalhadas sobre esta convenção de sinais.
Aplicando a Equação 32 na Equação 31 resulta:
⎛ 6 12 ⎞
M = ⎜ 2 − 3 x ⎟ EI
L
⎝L
⎠
Como EI
12 EI
d 3w
= ±V , então: V = ± 3
3
L
dx
V
V
Portanto:
x = 0 ⇒ M AB = +
x = L ⇒ MBA = −
6EI
L2
6EI
L2
Notar que os valores de MAB e MBA aqui indicados coincidem com os valores apresentados na
Figura 17, pois como dito anteriormente, pela convenção usada para o emprego da Equação
32, estes momentos fletores tracionam as fibras superiores. Como a convenção usada para os
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
24
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
momentos fletores nas extremidades da barra é conforme apresentado na Figura 18, então se
mantém o sinal de MAB e troca-se o de MBA, portanto resultam momentos fletores totalmente
iguais, bem como para as forças cortantes.
8 FUNÇÕES DE INTERPOLAÇÃO E FUNÇÕES DE FORMA
Em formulações usando o MEF utilizam-se funções aproximadoras para estabelecer o
comportamento do elemento. As incógnitas sempre são os deslocamentos nodais
representados pelas translações e rotações.
Para introduzir a idéia das funções de interpolação e funções de forma, será usado o caso de
uma barra simples que absorve somente forças axiais (caso de barras de treliças). Portanto, a
barra fica representada por apenas duas coordenadas paralelas ao eixo da barra, posicionadas
nas extremidades da barra, conforme ilustra a Figura 19. Note que as coordenadas estão
representadas ligeiramente acima do eixo da barra apenas por questão visual, pois caso
contrário a coordenada 1 ficaria invisível ao estar coincidindo com o eixo da barra.
1
1
E, A
x
2
2
L
Figura 19 – Barra de treliça
Adotando que os deslocamentos ao longo do elemento sejam dados por uma função do tipo: u
= a + b x, é possível observar que para este caso simples esta função representará de forma
exata os deslocamentos ao longo de L.
O objetivo aqui é encontrar uma função u escrita em termos de u1 e u2 que, respectivamente,
são os deslocamentos nas extremidades 1 e 2. Estes representam as condições de contorno do
problema. Assim: u = ϕ1 u1 + ϕ2 u2. As parcelas ϕ1 e ϕ2 são chamadas de funções de forma.
Aplicando as condições de contorno tem-se:
a) para x = 0 ⇒ u = u1 ∴ u = a + b⋅0 = u1 ⇒ a = u1
b) para x = L ⇒ u = u2 ∴ u = a + b⋅L = u2 ⇒ b = (u2 - u1)/L
Ou seja:
⎛u −u ⎞
u ( x) = a + b ⋅ x = u1 + ⎜ 2 1 ⎟ x ⇒
⎝ L ⎠
x
⎛L−x⎞
u =⎜
⎟ ⋅ u1 + ⋅ u2
L
L4
⎝1
{
42
3⎠
ϕ
ϕ1
2
A equação dos deslocamentos pode, então, ser escrita de forma matricial:
⎧u ⎫
u ( x) = [ϕ1 ϕ 2 ]⋅ ⎨ 1 ⎬
⎩u 2 ⎭
Graficamente as funções ϕ1 e ϕ2 podem ser representadas conforme ilustrado na Figura 20.
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
25
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
1
1
0
l−x
ϕ1 =
l
0
l
x
ϕ2 =
l
l
Figura 20 – Funções de forma ϕ1 e ϕ2 (lineares)
Pode-se, portanto estabelecer as seguintes propriedades para ϕ1 e ϕ2:
•
têm sempre a mesma característica da função de interpolação: são sempre
polinômios do mesmo tipo;
•
uma das extremidades tem valor nulo e a outra vale 1;
•
em qualquer ponto dentro de intervalo a soma das duas funções ϕ1 + ϕ2 resulta em
valor igual a 1;
•
a soma das derivadas de ϕ1 + ϕ2 em relação a “x” vale zero.
ue
ϕ2u2
ϕ1u1
u2
u1
1
1
0
l
Figura 21 – Deslocamentos dentro de um elemento associados às funções de forma para o caso linear
9 FORMULAÇÃO GERAL USANDO O PRINCÍPIO DOS TRABALHOS VIRTUAIS
A formulação apresentada é similar ao caso da mínima energia, embora mais simples em
termos de entendimento. A formulação do princípio é baseada na afirmação de que “para um
sistema em equilíbrio ao qual se impõem deformações e deslocamentos compatíveis com suas
vinculações, pode-se afirmar que o trabalho virtual produzido pelas forças externas nos
deslocamentos é igual ao trabalho interno das tensões nas deformações internas”.
Resumidamente pode ser afirmado que:
Trabalho virtual externo = trabalho virtual interno (energia de deformação interna virtual)
Mais especificamente, para um sistema com n forças aplicadas:
n
∑ (forças reais) × (deslocamentos virtuais) = ∫ (tensões reais) × (deformações virtuais)dV
i =1
V
Nesta aplicação é importante lembrar que os deslocamentos e as deformações são impostas. A
força (ou a tensão) realiza trabalho em deslocamentos (ou deformações) pré-existentes, ou
seja, solicitação e a mudança de forma não ocorrem simultaneamente e, portanto, o trabalho
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
26
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
realizado corresponde à solicitação multiplicada pela correspondente mudança de forma, não
existindo o termo ½ na expressão do trabalho.
10 FORMULAÇÃO
USANDO O PRINCÍPIO DOS
BARRAS DE TRELIÇA PLANA
TRABALHOS VIRTUAIS
PARA O CASO DE
Adotando um polinômio linear para representar os deslocamentos ao longo da barra (Figura
22), tem-se a Equação (33).
2
1
nó j
nó k
Figura 22
u(x) = a + b x
(33)
Impondo as condições de contorno:
⎧ u1 = a + b ⋅ x j
⎨
⎩u 2 = a + b ⋅ xk
Se xj = 0 e xk = L ⇒ a = u1 e b =
Portanto: u ( x) = u1 +
u 2 − u1
L
u 2 − u1
x
x
⋅ x = (1 − ) ⋅ u1 + ⋅u 2
L
L
L
Matricialmente:
x
⎡
u ( x) = ⎢(1 − )
L
⎣
x ⎤ ⎧ u1 ⎫
e
⎨ ⎬ = [ϕ ] {u }
⎥
L ⎦ ⎩u 2 ⎭
matriz das funções de forma
[ B]
}
du
d [ϕ ] e ⎡ 1
Sabendo que ε =
, então: ε =
{u } = ⎢−
dx
dx
⎣ L
1 ⎤ ⎧u1 ⎫
e
⎨ ⎬ = [ B]{u }
⎥
L ⎦ ⎩u 2 ⎭
Usando a relação tensão × deformação dada pela equação de Hooke:
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
27
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
σ = E ⋅ ε = E ⋅ [B ] ⋅ {u e }
(34)
Por se tratar de uma análise para um caso específico de elemento, o sistema global coincide
com o local, ou seja, {u} = {ue} e {F} = {Fe}.
A aplicação do PTV significa impor uma variação de deformações e deslocamentos, ou seja, o
trabalho externo é gerado por uma variação de deslocamentos δue e o trabalho interno é
gerado por uma variação de deformação δε.
Portanto, a aplicação do PTV (Text=Tint) leva a:
2
∑F
i =1
e
i
⋅ δuie = ∫ σ ⋅ δε ⋅ dV
V
Reescrita de outra forma: {δu e } {F e } = ∫ δε ⋅ σ ⋅ dV
T
V
Substituindo o valor de σ dado na Equação (34) e observando que δε vale [B]{δue} ou
δε={δue}T[B]T, então:
{δu } {F }= ∫ {δu } [B]
e T
e T
e
T
{ }
⋅ E ⋅ [B ] u e dV
V
Como {δue}T é uma constante (valores de deslocamentos impostos nas extremidades), este
pode ser colocado fora da integral. Consequentemente, será cancelado com o mesmo valor
que aparece no primeiro membro. Assim:
{F }= ⎛⎜ ∫ [B]
T
e
⎝V
{ }
⎞
⋅ E ⋅ [B ]⋅ dV ⎟ ⋅ u e
⎠
Ou seja:
{F }= [k ]⋅ {u },
e
e
e
Sabendo que [B ] =
onde :
[k ] = ∫ [B]
T
e
V
1
[− 1 1] , então:
L
⋅ E ⋅ [B ]⋅ dV
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
28
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
[k ] = ∫ L1 ⎡⎢−11⎤⎥ ⋅ E ⋅ L1 [− 1 1]dV = LE ∫ ⎡⎢−11
e
⎣
V
2
⎦
V
Como: ∫ dV = A ⋅ L,
⎣
− 1⎤
E
dV = 2
⎥
1⎦
L
⎡ 1
[k ] = EA
L ⎢− 1
e
então :
⎣
V
⎡ 1 − 1⎤
⎢− 1 1 ⎥ ∫ dV
⎣
⎦v
− 1⎤
1 ⎥⎦
11 EXEMPLO NUMÉRICO
Proposta: obter a matriz de rigidez para a barra de treliça dada na Figura 19.
Solução:
A matriz [K] será 2×2, pois a barra tem dois graus de liberdade associado à mesma. A barra
suporta apenas forças axiais, então a relação entre força e deslocamento pode ser estabelecida
a partir da lei de Hooke, ou seja, a imposição de um deslocamento u gera uma força na barra
igual a k⋅u, onde k vale EA/L, conforme dedução a seguir.
σ =Eε
F
u
= E⋅
A
L
⇒
∴
F=
EA
⋅u
L
→
k=
EA
L
Tabela 2 – Obtenção da matriz de rigidez para a barra de treliça da Figura 19
1
Impondo
2
EA
= k 11
L
EA
F2 = −
= k 21
L
F1 = +
1⎫
⎬
⎩0⎭
{u} = ⎧⎨
F1
1
F2
L
2
1
Impondo
0⎫
⎬
⎩1⎭
{u} = ⎧⎨
F1
EA
= k 12
L
EA
F2 = +
= k 22
L
F1 = −
L
1
F2
Portanto, a matriz de rigidez de uma barra de treliça é dada pela Equação 35.
⎡ EA
⎢+ L
[K] = ⎢ EA
⎢−
⎣ L
EA ⎤
L ⎥
EA ⎥
⎥
+
L ⎦
−
12 FORMULAÇÃO
USANDO O PRINCÍPIO DOS
BARRAS DE PÓRTICO PLANO
Equação 35
TRABALHOS VIRTUAIS
PARA O CASO DE
Será adotado o polinômio de 3º. grau u(x) = a1 + a2 x + a3 x2 + a4 x3 para representar a
elástica da barra. São usados apenas quatro termos, pois efetivamente têm-se apenas quatro
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
29
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
graus de liberdade (Figura 23b), uma vez que dois graus de liberdade – os definidos na
direção axial – são independentes dos demais, que separadamente já foram tratados.
5
2
3
4
6
1
4
1
3
2
b) barra de pórtico – com apenas coordenadas de
flexão
a) barra de pórtico – completa
Figura 23 – Sistema de coordenadas para barra de pórtico plano.
Com a exclusão das coordenadas 1 e 4, por facilidade, os deslocamentos nodais receberão
índices seqüências de 1 a 4, respectivamente, para dois pares de translações perpendiculares
ao eixo da barra (coordenadas 2 e 5) e rotações de extremidades.
Lembrar, como já dito, que as coordenadas na direção axial da barra são independentes destas
quatro que estão vinculadas aos efeitos de flexão.
A aplicação do Princípio dos Trabalhos Virtuais para o caso da barra em questão admitindo
variações de deslocamentos (δue) e deformações (δε) resulta na Equação 36.
{δ u e }T { f } = ∫ (δε )σ dV
Equação 36
V
onde:
δue: variação dos deslocamentos nas coordenadas locais (do elemento)
δε: variação das deformações no elemento
Impondo-se as condições de contorno, pode-se escrever ue(x) em função dos deslocamentos
nodais da barra {ue}T = {ue1 ue2 ue3 ue4}, na forma u e ( x) = ∑ ϕiuie ou u ( x) = [ϕ ]{u e } . Assim,
obtem-se a Equação 37.
ϕ
ϕ
u e ( x) =
ϕ
1448
6447
3
2
(1 − 2 x 2 + 3 x 3 ) u1e +
L
L
344
6447
8
3 2 2 3 e
( 2 x − 3 x ) u3
L
L
+
ϕ
2444
644
47
8
2 2 1 3 e
( x − x + 2 x ) u2 +
L
L
Equação 37
4448
6447
1
1
(− x 2 + 2 x 3 ) u4e
L
L
{
u e ( x) = ϕ1u1e + ϕ 2u e2 + ϕ3u e3 + ϕ 4u e4 = [ϕ 1 ϕ 2 ϕ 3 ϕ 4 ] ⋅ u1e
u2e
u3e
u4e
}
T
Equação 38
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
30
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
O
dϕ
r
M
M
B
A
y
C D
Δdx
dx
Figura 24
Valendo-se das relações da elasticidade na flexão, obtidas com auxílio da Figura 24, obtém-se
o valor da deformação específica pela Equação 39 − (σ = E⋅ε, então, ε = σ/E = (M⋅y/EI).
ε=
Δdx M
=
y
dx
EI
Equação 39
Substituindo a equação diferencial
d 2u e
M
=−
aplicada à Equação 39, resulta na Equação
2
dx
EI
40.
ε=
⎞
M
y ⎛ d 2u e
⎜⎜ −
y=
EI ⎟⎟
2
EI
EI ⎝ dx
⎠
Equação 40
Lembrando que a variação de ε (deformação) implica na variação de ue (deslocamento), ou
vice-versa, então, a variação de ε, ou seja, δε, é dado pela Equação 41.
δε = − y
d 2 (δ u e )
dx 2
Equação 41
Calculando-se a tensão normal (σ) pela expressão clássica da Resistência dos Materiais e
substituindo o valor de M obtido pela Equação 40:
⎞
M
y ⎛ d 2u e
d 2u e
⎟
⎜
σ=
y = ⎜−
EI ⎟ = − E y
I
I ⎝ dx 2
dx 2
⎠
Equação 42
Substituindo os valores de δε (Equação 41) e σ (Equação 42), respectivamente na equação do
PTV (Equação 36), tem-se:
{δ u e }T { f }
=
∫V (δε )σ dV
=
⎛
d 2 δ ue
⎜
y
−
∫V ⎜
dx 2
⎝
(
) ⎞⎟⎛⎜ − E y d u
⎞
⎟ dV
dx ⎟⎠
2 e
⎟⎜
⎠⎝
2
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
31
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
Neste caso, o valor de E é constante e pode ser retirado da integral. Rearranjando o segundo
membro da equação anterior e observando que dV = dA dx, resulta a Equação 43.
⎛
⎞
L⎜
L
⎟
⎛ d 2 δ u e ⎞⎛ d 2u e ⎞ 2
2
⎟⎟⎜⎜
⎟ y dV = E ∫ ⎜ [ B' ] ∫ y dA ⎟dx = EI ∫ [ B' ]dx
{δ u } { f } = E ∫ ⎜⎜
2
2 ⎟
Equaçã
V
A
0
0
dx ⎠
123 ⎟
⎝14dx
⎝4
⎜
442⎠4
43
o 43
mom. inércia I ⎠
⎝
[ B ']
e T
(
)
Lembrar que I (momento de inércia) também é constante, pois o estudo restringe-se ao caso
de vigas de inércia constante.
Retornando à matriz [B’], esta pode ser reescrita com auxílio da relação a seguir, obtida a
partir da Equação 38, com a introdução da matriz chamada de [B]:
⎧u1e ⎫
⎪ e⎪
d 2u e
d2
⎪u ⎪
= 2 [ϕ1 ϕ 2 ϕ3 ϕ 4 ] ⎨ 2e ⎬ = [ B ]{u e }
2
dx
dx44424443 ⎪u3 ⎪
1
[ B]
⎪⎩u4e ⎪⎭
( )
d 2u e
d 2 δu e
e
e
=
[
B
]{
u
}
,
então
para
a
variação
de
u
pode-se
escrever:
= [ B]{δu e } . O
dx 2
dx 2
d 2 δu e
= {δu e }T [ B]T
mesmo para a forma transposta:
2
dx
Se
( )
Assim, a matriz [B’], com base na equação anterior, torna-se:
(
⎛ d 2 δ ue
[ B' ] = ⎜⎜
2
⎝ dx
) ⎞⎟⎛⎜ d u
⎞
δ 42
u e }T [ B]T [ B]{u e }
⎟⎜ dx 2 ⎟⎟ = {1
43
⎠
⎠⎝
[ B ]{δu e }
2
e
Substituindo estes valores na Equação 43 tem-se:
{δ u e }T { f } = EI ∫0 {δ u e }T [ B ]T [ B ]{u e }dx ou
L
L
{ f } = EI ∫0 [ B]T [ B]dx {u e } ⇒ { f } = [k e ]{u e }
Como {δ ue}T e {ue} são vetores formados por constantes (não dependem de x), estes podem
ser colocados fora da integral. Conseqüentemente, {δue}T será cancelado. Resulta, portanto, o
sistema de equações [ke]{ue} = {f}, onde a matriz [B] vale:
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
32
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
⎡⎛ 6 12 ⎞ ⎛ 4 6 ⎞ ⎛ 6 12 ⎞ ⎛ 2 6 ⎞⎤
[ B] = ⎢⎜ − 2 + 3 x ⎟ ⎜ − + 2 x ⎟ ⎜ 2 − 3 x ⎟ ⎜ − + 2 x ⎟⎥
L ⎠ ⎝ L L ⎠ ⎝L
L ⎠ ⎝ L L ⎠⎦
⎣⎝ L
Finalmente tem-se a matriz de rigidez no sistema local de um elemento de pórtico plano:
6 L − 12 6 L ⎤
⎡ 12
⎢ 6 L 4 L2 − 6 L 2 L2 ⎥
EI
⎥
[k e ] = EI ∫ [ B]T [ B] dx = 3 ⎢
⎢
⎥
−
−
−
12
6
12
6
L
L
L
⎢
⎥
2
− 6 L 4 L2 ⎦
⎣ 6L 2L
Esta matriz é dada no sistema de coordenadas locais (da barra – coordenadas paralelas e
perpendiculares ao eixo da barra). Portanto, é necessário fazer a transformação de
T
informações do sistema de coordenadas local para o global, ou seja, K e = [T ] k e [T ] . A
matriz resultante é fornecida no Apêndice I.
[ ]
[ ]
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
33
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
13 FORMULAÇÃO PELO MÉTODO DOS ELEMENTOS FINITOS
13.1 Generalidades
A solução de problemas usando o método dos elementos finitos sempre apresentará como
solução um sistema de equações do tipo [K]{u} = {F}, ou seja, (matriz de rigidez) × (vetor de
deslocamentos) = (vetor de forças). Esta matriz de rigidez pode ser obtida usando diferentes
procedimentos, basicamente de três formas usando o chamado método direto, método
variacional e método dos resíduos ponderados.
O Método Direto é usado para problemas que envolvem elementos mais simples onde as
relações entre parâmetros são facilmente conhecidas, e as equações da estática são suficientes
para solução do problema. O Método Variacional é possível de ser usado, desde que o
problema apresente um funcional, baseado no princípio da mínima energia. O Método de
Rayleigh-Ritz está diretamente ligado à solução variacional. O Método dos Resíduos
Ponderados é mais abrangente e baseia-se na diferença entre os valores exatos e aproximados,
ou seja, o erro é ponderado no domínio usando diferentes critérios de distribuição dentre os
quais o Método de Galerkin é o mais empregado.
13.2 Método direto
No chamado método direto, a matriz de rigidez é gerada de forma direta – como o próprio
nome sugere, isto é, não é utilizado nenhum procedimento matemático sofisticado associado a
equações diferencias complexas. O simples equilíbrio estático e o uso de relações da teoria da
elasticidade serão suficientes para resolver o problema. Obviamente que o método é limitado
por abranger apenas problemas mais simples. No entanto, este procedimento servirá como
base para a montagem de problemas mais complexos, especialmente quanto ao
posicionamento dos coeficientes de rigidez na matriz [K], obtidos a partir da matriz de rigidez
dos elementos individualmente. Os exemplos apresentados anteriormente representam
aplicações do método direto.
Para ilustrar a sua aplicação será apresentado mais um exemplo de uma estrutura simples
solicitada apenas por forças axiais, conforme a Figura 25, para a qual se calculam os
deslocamentos e as tensões em alguns pontos da barra.
2P
A1
240 cm
E = 2100 kN/cm2
P
A2
A1 = 1200 cm2
100 cm
A2 = 400 cm2
Figura 25
Para solução será adotado que a peça é dividida em três elementos (4 nós), como indica a
Figura 26.
1
2
120
3
120
Figura 26
Considerando um material que segue a lei de Hooke:
4
100 cm
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
34
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
σ =Eε
⇒
F
EA
Δl
=E
⇒ F=
Δl
l
l
A
Observar que EA/ℓ é o coeficiente de rigidez que relaciona deslocamento com força em um
mesmo ponto e passará a ser representado pela letra k. Usando a notação usual para
deslocamento através da letra u e entendendo que neste caso u representa a variação de
comprimento Δ ℓ e vale |uj – uk| (j é o nó inicial e k é o nó final do elemento i), pode ser
escrito:
F =ku
∴
k=
EA
l
Isolando cada elemento e avaliando os esforços atuantes (Fji e Fki) e os deslocamentos (uj e uk)
nas extremidades, tem-se a representação dada na Figura 27.
Fji
nó k
nó j
elemento i
uj
Fki
uk
Figura 27
Pelo equilíbrio estático do elemento i da Figura 27 têm-se duas equações que relacionam as
forças atuantes nas extremidades com os deslocamentos:
Pelo equilíbrio de forças sobre cada elemento (somatório de forças na direção do elemento),
pode-se concluir que Fji = − Fki.
Para cada extremidade da barra i pode-se afirmar que Fji = ki (uj – uk) e Fki = ki (uk – uj). As
duas expressões podem ser escritas matricialmente na forma:
⎡ ki
⎢
⎣− ki
− ki ⎤ ⎧u j ⎫ ⎧ F ji ⎫
⎥⎨ ⎬ = ⎨ ⎬
ki ⎦ ⎩uk ⎭ ⎩ Fki ⎭
Agora considerando o domínio completo, representado na Figura 28, onde toda a estrutura é
subdividida em elementos e estão indicados os esforços internos atuando sobre os nós.
Observar que estas forças estão no sentido contrário àquelas aplicadas sobre os elementos da
Figura 27, pois nesta representação atuam sobre os nós.
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
35
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
nó 1
nó 2
nó 3
elemento 2
elemento 1
R
Fk1
Fj1
elemento 3
Fk2
Fj2
nó 4
2P
Fj3
Fk3
P
x
Figura 28 – Solicitações sobre os nós ao longo da estrutura representada pelos elementos
Com esta discretização têm-se quatro valores de u (um valor para cada nó), que podem ser
representados pelo vetor: {u}T = {u1 u2 u3 u4}.
Para a estrutura podem ser estabelecidas quatro equações de equilíbrio estático (uma para
cada nó), seguindo a orientação dada pela coordenada x:
nó
Equilíbrio de forças
Equilíbrio de forças em função de un
1
R – Fj1 = 0
k1 (u1 – u2) = R
2
-Fk1 – Fj2 = 0
3
-Fk2 – Fj3 + 2P = 0
4
-Fk3 + P = 0
k1 (u2 – u1) + k2 (u2 – u3) = 0
⇒
k2 (u3 – u2) + k3 (u3 – u4) = 2P
k3 (u4 – u3) = P
Reescrevendo estas equações de forma matricial, tem-se:
⎡ k1
⎢
⎢− k1
⎢ 0
⎢
⎢⎣ 0
− k1
0
k1 + k 2
− k2
− k2
k 2 + k3
0
− k3
0 ⎤ ⎧u1 ⎫ ⎧ R ⎫
⎥⎪ ⎪ ⎪ ⎪
0 ⎥ ⎪u 2 ⎪ ⎪ 0 ⎪
⎨ ⎬=⎨ ⎬
− k 3 ⎥ ⎪u 3 ⎪ ⎪2 P ⎪
⎥
k 3 ⎥⎦ ⎪⎩u 4 ⎪⎭ ⎪⎩ P ⎪⎭
O sistema de equações anterior representa toda a estrutura independentemente das condições
de contorno. Para a solução efetiva é necessário impor que u1 vale zero. Em outras palavras, a
solução vem de um sistema de equações formado por três equações e três incógnitas pela
eliminação da primeira linha e primeira coluna do sistema (u1 = 0).
Outra forma de montar este sistema de equações, tendo em vista a automatização
computacional, é usar os fundamentos do MEF considerando a discretização estrutural onde u
é representado por:
m
u = ∑ u ( e ) = u (1) + u ( 2 ) + u ( 3)
e =1
Esta é a equação de compatibilidade de deslocamentos dos nós. Desta forma, a matriz de
rigidez é gerada pela sobreposição das contribuições de cada elemento, o que torna o processo
mais padronizado (automatizado) facilitando a manipulação computacional.
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
36
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
13.3 Método da energia
Basicamente duas formas de energia se manifestam no estudo da estática e da estabilidade:
energia de deformação e energia potencial das cargas (ações).
Pela conservação da energia, uma das manifestações de energia se transforma na outra. A
aplicação de forças produz um trabalho (redução do potencial das forças) que se acumula na
forma de energia de deformação.
O trabalho das forças representa uma transformação de energia. No conceito de trabalho
positivo, ou seja, força e deslocamento no mesmo sentido, significa produto escalar positivo.
Porém, no cálculo do potencial de forças deve-se interpretar isto com um sinal negativo, tendo
em vista que isto representa uma redução do potencial da força, ou seja, é a passagem da força
na posição deslocada para a posição indeformada. Representa a transformação da energia de
deformação em energia potencial.
Energia de deformação ⇒ U
Potencial de forças atuantes ⇒ Ω
Π=U+Ω
Enegia potencial total ⇒ Π
Pelo princípio da conservação de energia, pode-se afirmar que:
δΠ = δU + δΩ ⇒ condição estacionária
Potencial de forças atuantes (Ω): corresponde ao somatório do produto escalar das forças
aplicadas sobre os nós pelos correspondentes deslocamentos. Lembrar que se a força estiver
aplicada sobre o elemento, deve-se fazer uma equivalência de trabalho transformando estas
ações em concentradas sobre os nós. Como as forças são constantes na passagem da posição
inicial para a final, o potencial das forças resultará em:
F
n
Ω = ∑ Fi × ui
F
i +1
u
Energia de deformação (U):
F
U=
1
(σ xε x + σ yε y + σ zε z + τ xyγ xy + τ yzγ yz + τ xzγ xz ) dV
2 ∫V
F
u
Observar que neste caso existe o fator ½, uma vez que a deformação e a aplicação das ações
ocorrem de forma gradativa.
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
37
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
De forma matricial:
U=
Se [σ] = [C] {ε} então:
1
[σ ] T ⋅ {ε } dV
∫
V
2
U=
1
{
ε }T [C ]{ε } dV
∫
2 V
Lembrando que {ε} = [B]{ue}, ou seja, [B] é a matriz que transforma deslocamentos em
deformações, então:
U=
1
{u e }T [ B]T [C ][ B]{u e } dV
∫
V
2
Minimizando a energia de deformação (condição estacionária) em relação aos deslocamentos
∂U
obtem-se a matriz [ke]:
{ue}, ou seja, fazendo-se
e
∂ui
[k e ] = ∫ [B]T [C ][ B]dV
V
Esta expressão é idêntica àquela encontrada quando se utilizou o princípio dos trabalhos
virtuais.
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
APÊNDICES
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
39
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
APÊNDICE A – Componentes de tensões para o caso geral Y
Estado plano de tensões
σyy
σyy
τyz
τzy
τyx
τyx
τzx
σzz
τxy
σxx
τxy
σxx
σxx
τxy
τxz
τyx
X
σyy
Z
Elemento deformado para o caso plano:
∂u
dy
∂y
Y, v
∂v
dy
∂y
γ2
dy
∂v
dx
∂x
γ1
∂u
dx
∂x
dx
X, u
Pode-se concluir pela figura anterior que:
εx =
∂u
∂x
εy =
∂v
∂y
γ xy = γ 1 + γ 2 =
∂u ∂v
+
∂y ∂x
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
40
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
APÊNDICE B – Relações entre tensões e deformações {σ} = [C] {ε}
⎧σ xx ⎫
⎪ ⎪
E
⎨σ yy ⎬ =
2
⎪ ⎪ 1−υ
⎩τ xy ⎭
Tensões planas
(σzz = τxz = τyz = 0)
Deformações
planas
(εzz = γxz = γyz = 0)
Tridimensional
y
⎡1 υ
0 ⎤ ⎧ε xx ⎫
⎢
⎥⎪ ⎪
⎢υ 1
0 ⎥ ⎨ε yy ⎬
⎢
1 - υ ⎥⎪ ⎪
⎢0 0
⎥ γ xy
2 ⎦⎩ ⎭
⎣
⎡1 − υ υ
⎧σ xx ⎫
0 ⎤ ⎧ε xx ⎫
⎢
⎥⎪ ⎪
⎪ ⎪
E
⎢
υ 1- υ
0 ⎥ ⎨ε yy ⎬
⎨σ yy ⎬ =
(
)(
)
⎢
+
−
1
υ
1
2
υ
1 - 2υ ⎥ ⎪ ⎪
⎪ ⎪
τ
0
0
⎢
⎥ γ xy
⎩ xy ⎭
2 ⎦⎩ ⎭
⎣
z
x
y
z
x
υ
υ
0
0
0 ⎤
⎡1 − υ
⎧σ xx ⎫
⎢
⎥ ⎧ε xx ⎫
υ
0
0
0 ⎥⎪ ⎪
⎪ ⎪
⎢ υ 1− υ
ε
⎪σ yy ⎪
⎢ υ
υ 1− υ
0
0
0 ⎥ ⎪ yy ⎪
⎪σ ⎪
⎢
⎥ ⎪ε ⎪
1 - 2υ
E
⎪ zz ⎪
⎢
⎥ ⎪⎨ zz ⎪⎬
0
0
0
0
0
=
⎨ ⎬
2
⎥ ⎪γ xy ⎪
⎪τ xy ⎪ (1 + υ )(1 − 2 υ ) ⎢
1
2υ
⎢
⎥
⎪τ ⎪
0
0
0
0
0 ⎥ ⎪γ ⎪
⎢
2
⎪ yz ⎪
⎪ yz ⎪
⎢
1 - 2υ ⎥ ⎪⎩γ xz ⎪⎭
⎪⎩ τ xz ⎪⎭
0
0
0
0
⎢ 0
⎥
2 ⎦
⎣
υ: coeficiente de Poisson
E: módulo de deformação longitudinal
∂u
∂x
∂u ∂v
=
+
∂y ∂x
∂v
∂y
∂v ∂w
=
+
∂z ∂y
∂w
∂z
∂u ∂w
=
+
∂z ∂x
ε xx =
ε yy =
ε zz =
γ xy
γ yz
γ xz
As relações entre deformações e tensões são dadas pela lei de Hooke generalizada (material
isótropo):
[
]
[
]
[
]
1
σ xx −ν (σ yy + σ zz )
E
1
ε yy = σ yy −ν (σ xx + σ zz )
E
1
ε zz = σ zz −ν (σ xx + σ yy )
E
E
1
G=
γ xy = τ xy
G
2(1 + ν )
ε xx =
γ yz =
1
τ yz
G
γ zx =
1
τ zx
G
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
41
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
APÊNDICE C – Deformações na flexão Será tomado como referência a Figura 29 para indicar os parâmetros utilizados nesta
demonstração.
O
dϕ
dϕ
dx
r
r
M
x
M
y
y=v
dx
ds
B
A
ϕ
C D
ϕ
ds
Δdx
v = f(x)
(b)
(a)
dx
dv
Figura 29 – Elemento de viga deformado por flexão.
Para a análise em questão considera-se que os ângulos são muito pequenos podendo-se
considerar tg ϕ ≅ ϕ.
Δdx
Por Hooke: σ = Eε ⇒ σ = E
dx
Substituindo na equação da flexão:
σ=
M
M × dx
Δdx
⇒ Δdx =
y=E
y
I
dx
EI
(I)
Pelo triângulo ΔOAB:
r dϕ = dx ⇒ dϕ =
dx
r
(II)
Pelo triângulo ΔBCD:
y dϕ = Δdx ⇒ dϕ =
Δdx
y
(III)
Substiuindo (I) em (III):
dϕ =
M × dx
EI
Igualando (II) e (IV) tem-se:
(IV)
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
42
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
dx M
=
dx
r
EI
A relação
⇒
1 M
=
r EI
(V)
1
é chamada de curvatura e será representada pela letra grega κ (kappa).
r
De (II) tem-se:
1 dϕ
=
, que igualada a (V) resulta:
r dx
dϕ M
=
dx EI
(VI)
Pelo detalhe mostrado na parte inferior esquerda da Figura 29 pode-se escrever que
dϕ d 2 v
dv
tgϕ ≅ ϕ =
, ou seja:
.
=
dx
dx dx 2
Igualando-se com (VI) tem-se:
d 2v
M
=−
2
dx
EI
(VII)
O sinal negativo introduzido na expressão acima é necessário para ajustar os sinais de acordo
com as explanações a seguir.
Este sinal é introduzido para compatibilizar o resultado envolvendo o momento fletor com o
sentido indicado na Figura 29(b) – tração embaixo. Para justificar este sinal tem-se como
referência a Figura 30 que representa o elemento de comprimento dx com momentos fletores
tracionando a parte inferior do elemento, de acordo com a referência estabelecida na Figura
29.
Figura 30 – Referência para avaliação das derivadas associadas à elástica de um elemento de viga de
comprimento dx, submetido à ação dos momentos fletores que causam tração na parte inferior.
A justificativa do sinal baseia-se na análise da segunda derivada da função v. Inicialmente é
considerado o trecho à esquerda do centro do elemento, onde a função v é positiva e
crescente. Sendo crescente, os valores dv/dx serão positivos variando de um valor máximo até
zero, o que indica que dv/dx é decrescente. Sendo decrescente, os valores de d 2v/dx2 serão
negativos. Assim, se o momento M que gera a curvatura indicada na Figura 30 foi
considerado positivo, então d2v/dx2 e M estão com sinais contrários. Para compensar tal fato,
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
43
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
adota-se o sinal negativo na expressão (VI). A Figura 31 ilustra as funções v, dv/dx e d 2v/dx2
no intervalo de 0 a dx/2.
0
dx/2
0
-
dx/2
0
+
+
v
dx/2
d2v/dx2
dv/dx
Figura 31 – Análise do sinal da 2a. derivada de v para x de 0 a dx/2.
A mesma análise quando aplicada ao intervalo dx/2 a dx revela que d 2v/dx2 também é
negativo como indicado na Figura 32. Portanto, a equação (VI) é válida para o trecho 0 a dx,
de fato contendo o sinal negativo.
dx/2
-
dx
dx/2
+
v
dx/2
dx
dx
d2v/dx2
dv/dx
Figura 32 – Análise do sinal da 2a. derivada de v para x de dx/2 a dx.
Da Figura 31 e da Figura 32 pode-se concluir que a terceira derivada de v em relação a x, que
corresponde à força cortante (V=dM/dx= d 3v/dx3) terá sinal positivo de 0 a dx/2 e sinal
negativo de dx/2 a dx. Isto significa que a força cortante é positiva à esquerda quando atua de
cima para baixo e, de baixo para cima, do lado direito do elemento, exigindo a troca do sinal.
V
V
v
M
M
M
M
V
v
a) sinal positivo de v para baixo
V
b) sinal positivo de v para cima
Figura 33 – Convenção de sinais de momentos fletores e forças cortantes quando se utiliza a expressão
d 2v
M
d 3v
V
=
−
=±
e
(+ à esquerda e – à direita)
2
3
dx
EI
dx
EI
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
44
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
Para ilustrar ou elucidar o exposto anteriormente será usada a equação da elástica de uma viga
bi-apoiada solicitada por força uniformemente distribuída, mostrada na Figura 34. No caso, é
considerado que a elástica é positiva quando desloca de cima para baixo. Na Figura 34 estão
mostrados os sentidos do momento fletor e da força cortante envolvidos no caso, conhecidos
da estática: tração positiva na fibra inferior e forças cortantes positivas se têm sentido horário.
q
q
v
M1
M2
V2
V1
Equação da elástica: v =
(
)
q
x 4 − 2 l x 3 + l 3 x (ver dedução no Erro! Fonte de referência não
24 ⋅ EI
encontrada.)
Figura 34 – Viga bi-apoiada com força uniformemente distribuída
Como é sabido, as derivadas de diferentes ordens da equação da elástica estão associadas a
parâmetros da viga, tal como se indica a seguir:
v → elástica;
v ' → rotação;
v ''' → força cortante;
v '' → momento fletor
v '''' → força distribuída
O objetivo agora é saber qual o verdadeiro sentido das solicitações ao se aplicar as relações
diferenciais. Será considerado o ponto central da viga como referência para o momento fletor
(x = ℓ/2) e o apoio da esquerda (x = 0) para a força cortante.
Calculando a segunda derivada da elástica para associá-la ao momento fletor (EIv’’= −M),
para x = ℓ/2 tem-se:
(
q
EIv = x 2 − l x
2
''
)
x =l / 2
ql 2
=−
8
Pelo resultado do estudo anterior, quando se troca o sinal obtido pela segunda derivada, este
provocará tração na fibra inferior. Então, trocando-se o sinal do resultado anterior para
adaptar-se à expressão EIv’’= −M, tem-se
M = − EIv'' =
ql 2
8
Como a convenção da elástica está como indicado na Figura 33a, então pode ser afirmado que
o momento fletor traciona em baixo.
Para a força cortante emprega-se a terceira derivada da elástica, no caso, calculando-a para
x=0, tem-se:
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
45
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
EIv ''' =
q
(2 x − l ) = − ql
2
2
x =0
No caso da avaliação da força cortante no apoio esquerdo (x=0), resultou em valor negativo
pela terceira derivada, devendo ser mantido pela conclusão ilustrada na Figura 33a. Então, se
a esquerda o sinal é negativo, isto significa que o binário produzido pela força cortante é
horário. A conclusão é a mesma quando se analisa pela direita. O sinal obtido é negativo e
precisa ser trocado pelas conclusões da análise de sinais. Isto implica, também, em binário
horário sobre o elemento, resultado na situação indicada na Figura 35.
V=q ℓ/2
V=q ℓ/2
x=0
v
Figura 35 – Força cortante no apoio da esquerda da viga (x=0)
Finalmente, é possível obter a força q uniformemente distribuída pela relação que envolve a
quarta derivada, ou seja, EIv '''' = q . Resulta em valor constante qualquer que seja x. Observe
que o valor de q positivo é de cima para baixo, o que coincide com o valor dado pela
expressão. Assim, não há troca de sinal.
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
46
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
APÊNDICE D – Elástica para viga bi‐apoiada com força uniformemente distribuída Para a viga da Figura 36 é possível expressar, por meio do simples equilíbrio estático, o
momento fletor pela equação:
M=
qLx qx 2
−
2
2
Como foi deduzido no APÊNDICE C, o momento fletor está relacionado com a segunda
derivada da elástica, ou seja, EI d 2v dx 2 = − M . Então, as duas expressões anteriores podem
ser igualadas:
EI
d 2v
qLx qx 2
=
=
−
+
EIv
"
dx 2
2
2
Figura 36 – Viga bi-apoiada com força uniformemente distribuída
Para conhecer a primeira derivada de v, faz-se a integração da equação anterior:
EIv' = −
qLx 2 qx 3
+
+ c1 .
4
6
Para a viga em questão, pode ser afirmado que a inclinação da curva elástica no centro é nula,
pois a viga é elástica e geometricamente simétrica e submetida a um carregamento simétrico.
Então, pode ser afirmado que para x = L/2, v’ = 0. Isto permite o cálculo de c1, ou seja:
qL3
c1 =
.
24
Este valor levado à equação EIv’, resulta em: EIv' = −
qLx 2 qx 3 qL3
+
+
.
4
6
24
Integrando a equação anterior, chega-se à curva elástica v: EIv = −
qLx 3 qx 4 qL3 x
+
+
+ c2 .
12
24
24
Impondo-se a condição de que a elástica é nula em x = 0, obtêm-se o valor de c2, que no caso
resulta igual a zero. O mesmo poderia ser obtido se fosse empregada a condição de v = 0 em x
= L. Portanto, resulta a equação da elástica:
v=
(
q
x 4 − 2 L x 3 + L3 x
24 ⋅ EI
)
Esta equação permite calcular o valor do deslocamento vertical em qualquer ponto da
estrutura ao longo do comprimento L. Em particular, no centro do vão tem-se o valor da
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
47
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
máxima flecha para a viga simplesmente apoiada com força uniformemente distribuída.
Então, para x = L/2 tem-se:
vmax =
4
3
⎞ 5qL4
q ⎛⎜ ⎛ L ⎞
⎛L⎞
3⎛ L ⎞⎟
L
L
−
+
2
⎜ ⎟⎟ =
⎜ ⎟
⎜ ⎟
24 ⋅ EI ⎜⎝ ⎝ 2 ⎠
⎝2⎠
⎝ 2 ⎠ ⎠ 384 EI
Exercício proposto: deduza a equação da linha elástica para uma viga engastada-livre com
qx 2
x 2 − 4 L x + 6 L2 .
força uniformemente distribuída. Resposta: v =
24 ⋅ EI
(
)
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
APÊNDICE E – Regras para diferenciais (Informações obtidas em Murray R. Spiegel)
d(u ± v ± w ± ...) = du ± dv ± dw ± ...
d(uv) = udv + vdu
d(u/v) = (vdu + udv)/v2
d(un) = n un-1 du
d(sen u) = cos u du
d(cos u) = -sen u du
d(ln u) = du/u
Regras de diferenciação:
dy dy du
=
×
dx du dx
du
1
=
dx dx / du
d
(ku ) = k du
dx
dx
dy dy / du
=
dx dx / du
d
(uv ) = u dv + v du
dx
dx
dx
d
(cx ) = c
dx
d ⎛ u ⎞ v(du / dx) − u (dv / dx)
⎜ ⎟=
dx ⎝ v ⎠
v2
d2
d 2v
du dv
d 2u
(
)
uv
=
u
+
2
+
v
dx 2
dx 2
dx dx
dx 2
∂2 f
∂ ⎛ ∂f ⎞
= ⎜ ⎟
2
∂x
∂x ⎝ ∂x ⎠
∂2 f
∂ ⎛ ∂f ⎞
= ⎜⎜ ⎟⎟
2
∂y
∂y ⎝ ∂y ⎠
∂2 f
∂ ⎛ ∂f ⎞
= ⎜⎜ ⎟⎟
∂x∂y ∂x ⎝ ∂y ⎠
∂2 f
∂ ⎛ ∂f ⎞
= ⎜ ⎟
∂x∂y ∂y ⎝ ∂x ⎠
df =
∂f
∂f
dx +
dy
∂x
∂y
n
df = ∑
i =1
∂f
dxi
∂xi
onde
dx = Δx e dy = Δy
f = f(x1 , x2 ,....xn )
48
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
49
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
APÊNDICE F – Um problema específico Dada a função: v(ξ) = a1 + a2 ξ+ a3 ξ2 + a4 ξ3 sendo ξ = x/L ∴ dξ = dx / L
Isto significa que:
Se: y =
dv dv dξ 1 dv
=
×
= ×
dx dξ dx L dξ
dv
dy
d ⎛ dv ⎞ d 2 v
então y ' =
=
⎜ ⎟=
dx
dx ⎝ dx ⎠ dx 2
dx
d 2 v d ⎛ dv ⎞
dv
1 dv
por
tem - se :
×
= ⎜ ⎟ . Substituindo dx por L × dξ e
2
dx
dx ⎝ dx ⎠
dx
L dξ
d 2v
d ⎛ 1 dv ⎞ 1 d 2 v
⎜ ×
⎟= ×
=
dx 2 L × dξ ⎜⎝ L dξ ⎟⎠ L2 dξ 2
Refazendo os mesmo cálculos anteriores, quanto à determinação da segunda derivada de v em
função de x, usando a regra da cadeia:
Regra da cadeia : (
d (uv )
dv
du
=u +v ):
dx
dx
dx
d 2 v d ⎛ dv ⎞ d ⎛ 1 dv ⎞
1
⎟⎟. Fazendo u =
= ⎜ ⎟ = ⎜⎜ ×
2
dx
dx ⎝ dx ⎠ dx ⎝ L dξ ⎠
L
e w=
dv
tem - se :
dξ
⎞ 1 d 2v
d 2 v ⎛ 1 d ⎛ dv ⎞ ⎞ ⎛ dv
⎜
⎟
⎜
⎟
⎜
⎟⎟ = 2 × 2 (lembrar que dx = L dξ )
=
×
+
×
0
dx 2 ⎜⎝ L dx ⎜⎝ dξ ⎟⎠ ⎟⎠ ⎜⎝ dξ
⎠ L dξ
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
APÊNDICE G – Quadratura de Gauss No de pontos
2
3
4
5
Posição (ξ)
Fator de ponderação (ω)
-0,577350269
1,00000000
0,577350269
1,00000000
-0,774596669
0,55555556
0
0,88888889
0,774596669
0,55555556
-0,861136312
0,3478548
-0,339981044
0,6521452
0,339981044
0,6521452
0,861136312
0,3478548
-0,906179846
0,2369269
-0,538469310
0,4786287
0
0,5688889
0,538469310
0,4786287
0,906179846
0,2369269
50
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
51
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
APÊNDICE H – Galerkin Dedução das matrizes dos elementos do exemplo utilizando o método dos resíduos
ponderados - Galerkin
⎛ d φr( e ) d φt( e )
⎞
K = ∫ ⎜⎜
⋅
+ φr φt ⎟⎟ i, j dξ
dξ
dξ
⎠
0 ⎝
dφ1( e )
dφ2( e ) 1
1
e
=−
=
dξ
h
dξ
h
Para r = 1 e t = 1 :
h
(e)
ij
(r, t = 1, 2)
ξ
ξ ⎞
1
⎛ 1
= ∫ ⎜ (− )(− ) + (1 − )(1 − ) ⎟ i, j dξ = ∫
h
h
h
h ⎠
0
0 ⎝
h
K
(e)
ij
h
⎛1
2ξ ξ 2 ⎞
1 h
⎜⎜ 2 + 1 −
+ 2 ⎟⎟ dξ = +
h h ⎠i, j
h 3
⎝h
Para r = 2 e t = 2 :
ξ ξ ⎞
⎛ 1 1
K ij(e) = ∫ ⎜ ( )( ) + ( )( ) ⎟ i, j dξ = ∫
h h
h h ⎠
0 ⎝
0
h
h
h
⎛ξ
⎛ 1 ξ2 ⎞
ξ3 ⎞
1 h
⎜⎜ 2 + 2 ⎟⎟ dξ = ⎜⎜ 2 + 2 ⎟⎟ = +
h ⎠i, j
⎝ h 3h ⎠ 0 h 3
⎝h
Para r = 1 e t = 2 :
ξ ξ ⎞
⎛ 1 1
= ∫ ⎜ (− )( ) + (1 − )( ) ⎟ i, j dξ = ∫
h h
h h ⎠
0
0 ⎝
h
K
(e)
ij
h
h
⎛ 1 ξ ξ2 ⎞
⎜⎜ − 2 + − 2 ⎟⎟ dξ =
h h ⎠i, j
⎝ h
⎛ ξ ξ2 ξ3 ⎞
1 h
= ⎜⎜ − 2 +
− 2 ⎟⎟ = − +
h 6
⎝ h 2h 3h ⎠ 0
Para r = 2 e t = 1 : semelhante ao caso anterior
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
52
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
APÊNDICE I – Uso do ANSYS Algumas dicas para combinar carregamentos com diferentes fatores de majoração:
1) Insira as condições de carregamento referentes ao peso próprio
2) Execute na linha de comando: solve
3) Insira as condições de carregamento referentes a sobrecarga
4) Execute na linha de comando: solve
5) Insira as condições de carregamento referentes ao vento
6) Execute na linha de comando: solve
Após estes passos, vc terá 3 load steps.
Siga os passos abaixo:
7) Main menu - General Postproc - Load case - Create load case - results - ok
Insira o nº do caso de carga = 2 e escolha loadstep = 2 - Apply
Repita o processo para o caso de carga 3 e de ok
8) Main menu - General Postproc - Load case - Calc Options - Scale factor
Escolha load case = 2 e insira o fator de escala para a sobrecarga - apply
Repita o processo para o load case 3 e de ok
9) Main menu - General Postproc - Read Results - By load step
Escolher o loadstep 1 e aplicar fator de escala para o peso próprio - ok
10) Main menu - General Postproc - Load case - Add
Escolher load case 2 - apply
Escolher load case 3 - ok
11) Main menu - General Postproc - Write results (usar write results para gerar o
Load step 4)
Escolher Load step 4 e ok
O load step 4 será os resultados combinados.
Notas de Aula: Método dos Elementos Finitos – Prof. Francisco A. R. Gesualdo
53
Programa de Pós-Graduação em Engenharia Civil – FECIV – UFU
APÊNDICE J – Matriz de Rigidez do elemento no Sistema de Coordenadas Global [Ke] [K ] = [T ] [k ][T ]
T
e
[K ]
e
⎡cos θ − senθ
⎢
⎢ senθ cos θ
⎢
⎢ 0
0
=⎢
⎢ 0
0
⎢
⎢ 0
0
⎢
⎢ 0
0
⎣
0
0
1
0
0
0
⎡ EA
0
0
0⎤ ⎢ l
⎥ ⎢⎢ 0
0
0
0⎥ ⎢
⎥ ⎢
0
0
0⎥ ⎢ 0
⎥× ⎢
cos θ − senθ 0⎥ ⎢ − EA
⎥ ⎢
l
senθ cos θ 0⎥ ⎢
⎥ ⎢ 0
0
0
1⎥⎦ ⎢
⎢ 0
⎣
0
0
12 EI
l3
6 EI
l2
6 EI
l2
4 EI
l
0
0
−
12 EI
l3
6 EI
l2
6 EI
l2
2 EI
l
−
e
−
EA
l
0
0
EA
l
0
0
0
0
12 EI
− 3
l
6 EI
− 2
l
6 EI
l2
2 EI
l
0
0
12 EI
l3
6 EI
− 2
l
6 EI
l2
4 EI
l
−
⎤
⎥ ⎡ cos θ senθ
⎥ ⎢
⎥ ⎢
⎥ ⎢− senθ cos θ
⎥ ⎢
0
⎥ ⎢ 0
×
⎥ ⎢
0
⎥ ⎢ 0
⎥ ⎢
0
⎥ ⎢ 0
⎥ ⎢
0
⎥ ⎣ 0
⎥
⎦
0
0
1
0
0
0
0⎤
⎥
0
0
0⎥
⎥
0
0
0⎥
⎥
cos θ senθ 0⎥
⎥
− senθ cos θ 0⎥
⎥
0
0
1⎥⎦
0
0
12 EI
⎡ EA
12 EI
⎛ EA 12 EI ⎞
⎞
⎛ EA 12 EI ⎞
⎛ EA
2
2
6 EI
cos 2 θ + 3 sen 2θ ⎟ − ⎜
− 3 ⎟ senθ cos θ
− 3 ⎟ senθ cos θ − 6 EI senθ ⎤
−⎜
⎜
sen
−
θ
⎢ l cos θ + l 3 sen θ
2
⎥
l
l
l
l ⎠
l
l ⎠
⎠
⎝
⎝
⎝
l
l2
⎢ ⎛ EA 12 EI ⎞
⎥
EA
EI
12
EA
EI
EA
EI
12
12
6 EI
⎛
⎞
⎛
⎞ 6 EI
⎢ ⎜
− 3 ⎟ senθ cos θ
sen 2θ + 3 cos 2 θ
sen 2θ + 3 cos 2 θ ⎟
−⎜
− 3 ⎟ senθ cos θ − ⎜
cos
⎥
θ
θ
cos
2
2
l ⎠
⎢ ⎝ l
l
l
l ⎠
l
l
l
⎝ l
⎝ l
⎠
⎥
⎢
6 EI
6 EI
4
EI
EI
2
EI
EI
6
6
⎥
cos θ
− 2 senθ
senθ
− 2 cos θ
⎢
2
2
⎥
l
l
e
l
l
l
l
K =⎢
⎥
12
EA
EI
12
EA
EI
EI
6
6
EI
⎛
⎞
⎛
⎞
EA
12 EI
⎛ EA 12 EI ⎞
⎢− ⎜
cos 2 θ + 3 sen 2θ ⎟ − ⎜
− 3 ⎟ senθ cos θ
senθ ⎥
senθ
cos 2 θ + 3 sen 2θ
sen
cos
θ
θ
−
⎜
⎟
2
2
⎢ ⎝ l
l
l ⎠
l
l
⎠
⎝ l
l
l
l3 ⎠
⎥
⎝ l
⎢ ⎛ EA 12 EI ⎞
EI
6
6
EI
EA
EI
12 EI
12
⎛
⎞
⎛ EA
2
2 ⎞ −
EA
EI
12
− 2 cos θ⎥
cos θ
− 3 ⎟ senθ cos θ
sen θ + 3 cos θ ⎟
− 3 ⎟ senθ cos θ − ⎜
⎢ −⎜
⎜
2
sen 2θ + 3 cos 2 θ
⎥
l
l
l ⎠
l
l ⎠
⎝ l
⎝ l
⎠
l
l
⎢ ⎝ l
⎥
EI
4
2
EI
6 EI
6 EI
6 EI
6 EI
⎢
⎥
cos
θ
θ
−
θ
sen
cos
−
θ
sen
l
l
⎦
⎢⎣
l2
l2
l2
l2
[ ]