Prof. Luiz Antonio Farani de Souza
49 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Conteúdo Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos .................................. 49
4.1 Introdução ......................................................................................................................... 50
4.2.1 Teoria de Euler-Bernoulli............................................................................................ 51
4.3 Entendendo a estrutura .................................................................................................... 54
4.4 Matriz de rigidez global da estrutura ................................................................................ 55
4.5 Problema estrutural .......................................................................................................... 57
4.6 Metodologia ...................................................................................................................... 57
4.7 Programa de Elementos Finitos para viga e pórtico 2D .................................................... 58
4.7.1 Passos para criação do programa de elementos finitos ............................................ 58
4.7.2 Passos para a análise da estrutura (execução do programa). .................................... 65
4.8 Exercícios resolvidos de problemas de vigas .................................................................... 66
Exemplo 1 - Viga simplesmente apoiada: teoria de viga de Euler-Bernoulli ...................... 66
Exemplo 2 - Viga simplesmente apoiada: teoria de viga de Timoshenko ........................... 68
Exemplo 3 - Viga biengastada com carga distribuída.......................................................... 69
Exemplo 4 - Viga biapoiada com balanço e uma mola ........................................................ 71
4.9 Exercícios propostos de vigas ............................................................................................ 74
Exercício 1 .......................................................................................................................... 74
Exercício 2 .......................................................................................................................... 74
Prof. Luiz Antonio Farani de Souza
50 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
4.1 Introdução
As estruturas lineares são constituídas por uma ou mais barras (vigas, pilares, arcos,
pórticos, grelhas, etc.). As vigas são estruturas lineares, dispostas horizontalmente ou inclinadas,
com um ou mais apoios sujeitos a cargas transversais (Figura 4.1). A viga é geralmente utilizada
no sistema laje-viga-pilar (Figura 4.2) para transferir os esforços verticais recebidos da laje para
o pilar ou para transmitir uma carga estrutural concentrada, caso sirva de apoio a um pilar. Pode
ser composta de madeira, ferro ou concreto (português brasileiro) ou betão (português europeu)
armado. A viga transfere o peso das lajes e dos demais elementos (paredes, portas, etc.) às colunas
(SOUZA; RODRIGUES; MASCIA, 2008).
As edificações basicamente apresentam três tipos de vigas, que diferem na forma como
são ligados aos seus apoios. Portanto, classificam-se em:
- Viga em balanço: é uma viga de edificação com um só apoio. Toda a carga recebida é transmitida
a um único ponto de fixação (Figura 4.1a).
- Viga biapoiada: diz-se das vigas com dois apoios, que podem ser simples e/ou engastados,
gerando-se vigas do tipo simplesmente apoiadas, vigas com apoio simples e engaste, vigas
biengastadas (Figura 4.1b).
- Viga contínua: diz-se da viga com múltiplos apoios (Figura 4.1c).
Figura 4. 1 Vigas: a) em balanço; b) biapoiada; e c) contínua. Fonte: Souza, Rodrigues e Mascia (2008).
Figura 4. 2 Sistema laje-viga-pilar. Fonte:
https://eduqc.com.br/concursos/engenharia/engenharia-civil-analise-estrutural/
Pórticos planos são estruturas lineares planas com solicitações coplanares (Figura 4.3).
Podemos generalizar os pórticos planos para pórticos espaciais, cujas barras se dispõem em planos
diversos. Os pórticos, em conjunto com os elementos secundários, formam o esqueleto resistente
do sistema construtivo, onde são fixados os elementos de cobertura e de vedação lateral.
Prof. Luiz Antonio Farani de Souza
51 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
analise-estrutural/
Há duas teorias de vigas: Euler-Bernoulli e Timoshenko. A teoria clássica de Euler-
Bernoulli é considerada mais simples, porque os efeitos devidos às deformações cisalhantes na
seção transversal da barra são desprezados no cálculo da rigidez da estrutura. Nesse caso, estima-
se que a seção transversal plana antes da deformação da viga continuará sendo plana e
perpendicular ao eixo neutro da viga (SILVA et al., 2016).
A outra teoria de viga é a de Timoshenko, em que a deformação cisalhante na seção
transversal da barra é levada em consideração no processo de solução do problema. Considera-
se, entretanto, que a deformação por cisalhamento é constante ao longo da seção transversal.
Segundo Oñate (2013), essa teoria de viga é aplicável em vigas com esbeltezas menores (λ = L0/h
< 10; em que L0 é o comprimento da barra e h é a altura da seção transversal), e que ela se torna
especialmente importante à medida que a altura da seção transversal aumenta. Considera-se que
as seções transversais planas e perpendiculares ao eixo neutro se mantêm planas, mas não
perpendiculares ao eixo neutro deformado. A teoria de Timoshenko é considerada mais acurada
segundo Timoshenko et al. (1974).
4.2 Formulação de Elementos Finitos para viga-coluna 2D
4.2.1 Teoria de Euler-Bernoulli
A formulação de Elementos Finitos para a viga-coluna 2D com dois nós e três graus de
liberdade/nó é descrita a seguir. Designa-se por Euler-Bernoulli a formulação do elemento finito
de viga-coluna em que se considera que as seções se mantêm planas e normais ao eixo da barra
após a deformação. Deste modo não é considerada a deformação devida ao esforço cortante
(KWON; BANG, 2000; ZIENKIEWICZ; TAYLOR; ZHU, 2005; BATHE, 2016).
Na configuração inicial, as coordenadas dos nós 1 e 2 do elemento de viga no sistema
global de coordenadas são (X1, Y1) e (X2, Y2), respectivamente. O comprimento original
(indeformado) L0 da viga é dado pela seguinte equação:
L0 = √(X2 − X1) 2 + (Y2 − Y1)
2 (4.1)
Para o elemento de viga na sua configuração corrente, as coordenadas nodais globais são
(X1 + u1, Y1 + v1) para o nó 1 e (X2 + u2, Y2 + v2) para o nó 2, em que ui é o deslocamento do nó
i na direção X e vi é o deslocamento do nó i na direção Y, com i = 1, 2. O comprimento corrente
L (deformado) é:
L = √(X2 + u2 − X1 − u1) 2 + (Y2 + v2 − Y1 − v1)
2 (4.2)
Prof. Luiz Antonio Farani de Souza
52 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
O vetor de deslocamentos no sistema global de coordenadas p do elemento m é dado por:
m = [u1 v1 θ1 u2 v2 θ2] T (4.3)
Assumindo que a rigidez à flexão EI seja constante no elemento, a matriz de rigidez
elementar kelem no sistema local de coordenadas é dada por (KWON; BANG, 2000):
=
−62 2 0 0
2 =
0 (4.6)
na qual A é a área e I é o momento de inércia da seção transversal e é o módulo de elasticidade.
Nas estruturas em geral os elementos constituintes não possuem uma mesma inclinação
(vigas e pilares, por exemplo) o que faz com que o sistema local de um elemento não coincida
com o sistema local do outro, sendo então necessário reescrever as matrizes de rigidez dos
elementos em função de um único sistema de coordenadas, chamado de sistema global de
coordenadas. Isso é feito com auxílio de uma matriz conhecida por matriz de rotação R. A matriz
de rotação R é dada por:
=
(4.7)
em que C e S são o cosseno e o seno do ângulo dados por, respectivamente:
C = cos (θ) = X2 − X1
L0 (4.8)
L0 (4.9)
A matriz de rigidez no sistema global de coordenadas Kelem é dada por:
= T (4.10)
O vetor de deslocamentos no sistema local de coordenadas ulocal é calculado por:
Prof. Luiz Antonio Farani de Souza
53 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
θ2]
(4.11)
O vetor de forças internas elementar no sistema local de coordenadas felem é dado por:
= =
[ N1
V1
M1
N2
V2
M2]
(4.12)
na qual Ni, Vi e Mi, com i = 1,2, são a força normal, a força cortante e o momento fletor no nó i,
respectivamente.
No sistema global de coordenadas, o vetor de forças internas elementar Felem é
determinado por:
4.2.2 Teoria de viga de Timoshenko
Na formulação do elemento de viga de Timoshenko é considerado que as seções planas
se mantêm planas. No entanto, supõe-se que uma seção normal ao eixo da viga não mantém essa
característica após a deformação. Desse modo é possível considerar a deformação devida à força
cortante.
A formulação do Método do Elementos Finitos para a viga-coluna 2D para a teoria de
Timoshenko é semelhante à da viga de Euler-Bernoulli, com exceção à determinação da matriz
de rigidez elementar kelem no sistema local de coordenadas, que é determinada em função de duas
parcelas (KWON; BANG, 2000):
= +
Prof. Luiz Antonio Farani de Souza
54 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
na qual G é módulo de elasticidade transversal e é o fator de correção do cisalhamento, cujo
valor normalmente é 5/6.
4.3 Entendendo a estrutura
O modelo estrutural de viga ou pórtico deve conter informações acerca das dimensões e
apoios, o carregamento aplicado e as propriedades físicas e geométricas. Na Figura 4.4 é ilustrado
um exemplo de um pórtico plano em L com essas informações.
Figura 4. 4 Modelo estrutural do pórtico plano em L.
As dimensões são definidas com base no sistema de coordenadas cartesiano
bidimensional global (eixo de coordenadas axiais X e eixo de coordenadas transversais Y).
Os apoios são os admitidos os usuais (rolete, pino ou engaste), que podem possuir
restrição de deslocamento axial (u), de deslocamento transversal (v) e de rotação (), sendo
possível a combinações dos mesmos.
Os carregamentos na estrutura podem ser forças concentradas (axial e transversal), forças
distribuídas e momentos localizados.
As propriedades geométricas da seção transversal do elemento são a área A e o momento
de inércia em relação ao eixo centroidal I. As propriedades materiais são o módulo de elasticidade
longitudinal E, o módulo de elasticidade transversal G, o coeficiente de Poisson e o fator de
correção do cisalhamento .
Com a definição do modelo da estrutura, a mesma é discretizada por elementos finitos,
isto é, é dividida em partes menores que serão chamadas de elementos, que possuem nós em suas
extremidades que definem o seu início e o seu final. Na Figura 4.5 aparece o pórtico em L
discretizado em 16 elementos finitos de viga-coluna e 17 nós.
Figura 4. 5 Malha de elementos finitos do pórtico plano em L.
Prof. Luiz Antonio Farani de Souza
55 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
Cada elemento possui dois nós - um nó inicial (i) e um nó final (j). O elemento 1 da malha,
por exemplo, tem o nó inicial 1 e o nó final 2. Nota-se que, em geral, um mesmo nó pertence a
mais de um elemento (por exemplo, o nó 5 pertence aos elementos 4 e 5 da malha), sendo eles
responsáveis pela conectividade da estrutura.
Cada elemento apresenta 6 graus de liberdade no total (3 graus de liberdade por nó -
translações horizontal e vertical e uma rotação). E cada elemento possui o seu próprio sistema
local de coordenadas x, y que, usualmente, difere do sistema global de coordenadas em que a
estrutura é definida (Figura 4.6).
Figura 4. 6 Elemento finito no sistema local de coordenadas.
Salienta-se que apesar dos elementos possuírem um sistema local de coordenadas próprio
(x, y), todos os dados no modelo estrutural devem ser definidos a partir do sistema global de
coordenadas (X, Y), uma vez que este sistema vale para toda a estrutura e é com base nele que os
dados de entrada do programa são criados.
4.4 Matriz de rigidez global da estrutura
Para a definição da matriz de rigidez global K de uma estrutura do tipo viga ou pórtico
plano, faz-se o uso de um processo denominado redistribuição ordenada das matrizes locais, ou
seja, todas as matrizes de rigidez locais dos elementos k são combinadas de tal forma a representar
o modelo matricial da estrutura.
Esse processo leva em consideração a conectividade dos nós da estrutura o que gera no
sistema global uma superposição das matrizes de rigidez local para descrever o comportamento
físico de um ponto nodal. Considerando o sistema global de coordenadas X, Y, cada nó que une
duas barras compartilha 3 (três) graus de liberdade, sendo uma translação horizontal na direção
de X, uma translação vertical na direção de Y e uma rotação no nó. Conforme o tamanho da
estrutura, pode-se obter a dimensão da matriz de rigidez que dependerá apenas do número de
elementos da estrutura.
A matriz de rigidez do sistema global da estrutura é formada a partir das matrizes
elementares de rigidez no sistema local, resultante de um processo ordenado de soma dessas
matrizes. Durante esse processo, o programa realiza a soma das contribuições de cada elemento,
isto é, cada matriz de rigidez elementar é combinada em uma só matriz para que seja realizada a
análise da estrutura em relação ao seu sistema global de coordenadas. A seguir é apresentado de
forma simples como ocorre o processo de montagem da matriz de rigidez global da estrutura para
um pórtico biapoiado com 3 elementos e 4 nós, mostrado na Figura 4.7. Na Figura 4.8 aparecem
as matrizes elementares de rigidez k e na Figura 4.9 a matriz de rigidez global K do sistema
estrutural.
Prof. Luiz Antonio Farani de Souza
56 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
Figura 4. 7 Pórtico biapoiado com 3 elementos e 4 nós.
Figura 4. 8 Matrizes elementares de rigidez k.
A matriz K na Figura 4.9 é obtida com a combinação das matrizes elementares da Figura
4.8, utilizando-se o princípio da superposição de efeitos para representar numericamente o que
seria a ação de um nó na estrutura quando interliga duas ou mais barras. A matriz de rigidez de
uma estrutura descreve o comportamento matemático do modelo analisado em termos de flexão
e esforços normais. Por meio de manipulações realizadas diretamente na matriz global, é possível
inserir algumas condições de contorno ao problema proposto que representam a presença de
apoios (rolete, pino e engaste). Esse procedimento é conhecido como prescrição da matriz de
rigidez global e pode ser realizado por meio de uma operação matricial simples, em que, quando
uma condição de apoio já conhecida pode ser inserida antes do início do processamento da
estrutura.
De acordo com Soriano (2005), desde que sejam prescritos deslocamentos em número
suficiente para impedir os deslocamentos do corpo rígido da estrutura e esta não tenha
mecanismos internos, a matriz de rigidez global K fica não singular (det(K) 0), permitindo a
resolução do sistema de equações Ku = Fext, com a obtenção do vetor de deslocamentos nodais u
(vetor de incógnitas).
Prof. Luiz Antonio Farani de Souza
57 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
4.5 Problema estrutural
O sistema de equações de equilíbrio gerado pelo Método dos Elementos Finitos que
determina os deslocamentos nodais u de uma viga ou pórtico com comportamento elástico linear
e supondo deformações infinitesimais é (ZIENKIEWICZ; TAYLOR; ZHU, 2005; BATHE,
2016):
= (4.17)
na qual K é a matriz de rigidez global da estrutura e Fext é o vetor de forças externas.
A matriz de rigidez K é obtida por:
= ∑
(4.18)
em que e é o elemento finito e n é o número total de elementos da estrutura.
4.6 Metodologia
Utilizando o programa livre Scilab versão 6.1.0 (SCILAB, 2020), realiza-se a
implementação do elemento de viga-coluna 2D para análise linear estática de vigas e pórticos
planos.
Inicialmente, o programa faz uma leitura dos dados de entrada (função
entrada_dados.sci). Assim, o usuário deve fornecer ao programa todas as informações
necessárias acerca da estrutura a ser analisada.
Em seguida, fazendo-se uso dos dados de entrada, o programa (principal.sce) processa-
os e calcula as matrizes de rigidez de cada elemento k. Essas matrizes são então usadas para
montar a matriz de rigidez global do sistema estrutural K.
O próximo passo é a aplicação das condições de contorno do problema. Essas condições
são pré-estabelecidas nos dados de entrada (vetor NOCC). Dessa forma, tem-se um sistema de
equações lineares (Ku = Fext) com uma matriz de rigidez não singular. Com a solução desse
sistema, obtém-se os valores dos deslocamentos nodais u. Posteriormente, o programa calcula a
os esforços internos (momento fletor M, força cortante V e força normal N) nos nós da malha e
as reações nos apoios. A Figura 4.10 apresenta o fluxograma de montagem do código
computacional.
Prof. Luiz Antonio Farani de Souza
58 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
4.7 Programa de Elementos Finitos para viga e pórtico 2D
4.7.1 Passos para criação do programa de elementos finitos
1. Fazer o download do programa Scilab, disponível no endereço
https://www.scilab.org/download/6.1.0.
2. Criar uma pasta com o nome do programa (por exemplo, "portico2D").
3. Dentro dessa pasta, criar no Scilab os arquivos referentes ao programa principal
(principal.sce) e às funções, utilizando o SciNotes.
Importante:
- As funções deverão ser salvas com a extensão .sci e o programa principal com a extensão
.sce; e
- Salvar os arquivos com o mesmo nome dentro da pasta "portico2D".
• programa principal.sce
//Teorias de viga de Euler-Bernoulli e Timoshenko
//____________________________________________
Prof. Luiz Antonio Farani de Souza
59 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
[K]=DK(NOCC,NNOSCC,NTGL,NTEL,dofno,inci,coord,E,A,I,itipo,G)
[Fint,Reacao,Ef]=Dfint
for i=1:NTNOS
end
n=length(NOCC)
for m=1:NTEL+1
M(m)=Ef(3*m);
V(m)=Ef(3*m-1);
N(m)=Ef(3*m-2);
for i=1:NTNOS
Prof. Luiz Antonio Farani de Souza
60 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
set(gca(),"auto_clear","off");
set(gca(),"auto_clear","off");
end
xlabel('x (m)','fontsize',2); //Eixo x
ylabel('y (m)','fontsize',2); //Eixo y
legend('posição inicial','posicão final',4); //Legenda
• função entrada_dados.sci
function [coord, inci, NTEL, NTNOS, NTGL, NNOSCC, dofno, E, A, I, NOCC, Fext, itipo, G,
fator, Fext2, q0]=entrada_dados()
//Entrada de dados (pré-processamento)
//NTGL -> NÚMERO TOTAL DE GRAUS DE LIBERDADE
//NNOSCC -> NÚMERO DE GRAUS RESTRITOS (CONDIÇÕES DE CONTORNO)
NTNOS=41;
NTEL=40;
//material 1 -> itipo(nel,3)==1
//material 2 -> itipo(nel,3)==2
itipo(i,3) = 1; //material 1
for i=1:NTEL
Prof. Luiz Antonio Farani de Souza
61 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
dofno(i,2)=inci(i,2)*3-1;
dofno(i,3)=inci(i,2)*3;
dofno(i,5)=inci(i,3)*3-1;
dofno(i,6)=inci(i,3)*3;
for m=1:NTEL
//%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
Fext=zeros(NTGL,1);
Fext(21*3-1,1)=-1000; //inserir o valor força concentrada
q0=-200; //inserir o valor da força distribuída (força/comprimento)
//cria o vetor de forças distribuídas
Fdistr=zeros(NTGL,1);
[Fdistr]=fdistr(NTEL,inci,q0,NTGL,dofno,itipo);
Fext2=Fext; //vetor utilizado para o cálculo das reações
[Fext]=contfg(NOCC,NNOSCC,Fext);
//determina o vetor de forças distribuídas
Fdistr=zeros(NTGL,1);
C=(X2-X1)/L0(m); //cosseno de teta
S=(Y2-Y1)/L0(m); //seno de teta
R=[C S 0 0 0 0; //matriz de transformação
-S C 0 0 0 0;
0 0 1 0 0 0;
0 0 0 C S 0;
0 0 0 -S C 0;
0 0 0 0 0 1];
Prof. Luiz Antonio Farani de Souza
62 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
[Fdistr]=ensamfg(m,fdistr,dofno,itipo,Fdistr);
//graus de liberdade nos nós do elemento
//ELEMENTO VIGA 2 NÓS COM 3GL/NÓ
if (itipo(m,2)==2) || (itipo(m,2)==3)
IPO(1)=1;
IPO(2)=2;
IPO(3)=3;
IPO(4)=4;
IPO(5)=5;
IPO(6)=6;
TAM=6;
function [K]=contkg(NOCC, NNOSCC, NTGL, K)
//impõe as condições de contorno na matriz de rigidez global K
for J=1:NNOSCC
for I=1:NTGL
• função Dfint.sci
function [Fint, Reacao, Ef]=Dfint (u, NOCC, NNOSCC, NTGL, NTEL, dofno, inci, coord, E, A, I,
itipo, G, fator, Fext2, q0)
//Determina o vetor global de forças internas
//Determina as reações nos apoios e os esforços internos (N,V,M)
Fint=zeros(NTGL,1);
Fint2=zeros(NTGL,1);
pm(i,1)=u(dofno(m,i),1); //seleciona os deslocamentos nos nós do elemento m
end
L0(m)=sqrt((X2-X1)^2+(Y2-Y1)^2); //comprimento indeformado L0 da viga
L(m)=sqrt((X2+u(4,1)-X1-u(1,1))^2 + (Y2+u(5,1)-Y1-u(2,1))^2); //comprimento deformado L da
viga
R=[C S 0 0 0 0; //matriz de transformação
-S C 0 0 0 0;
0 0 1 0 0 0;
0 0 0 C S 0;
0 0 0 -S C 0;
Prof. Luiz Antonio Farani de Souza
63 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
u_local=R*pm;
//matriz de rigidez elementar no sistema local de coordenadas
c1=I(m)/L0(m)^2;
0 12*c1 6*c2 0 -12*c1 6*c2;
0 6*c2 4*I(m) 0 -6*c2 2*I(m);
-A(m) 0 0 A(m) 0 0;
0 -12*c1 -6*c2 0 12*c1 -6*c2;
0 6*c2 2*I(m) 0 -6*c2 4*I(m)];
end
//matriz de rigidez elementar no sistema local de coordenadas
kelemb=E(m)/L0(m)*[ A(m) 0 0 -A(m) 0 0;
0 0 0 0 0 0;
0 0 I(m) 0 0 -I(m);
-A(m) 0 0 A(m) 0 0;
0 0 0 0 0 0;
0 0 -I(m) 0 0 I(m)];
kelems=G(m)*fator(m)*A(m)/(4*L0(m))*[0 0 0 0 0 0;
0 4 2*L0(m) 0 -4 2*L0(m);
0 2*L0(m) L0(m)^2 0 -2*L0(m) L0(m)^2
0 0 0 0 0 0;
0 -4 -2*L0(m) 0 4 -2*L0(m);
0 2*L0(m) L0(m)^2 0 -2*L0(m) L0(m)^2];
kelem = kelemb + kelems;
felem=kelem*u_local; //vetor forças internas elementar no sistema local
Ef(3*m-2)=-felem(1,1); //força normal N no nó i do elemento m
Ef(3*m-1)=felem(2,1); //força cortante V no nó i do elemento m
Ef(3*m)=-felem(3,1); //momento fletor M no nó i do elemento m
FELEM=R'*(felem); //vetor forças internas elementar no sistema global
[Fint]=ensamfg(m,FELEM,dofno,itipo,Fint); //vetor forças internas global
[Fint2]=ensamfg(m,fdistr,dofno,itipo,Fint2); //vetor força distribuída global
end
Ef(3*(m+1)-2)=felem(4,1); //força normal N no nó j do elemento m+1
Ef(3*(m+1)-1)=-felem(5,1); //força cortante V no nó j do elemento m+1
Ef(3*(m+1))=felem(6,1); //momento fletor M no nó j do elemento m+1
//Reações nos graus de liberdade (apoios)
n=length(NOCC)
for m=1:NTEL
end
Ef(3*(m+1))=(Ef(3*(m+1))-Fint2(3*(m+1)));
[Fint]=contfg(NOCC,NNOSCC,Fint); //impõe as condições de contorno
endfunction
• função DK.sci
function [K]=DK(NOCC, NNOSCC, NTGL, NTEL, dofno, inci, coord, E, A, I, itipo, G, fator)
Prof. Luiz Antonio Farani de Souza
64 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
K=zeros(NTGL,NTGL);
C=(X2-X1)/L0(m); //cosseno de teta
S=(Y2-Y1)/L0(m); //seno de teta
R=[C S 0 0 0 0; //matriz de transformação
-S C 0 0 0 0;
0 0 1 0 0 0;
0 0 0 C S 0;
0 0 0 -S C 0;
0 0 0 0 0 1];
KELEM=zeros(6,6);
//matriz de rigidez elementar no sistema local de coordenadas
c1=I(m)/L0(m)^2;
0 12*c1 6*c2 0 -12*c1 6*c2;
0 6*c2 4*I(m) 0 -6*c2 2*I(m);
-A(m) 0 0 A(m) 0 0;
0 -12*c1 -6*c2 0 12*c1 -6*c2;
0 6*c2 2*I(m) 0 -6*c2 4*I(m)];
end
//matriz de rigidez elementar no sistema local de coordenadas
kelemb=E(m)/L0(m)*[ A(m) 0 0 -A(m) 0 0;
0 0 0 0 0 0;
0 0 I(m) 0 0 -I(m);
-A(m) 0 0 A(m) 0 0;
0 0 0 0 0 0;
0 0 -I(m) 0 0 I(m)];
kelems=G(m)*fator(m)*A(m)/(4*L0(m))*[0 0 0 0 0 0;
0 4 2*L0(m) 0 -4 2*L0(m);
0 2*L0(m) L0(m)^2 0 -2*L0(m) L0(m)^2
0 0 0 0 0 0;
0 -4 -2*L0(m) 0 4 -2*L0(m);
0 2*L0(m) L0(m)^2 0 -2*L0(m) L0(m)^2];
kelem = kelemb + kelems;
end
KELEM=R'*kelem*R; //matriz de rigidez elementar no sistema global de coordenadas
[K]=ensamkg(m,KELEM,dofno,itipo,K); //monta a matriz de rigidez global K do sistema
end
endfunction
[IPO,TAM]=apontador(m,itipo);
Prof. Luiz Antonio Farani de Souza
65 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
K(P,Q)=K(P,Q)+KELEM(I,J);
//monta o vetor de força interna global Fint
[IPO,TAM]=apontador(m,itipo);
function [Fint]=contfg(NOCC, NNOSCC, Fint)
//impõe as condições de contorno no vetor de força interna global Fint
for I=1:NNOSCC
4.7.2 Passos para a análise da estrutura (execução do programa).
1. Executar o programa Scilab .
2. Selecionar a pasta que contém os arquivos do programa de EF criado (botão "Selecione
um diretório").
3. Inserir os dados de entrada da estrutura (malha de elementos finitos, propriedades
geométricas e materiais das barras e vetor de força externa) na função entrada_dados.sci. As
seguintes informações devem ser colocadas nessa função:
a) Informar:
NNOSCC → número de graus de liberdade restringidos (condições de contorno)
b) posição dos nós no sistema global de coordenadas X, Y (matriz coord);
c) conectividades dos nós com as barras (incidência dos elementos) (matriz inci);
d) tipo de elemento (escolher a teoria de viga a ser adotada - Euler-Bernoulli ou
Timoshenko e o tipo de material) (matriz itipo);
Prof. Luiz Antonio Farani de Souza
66 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
e) propriedades geométricas e materiais de cada barra: área A e momento de inércia I da
seção transversal, o módulo de elasticidade longitudinal E, módulo de elasticidade transversal G
e fator de cisalhamento ;
f) vetor de força externa (Fext) - indicar as forças concentradas e momentos nos graus de
liberdade e a força distribuída q0; e
g) informar os graus de liberdade restringidos (apoios) (vetor NOCC).
4. Executar o programa.
5. Os resultados numéricos (deslocamentos nos nós, esforços internos (N, V e M) nas
barras e reações nos apoios) são imprimidos no console e a configuração deformada da estrutura
é apresentada numa figura.
4.8 Exercícios resolvidos de problemas de vigas
Nesta seção são apresentados quatro exemplos numéricos de vigas, em que são realizadas
análises estáticas lineares (material no regime elástico linear e hipótese de deformações
infinitesimais) por meio do programa de Elementos Finitos desenvolvido no programa Scilab
versão 6.1.0 (ver seção 4.7).
Exemplo 1 - Viga simplesmente apoiada: teoria de viga de Euler-Bernoulli
A Figura 4.11 mostra uma viga simplesmente apoiada com comprimento 20 in e módulo
de elasticidade E = 1,0 × 107 lb/in2. A seção transversal é quadrada com dimensões (1,0 x 1,0)
in2. A viga está sujeita a uma força vertical de 100 lb aplicada no seu centro. A malha de
Elementos Finitos consiste de 11 nós e 10 elementos de viga. Esta viga foi estudada por Kwon e
Bang (2000). Foi utilizada na simulação a teoria de viga de Euler-Bernoulli.
Nó X Y Elemento i j
1 0 0 1 1 2
2 2 0 2 2 3
3 4 0 3 3 4
4 6 0 4 4 5
5 8 0 5 5 6
6 10 0 6 6 7
Prof. Luiz Antonio Farani de Souza
67 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
7 12 0 7 7 8
8 14 0 8 8 9
9 16 0 9 9 10
10 18 0 10 10 11
11 20 0 Figura 4. 11 Viga simplesmente apoiada sujeita a uma força. Fonte: Kwon e Bang (2000).
• Dados de entrada (função entrada_dados.sci):
Dados da malha Coordenadas
restringidos
Solução:
A configuração deformada da viga é mostrada na Figura 4.12. Na Tabela 4.1 são
apresentados os deslocamentos e rotações nodais obtidos com o programa e por Kwon e Bang
(2000). As reações nos apoios aparecem na Tabela 4.2.
Tabela 4. 1 Resultados numéricos (deslocamentos e rotações nos nós) obtidos com o programa e por Kwon e Bang (2000).
Programa Kwon e Bang (2000)
Nó u (in) v (in) (rad) Nó u (in) v (in) (rad)
1 0 0 0,003 1 0 0 0,003
2 0 0,00592 0,00288 2 0 0,0059 0,0029
3 0 0,01136 0,00252 3 0 0,0113 0,0025
4 0 0,01584 0,00192 4 0 0,0158 0,0019
5 0 0,01888 0,00108 5 0 0,0188 0,0011
6 0 0,02 -1,154963532D-17 6 0 0,02 0
Tabela 4. 2 Resultados numéricos (reações nos apoios) obtidos com o programa.
Grau de liberdade Reação
Prof. Luiz Antonio Farani de Souza
68 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
Forças:
Momento:
Exemplo 2 - Viga simplesmente apoiada: teoria de viga de Timoshenko
Neste exemplo é realizada análises com a teoria de Timoshenko para a viga do Exemplo
1, variando-se somente a altura h da seção transversal (h = 1,0 in, 1,5 in e 2,0 in). Considera-se o
módulo transversal de cisalhamento G = 3,8 106 lb/in2. A solução numérica é comparada com
a solução analítica (ver tabela no Anexo 4.1) e com a teoria de Euler-Bernoulli. Na entrada de
dados, são alterados o tipo de elemento e as propriedades geométricas A e I da seção transversal
da viga. Foram considerados 40 elementos de viga e 41 nós para a malha.
• Dados de entrada (função entrada_dados.sci):
Dados da malha Coordenada
Fext=zeros(NTGL,1);
Prof. Luiz Antonio Farani de Souza
69 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
Solução:
A configuração deformada da viga para h = 2 in, considerando a teoria de Timoshenko, é
mostrada na Figura 4.13. Na Tabela 4.3 são apresentados os deslocamentos verticais no centro da
viga, variando-se a altura h da seção, obtidas com o programa considerando as teorias de viga de
Euler-Bernoulli e Timoshenko e a solução analítica (viga biapoiada espessa).
Tabela 4. 3 Deslocamentos verticais no centro da viga, variando-se a altura h da seção, obtidas com o programa considerando as teorias de viga de Euler-Bernoulli e Timoshenko, e a solução analítica.
h A (in2) I (in4) Programa Solução analítica
(viga espessa) Euler-Bernoulli Timoshenko
1 1 0,083333333 0,02 0,020145395 0,020157895
1,5 1,5 0,28125 0,005925926 0,006027485 0,006031189
2 2 0,666666667 0,0025 0,002577385 0,002578947
A solução analítica para o deslocamento máximo vertical para viga espessa é dada pela
expressão:
Figura 4. 13 Configuração deformada da viga.
Exemplo 3 - Viga biengastada com carga distribuída
Seja a viga biengastada sujeita ao carregamento mostrado na Figura 4.14. O módulo de
elasticidade da viga é E = 1,18 108 kN/m2. A seção transversal é em T, cuja área é A = 0,09 m2
e o momento de inércia I = 1,9639 10-3 m4. A malha de elementos finitos consiste de 41 nós e
40 elementos de viga. Utiliza-se a teoria de viga de Euler-Bernoulli.
Prof. Luiz Antonio Farani de Souza
70 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
Figura 4. 14 Modelo estrutural da viga biengastada com carga distribuída.
• Dados de entrada (função entrada_dados.sci):
Dados da malha Coordenada
//inserir a força concentrada
q0=-200;
NOCC=[1 2 3 3*41-2 3*41-1 3*41]
Solução:
Na Tabela 4.4 são apresentados os deslocamentos e a rotação no centro da viga obtidos
com o programa desenvolvido e o Ftool. As reações nos engastes (nós 1 e 41) são apresentadas
na Tabela 4.5. O diagrama de momento fletor aparece na Figura 4.15 e a configuração deformada
da viga é mostrada na Figura 4.16.
Tabela 4. 4 Deslocamentos no nó central da viga (nó 21).
Programa Ftool
Nó u (m) v (m) (rad) Nó u (m) v (m) (rad)
21 0,0000094161959 0,0002157588541 0 21 9,416e-06 -2,158e-04 0
Tabela 4. 5 Reações nos engastes (nós 1 e 41).
Nó
1
3 (Mz) 316.6666666638 kN m 316.6666666638 kN m
41
123 (Mz) -316.6666666638 kN m -316.6666666638 kN m
Prof. Luiz Antonio Farani de Souza
71 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
Figura 4. 15 Diagrama de momento fletor.
Figura 4. 16 Configuração deformada da viga.
Código para a geração do diagrama de momento fletor (programa principal em "Saída de Dados"): y=zeros(NTNOS,1)
x=[coord(:,1)];
Exemplo 4 - Viga biapoiada com balanço e uma mola
Seja a viga biapoiada com um balanço e uma mola (rigidez da mola kmola = 200 MN/m),
sujeita ao carregamento mostrado na Figura 4.17. O módulo de elasticidade da viga é E = 1,18
108 kN/m2. A seção transversal é retangular com dimensões (0,1 0,3) m2. A malha de elementos
finitos consiste de 41 nós e 40 elementos de viga. Utiliza-se a teoria de viga de Euler-Bernoulli.
Prof. Luiz Antonio Farani de Souza
72 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
Dados da
Fext(41*3-1,1)=-1000; //inserir o
valor da força concentrada
Fext(41*3,1)=-500; //inserir o
distribuída (força/comprimento)
Solução:
Na Tabela 4.6 são apresentados os deslocamentos e as rotações nos nós 1, 21 e 41 obtidos
com o programa desenvolvido e o Ftool. O diagrama de momento fletor aparece na Figura 4.17 e
a configuração deformada da viga é mostrada na Figura 4.18.
//______________________________
[K]=DK(NOCC,NNOSCC,NTGL,NTEL,dofno,inci,coord,E,A,I,itipo,G)
//K(122,122)=K(122,122)+kmola
K(41*3-1,41*3-1) = K(41*3-1,41*3-1) + kmola; //soma a rigidez da mola no grau de liberdade
vertical do nó 41 (41*3-1=122)
//Deslocamentos nodais u
[Fint,Reacao,Ef]=Dfint2(u,NOCC,NNOSCC,NTGL,NTEL,dofno,inci,coord,E,A,I,itipo,G,fator,Fext2,
q0)
Tabela 4. 6 Deslocamentos e rotações nos nós 1, 21 e 41 da viga.
Programa Ftool
u
21 0 -0,0013039203079 21 0 0 -1,304e-03
41 0 -0,0029674387059 -0,0043385950402 41 0 -2,967e-03 -4,339e-03
Prof. Luiz Antonio Farani de Souza
73 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
Cálculo da força na mola:
--> Fmola=2*10^5*(0.0029674387059)
Fmola =
593.48774118000
Prof. Luiz Antonio Farani de Souza
74 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
4.9 Exercícios propostos de vigas
Resolver os exercícios de vigas a seguir utilizando o código de elementos finitos
desenvolvido no programa Scilab (ver seção 4.7).
Exercício 1 A viga em T de madeira é feita de duas tábuas de 200 mm x 30 mm. Considerar o
módulo de elasticidade da madeira E = 13,1 GPa. Este exercício foi adaptado de Hibbeler (2004).
Desprezar o peso próprio. Considerar a teoria de Euler-Benoulli. Considerar a malha com 10
elementos de viga. Determinar:
b) As rotações (inclinações) nos apoios;
c) Se a tensão de flexão admissível for adm = 12 MPa e a tensão de cisalhamento for adm = 0,8
MPa, a viga suportará com segurança o carregamento mostrado? Justificar calculando as tensões
normal (x) e de cisalhamento () máximas na viga.
Exercício 2 Seja uma viga biengastada sujeita à força distribuída q = 500 kN/m, comprimento
indeformado L0 = 1,0 m e de seção transversal retangular com base b = 0,1 m e altura h. Considere
o material concreto (E = 22,1 GPa e o coeficiente de Poisson = 0,15). O módulo de elasticidade
transversal G é obtido pela expressão:
G =
2(1 + )
(4.20)
Calcular o deslocamento máximo na viga com as teorias de viga de Euler-Bernoulli e
Timoshenko, variando-se a altura da seção h = 0,1 m, 0,5 m, 1,0 m e 1,5 m, e comparar com as
soluções analíticas. Também, obter os deslocamentos com diferentes malhas de elementos finitos:
4, 8, 16 e 32 elementos de viga, e determinar o erro relativo em porcentagem:
erro (%) = | desl. aproximado − desl. analítico
desl. analítico |
para as duas teorias de viga.
Prof. Luiz Antonio Farani de Souza
75 Capítulo 4 - Método dos Elementos Finitos para vigas e pórticos planos
Programação no Scilab
Referências
BATHE, K. J. Finite element procedures. 2ª Ed. Watertown, MA: Klaus-Jurgen Bathe, 2016.
FTOOL, version 4.00.04. Pontifícia Universidade Católica, Tecgraf/PUC-Rio, Rio de Janeiro,
2018.
HIBBELER, R. C. Resistência dos Materiais. São Paulo: Pearson Education do Brasil, 2004.
KWON, Y. W.; BANG, H. The Finite Element Method using MATLAB. Boca Raton, FL:
CRC 672 Press, 2000.
OÑATE, E. Structural Analysis with the Finite Element Method. Vol. 1, International Center
for Numerical Methods in Engineering (CIMNE), Spain, 2013.
SCILAB, version 6.1.0. ESI Group, 2020.
SILVA, J. L.; LEMES, Í. J. M.; SILVEIRA, R. A. M.; SILVA, A. R. D. Influência da teoria de
viga na análise geometricamente não linear de estruturas reticuladas. In: XXXVII Iberian Latin
American Congress on Computational Methods in Engineering, Brasília, Brasil, 2016. Anais...
Brasília: XXXVII CILAMCE, 2016.
SORIANO, H. L. Análise de estruturas - Formulação matricial e implementação computacional.
Rio de Janeiro: Editora Ciência Moderna, 2005.
SOUZA, M. F. S. M.; RODRIGUES, R. B.; MASCIA, N. T. Sistemas estruturais de edificações
e exemplos. Departamento de Estruturas – DES, FEC - UNICAMP, Campinas, 2008.
TIMOSHENKO, S.; YOUNG, D. H.; WEAVER, J. W. Vibration Problems in Engineering. 1a
Ed. Rio de Janeiro: Guanabara Dois, 1974.
ZIENKIEWICZ, O. C.; TAYLOR, R. L.; ZHU, J. Z. The finite element method: its basis and
fundamentals. Elsevier, 2005.