Trabalho 1 Portainer

Prévia do material em texto

APLICAÇÃO DO MÉTODO DE ELEMENTOS FINITOS A UM GUINDASTE PÓRTICO
B. L. R. Hoffmann, G. S. de M. Martins, J. F. Kopetski, M. de O. B. Costa, V. H. Ilgenfritz, W. G. de Andrada 
Universidade Federal de Santa Catarina; Engenharia Naval; Disciplina de Métodos Computacionais para Engenharia; Professor Vitor Takashi Endo
Resumo
Este trabalho preconiza a aplicação do Método de Elementos finitos a um guindaste pórtico, cuja geometria é dada na forma de uma treliça. A metodologia implementada tem como enfoque a simulação computacional bidimensional, por meio de dois programas distintos, ABAQUSTM e SciLabTM. Portanto, o objetivo consiste na análise de ambos os resultados e na ponderação e modificação dos parâmetros, com o propósito de obter benefícios relevantes no projeto estrutural.
Keywords: Treliça, Elementos Finitos, Computacional, Análise, Resistência estrutural.
�
Introdução
O Método de Elementos Finitos (MEF) aplicado ao âmbito de engenharia de estruturas tem por objetivo, determinar o estado de tensão e deformação de uma estrutura submetida a solicitações externas. A representação da estrutura contínua é simplificada a partir do processo de discretização, ou seja, a mesma é repartida em inúmeros elementos finitos de análise [1].
Em certos estudos, a complexidade do problema caracterizada por aspectos como: o comportamento do material; da geometria; das condições de fronteira ou até mesmo o número de elementos; inviabilizam a solução analítica. Em tais situações recorre-se a utilização de recursos computacionais para possibilitar a resolução [2][3].
A competência de um engenheiro de estruturas abrange desde a introdução de dados no programa computacional até a interpretação dos resultados. Sobretudo, é oportuno verificar diferentes situações de análise para criar uma concepção crítica e ponderar os resultados oriundos do meio computacional [2].
Sendo assim, este trabalho propõe a aplicação bidimensional do MEF a um guindaste pórtico, cuja a geometria é dada na forma de treliça, com o intuito de analisar os resultados obtidos (como a distribuição de tensão e deslocamentos nodais) a partir de diferentes programas computacionais: ABAQUSTM e SciLabTM. Além disso, será realizado uma análise crítica através da modificação da área de secção com o propósito de otimizar o projeto estrutural.
 
Metodologia
Descrição da estrutura de análise
A estrutura escolhida para análise consiste num guindaste pórtico, equipamento utilizado em atividades portuárias para a movimentação de cargas, em geral, os contêineres [4]. Os critérios utilizados na escolha do objeto de estudo levaram em consideração a empregabilidade no âmbito naval e também, na caracterização da geometria por meio de treliças. A figura 1 destaca a estrutura real de análise. 
Fig. 1. Guindaste pórtico [4]
Definição dos parâmetros de análise
Em ambas as análises computacionais é necessário atribuir certos parâmetros durante a simulação e principalmente, garantir a atribuição idêntica dos dados para não interferir nos resultados. 
Dentre estes parâmetros, a solicitação de carga é um aspecto evidente. Segundo o fabricante Liebherr [4], a estrutura é dimensionada para suportar uma carga de 400000 N, cuja distribuição do carregamento admitido pelos projetistas foi de carga concentrada, visto que esta representa a pior condição de análise.
Outras considerações atribuídas foram as caraterísticas do material de fabricação da estrutura adotado como aço A36, cuja a uma tensão de escoamento; a densidade; o módulo de elasticidade e o coeficiente de Poisson são respectivamente, 250 MPa; 7,8E-12 ton/mm³; 200 MPa; e 0,3.
A descrição da geometria, bem como as dimensões empregadas foram obtidas através de dados fornecidos pelo fabricante e normas comercias. Somente a discretização dos elementos, definida em 77 elementos, tal como suas respectivas medidas foram decisões tomadas pelos projetistas. Para realizar a análise de tensão e deslocamentos, considerando os critérios de falha, foram utilizadas as seguintes áreas de seção transversal, conforme a tabela 1.
Tabela 1
Descrição das áreas dos elementos (mm²)
	Elementos
	Área 
	Treliça 1 (Superiores)
	1500 mm²
	Treliça 2 (Inferiores)
	20000 mm²
	Treliça 3 (Verticais de apoio)
	40000 mm²
Tais valores foram dimensionados com base em estruturas reais e levando em consideração a norma da ABNT NBR 8800/2008 a qual define máximos deslocamentos verticais para estruturas feitas em aço [5]. De acordo com o Anexo A, referente à norma, o deslocamento máximo para este tipo de estrutura é definido pela seguinte equação.
δmáx = L/800 (1)
Onde L é o vão livre definido em 20 metros. Além disso, nestas simulações foram utilizadas duas condições de contorno, sendo engastes nos pontos em vermelho, e pinos nos pontos em verde, conforme demonstrado na figura 2.
Fig. 2. Geometria CAE e dimensões da estrutura em mm
Resolução no ABAQUSTM
Primeiramente, para a descrição computacional da geometria foram testados dois modelos de treliça, Pratt e Warren sob as mesmas condições (Apêndice A). O objetivo deste teste prevê a seleção da melhor modelagem para análise, ou seja, aquela que fornece maior resistência estrutural.
A treliça Pratt é composta por elementos diagonais que apontam para o centro inferior do vão. Exceto pelos elementos centrais, os demais elementos diagonais estão sujeitos somente à tração, enquanto os elementos verticais suportam as forças de compressão [6]. Já a tipo Warren é uma estrutura simples e contínua caracterizada por parte das diagonais estarem sujeitas à tração e a outra à compressão. É utilizada em vãos menores, quando não há necessidade de utilizar elementos verticais como amarra [6][7].
Contudo, ao final dos ensaios, o modelo escolhido para prosseguir a análise foi do tipo Warren, já que os valores de tensão atuantes foram menores. Além disso, outro critério adotado na seleção da modelagem se deu àquela que tivesse um número menor de elementos submetidos à compressão expressiva. Este fato facilita o estudo do critério de falha por flambagem, já que o tipo Warren possui somente dois elementos com valores significativos submetidos à compressão. Portanto, na análise de falha, caso os mesmos não sofram o efeito da flambagem é possível admitir que toda a estrutura estará segura. 
A partir da definição do tipo de treliça e da implementação das condições de contorno e parâmetros de análise discutidos anteriormente, é possível simular o comportamento da estrutura, em termos de deslocamento nodal e tensão, demonstrados nas figuras 3 e 4.
Fig. 3. Magnitude dos deslocamentos em mm
Fig. 4. Magnitude das tensões em MPa
Contudo, é necessário realizar uma análise de falha por flambagem. Sendo assim, foi realizado através do Excel, o cálculo analítico da tensão crítica para o efeito de flambagem, conforme a relatado na equação seguinte. 
σcr = π²EI/ALc² (2)
Onde Lc é o comprimento crítico, adotado como sendo igual ao comprimento das barras (1m). Ao observar as figuras 3 e 4, nota-se que somente os elementos superiores (Treliça 1), são comprimidos, com tensões próximas a 120 MPa. Portanto, não ocorrerá falha por flambagem, visto que a tensão crítica calculada anteriormente é de 2356,2 MPa.
Ainda foi realizada uma análise de carregamento unidimensional na direção x para observar a influência da mesma no estado de tensão da estrutura. Aplicou-se a mesma magnitude de carga correspondente a 400000 N na direção x. Conforme a simulação computacional disposta nas figuras 5 e 6, percebe-se que a força aplicada na direção x somente influencia a tensão de tração submetida à estrutura, para o caso de análise unidimensional.
Fig. 5. Magnitude dos deslocamentos em mm
Fig. 6. Magnitude da tensão em Mpa
Ainda julgou-se necessário realizar uma análise para a condição de carregamento distribuídona estrutura, uma vez que na maioria das aplicações esta apresenta um acoplamento através do qual é conectada a carga, conforme explicitado na figura 7.
Fig. 7. Destaque do acoplamento da carga na estrutura
Atribui-se a largura de carga corresponde ao comprimento do carro de deslocamento do contêiner, que compreende valores de 2 a 5 m de comprimento. Por fim, repetiu-se a análise de tensões e deslocamentos, utilizando 2m de largura de carga. No qual os resultados da simulação estão dispostos nas figuras 8 e 9. 
Fig. 8. Magnitude dos deslocamentos em mm
Fig. 9. Magnitude das tensões em MPa
Resolução no SciLabTM
Com o intuito de verificar e confirmar o resultado obtido no software ABAQUSTM fez-se uma simulação via software SciLabTM. Os dados de entrada considerados no programa são oriundos da simulação via ABAQUSTM, como a posição dos nós e os nós associados a cada elemento. Além disso, também foi ponderado o módulo de elasticidade do material e a área da seção transversal de cada barra. 
Para implementar a condição de contorno como pinada, restringiu-se o movimento na direção x e y do nó 2 (GL 3 e 4) e do nó 40 (GL 79 e 80). A aplicação da força se deu no nó 19 na direção y (GL 38) no sentido negativo. Contudo, é possível aderir informações ao código computacional disponibilizado pelo professor Vitor Endo e executar o programa. Por fim, são obtidos os valores de tensão e deslocamento correspondente a cada um dos nós. O programa em SciLab pode ser observado no Apêndice B.
Análise dos Resultados
Ao comparar os resultados oriundos de ambos os programas computacionais, percebe-se que os resultados estão bem próximos conforme esperado, uma vez que a simulação no SciLabTM utiliza o artifício advindo do ABAQUSTM. A única dificuldade surge na implementação correta dos dados de entrada e condições de análise. A tabela 2 demostra o comparativo entre os resultados.
Tabela 2
Comparativo dos resultados computacionais
	
	Tensão (MPa)
	Deslocamentos (mm)
	
	Tração
	Compressão
	
	ABAQUSTM
	35,56
	-120
	21,47
	SciLabTM
	35,575258 
	-120
	21,478973
Otimização
O processo de otimização é um enfoque recorrente em projetos de engenharia, visto que quando possível, a redução do peso estrutural agrega benefícios relevantes, sobretudo, relacionados ao setor financeiro.
Logo, conforme a proposta enunciada neste trabalho, após a análise da estrutura com as devidas dimensões comerciais, é possível otimizar a estrutura de modo a reduzir o peso estrutural e por sua vez, o uso de material.
Sendo assim, uma nova simulação computacional no ABAQUSTM foi desenvolvida a partir da redução das áreas de seção transversal dos elementos, levando em consideração apenas a tensão de escoamento do material. Como resultado do ensaio, as áreas modificadas estão descritas na sequência.
Treliça 1: 7250 mm².
Treliça 2: 2850 mm².
Treliça 3: 20000 mm².
Percebe-se então que é possível reduzir ainda mais a área das barras verticais de apoio, porém, este fato não influenciaria o estado de tensão da estrutura. No entanto, a diminuição das áreas acarretou um aumento considerável no deslocamento, excedendo o valor permitido em norma e comprometendo a resistência estrutural da estrutura
Além disso, realiza-se a análise de flambagem para os novos valores de área de seção transversal, culminando nos valores de tensão crítica em 1138,8 MPa e 447,7 MPa para os elementos superiores e inferiores, respectivamente. Logo, a estrutura otimizada não irá flambar.
Os resultados da simulação em relação ao deslocamento e tensão estão demonstrados nas figuras 10 e 11.
Fig. 10. Magnitude da tensão
Fig. 11. Magnitude dos deslocamentos em mm
Referências
[1] ALVES FILHO, Alvelino. ELEMENTOS FINITOS A BASE DA TECNOLOGIA CAE (2013)
 [2] AZEVEDO, Álvaro. MÉTODO DOS ELEMENTOS FINITOS (2003). Disponível em: <http://www.arquivoescolar.org/bitstream/arquivo-e/117/1/Livro_MEF.pdf>
[3] LEET, Kenneth; UANG, Chia-Ming Uang; GILBERT, Anna M. FUNDAMENTO DA ANÁLISE ESTRUTURAL 3º edição (2010)
[4] TECHNICAL DESCRIPTION RAIL MOUNTED GANTRY. Disponível em: <https://www.liebherr.com/shared/media/maritime-cranes/downloads-and-brochures/brochures/lcc/liebherr-rmg-technical-description.pdf>
[5] ABNT NBR 8800. PROJETO DE ESTRUTURAS DE AÇO E DE ESTRUTURAS MISTAS DE AÇO E CONCRETO DE EDIFÍCIOS (2008) Disponível em: <http://www.inti.gob.ar/cirsoc/pdf/acero/NBR8800_2008_1.pdf>
[6] RLT Engenharia. TRELIÇAS: DEFINIÇÃO, APLICAÇÃO E CLASSIFICAÇÃO (2015) 
Disponível em: <http://rltengenharia.blogspot.com.br/2015/11/trelicas-classificacao.html>
[7] FORNARI, José. DESENVOLVIMENTO DE UM SOFTWARE DE APOIO AO ESTUDO DE TRELIÇAS EM PONTES E TELHADOS (2002)
Disponível em: <https://repositorio.ufsc.br/bitstream/handle/123456789/82947/189264.pdf?sequence=1>
Anexos
Anexo A – Deslocamentos máximos verticais para estruturas de aço conforme a norma NBR (8800/2008) [5] 
Apêndices
Apêndice A – Comparação entre os modelos de treliça Pratt versus Warren
Pratt
Warren
Apêndice B – Programa Scilab
clc // Clear screen
clear all// Clear all variables in memory
// -----------------------------------------------------------------
// Dados de entrada do problema: otimizacao trabalho
// Definicao dos nos (coordenadas)
No=[9000 1000;
9000 0;
8000 0;
8000 1000;
10000 0;
7000 1000;
6000 1000;
5000 1000;
6000 0;
5000 0;
4000 0;
7000 0;
4000 1000;
3000 1000;
2000 1000;
1000 1000;
2000 0;
1000 0;
0 0;
3000 0;
0 1000;
-1000 1000;
-2000 1000;
-3000 1000;
-2000 0;
-3000 0;
-4000 0;
-1000 0;
-4000 1000;
-5000 1000;
-6000 1000;
-7000 1000;
-7000 0;
-8000 0;
-5000 0;
-6000 0;
-8000 1000;
-9000 1000;
-10000 0;
-9000 0] 
// Definicao dos elementos (conectividade)
El=[1 2;
3 2;
3 4;
1 4;
3 1;
2 5;
6 3;
6 7;
7 8;
8 9;
10 9;
8 10;
11 8;
11 10;
4 6;
5 1;
12 3;
6 12;
9 6;
9 12;
9 7;
8 13;
11 13;
14 11;
14 15;
15 16;
16 17;
18 17;
16 18;
19 16;
19 18;
13 14;
20 11;
14 20;
17 14;
17 20;
17 15;
16 21;
19 21;
22 19;
22 23;
23 24;
24 25;
26 25;
24 26;
27 24;
27 26;
21 22;
28 19;
22 28;
25 22;
25 28;
25 23;
24 29;
27 29;
30 27;
30 31;
31 32;
32 33;
34 33;
34 32;
29 30;
35 27;
30 35;
36 30;
32 36;
36 31;
36 35;
32 37;
33 36;
38 34;
38 39;
39 40;
40 34;
34 37;
37 38;
38 40]
 // Elemento 1: nos 1 e 2
// Definicao das propriedades da barra
Prop=[200000 15000; // E[MPa] A[mm^2]
200000 20000;
200000 15000;
200000 15000;
200000 15000;
200000 20000;
200000 15000;
200000 15000;
200000 15000;
200000 15000;
200000 20000;
200000 15000;
200000 15000;
200000 20000;
200000 15000;
200000 15000;
200000 20000;
200000 15000;
200000 15000;
200000 20000;
200000 15000;
200000 15000;
200000 15000;
200000 15000;
200000 15000;
200000 15000;
200000 15000;
200000 20000;
200000 15000;
200000 15000;
200000 20000;
200000 15000;
200000 20000;
200000 15000;
200000 15000;
200000 20000;
200000 15000;
200000 15000;
200000 15000;
200000 15000;
200000 15000;
200000 15000;
200000 15000;
200000 20000;
200000 15000;
200000 15000;
200000 20000;
200000 15000;
200000 20000;
200000 15000;
200000 15000;
200000 20000;
200000 15000;
200000 15000;
200000 15000;
200000 15000;
200000 15000;
200000 15000;
200000 15000;
200000 20000;
200000 15000;
200000 15000;
200000 20000;
200000 15000;
200000 15000;
200000 15000;
200000 15000;
200000 20000;
200000 15000;
200000 20000;
200000 15000;
200000 15000;
200000 20000;
200000 20000;
200000 15000;
200000 15000;
200000 15000] // E[MPa] A[mm^2]
// Condicao de contorno dp=[gl1, gl2, gl4]
dp=[3,4,79,80] // Restricao dos graus de liberdade aonde está pinado 2 de cada xy
// Forcas aplicadas
// -----------------------------------------------------------------
// Contar quantidade de nos do elemento
Num_nos_el=size(El,"c") // Contarnumero de colunas da matriz El
// Contar quantidade de elementos
Num_el=size(El,"r") // Contar numero de linhas da matriz El
// Contar quantidade de nos
Num_nos=size(No,"r") // contar numero de linhas da matriz No
// Quantidade de graus de liberdade por no
Num_gl_no=2 // Deslocamento na direcao x e y para elemento de barra
// Quantidade total de graus de liberdade
Num_gl=Num_nos*Num_gl_no
// Matriz de rigidez global (mesmo tamanho dos graus de liberdade total)
KK=zeros(Num_gl,Num_gl)
// Matriz para definir graus de liberdade para cada elemento
Edof=zeros(Num_el,4);
for i=1:Num_el
 Edof(i,1)=El(i,1)*2-1;
 Edof(i,2)=El(i,1)*2;
 Edof(i,3)=El(i,2)*2-1;
 Edof(i,4)=El(i,2)*2;
end
// Montagem da matriz de forca
F=zeros(Num_gl,1); // Criando matriz forca nos nos
F(38,1) = -400000 //N GL xy
// ------------------------------------------------------------------
// Calculos para cada elemento
for i=1:Num_el
 
 // Identificar a numeracao do no nesse elemento
 no_1=El(i,1)
 no_2=El(i,2)
 
 // Identificar a posicao no espaco
 x1=No(no_1,1); y1=No(no_1,2)
 x2=No(no_2,1); y2=No(no_2,2)
 
 // Calcular angulo em relacao ao eixo X global
 if(x2-x1)==0 // Elemento na vertical
 if(y2>y1) 
 theta(i)=90; // [degrees]
 else // y2<y1
 theta(i)=-90; // [degrees]
 end
 elseif(y2-y1)==0 // Elemento na horizontal
 if(x2>x1)
 theta(i)=0; // [degrees]
 else // x2<x1
 theta(i)=180; // [degrees]
 end
 elseif x2<x1 
 if y2<y1 then
 theta(i)=180+atand((y2-y1)/(x2-x1))
 else //y2>y1
 theta(i)=180+atand((y2-y1)/(x2-x1))
 end
 elseif x2>x1
 if y2<y1 then
 theta(i)=360+atand((y2-y1)/(x2-x1))
 else //y2>y1
 theta(i)=atand((y2-y1)/(x2-x1))
 end
 end
 // Comprimento de cada elemento como matriz L=[L1; L2; L3]
 L(i,1)=sqrt((x2-x1)^2+(y2-y1)^2)
 
 // Matriz de transformacao
 c=cosd(theta(i))
 s=sind(theta(i))
 // Matriz de transformacao no formato 4x4
 // T=[c -s 0 0;s c 0 0; 0 0 c -s; 0 0 s c] 
 T=[c s 0 0;0 0 c s]
// Mostrar matriz de transformacao
disp("Elemento:"); disp(i)
disp("Matriz de transformacao:"); disp(T)
// ------------------------------------------------------------------
 // Matriz de rigidez do elemento de barra (quadratura de gauss)
 
 // Funcao de interpolacao linear no dominio [-1,1]
 npi=4 // Numero de pontos de integracao
 k_l=0 // Valor inicial, em formato de matriz quad (ngliberdade)
 
 if npi==1 then
 Gauss=[0 2]
 elseif npi==2 then
 Gauss=[-sqrt(1/3) 1;
 +sqrt(1/3) 1]
 elseif npi==3 then 
 Gauss=[-sqrt(3/5) (5/9);
 0 (8/9);
 +sqrt(3/5) (5/9)]
 elseif npi==4 then 
 Gauss=[-sqrt((3/7)-(2/7)*sqrt(6/5)) (18+sqrt(30))/36;
 +sqrt((3/7)-(2/7)*sqrt(6/5)) (18+sqrt(30))/36;
 -sqrt((3/7)+(2/7)*sqrt(6/5)) (18-sqrt(30))/36;
 +sqrt((3/7)+(2/7)*sqrt(6/5)) (18-sqrt(30))/36] 
 elseif npi==5 then 
 Gauss=[0 (128/225);
 -(1/3)*sqrt(5-2*sqrt(10/7)) (322+13*sqrt(70))/900; 
 +(1/3)*sqrt(5-2*sqrt(10/7)) (322+13*sqrt(70))/900;
 -(1/3)*sqrt(5+2*sqrt(10/7)) (322-13*sqrt(70))/900;
 +(1/3)*sqrt(5+2*sqrt(10/7)) (322-13*sqrt(70))/900] 
 end
 
 for p=1:npi // De 1 ate o numero de pontos de integracao
 
 ksi(p)=Gauss(p,1) // Localizacao do ponto (Coluna 1)
 w(p)=Gauss(p,2) // Peso (Coluna 2)
 
 // Apresentacao das Funcoes de forma (interpolacao)
 N=[(1/2)*(1-ksi(p)) (1/2)*(1+ksi(p))]
 
 // Derivada das funcoes de forma (interpolacao)
 dNdksi=[(1/2)*(-1) (1/2)*(+1)]
 
 // Criar um vetor com coordenadas do elemento, sistema local
 x=[0; L(i)] // x=[0; comprimento do elemento]
 
 // Jacobiano
 J=dNdksi*x //J=L(i)/2 
 
 // Matriz deformacao-deslocamento
 B=[dNdksi(1)*(1/J) dNdksi(2)*(1/J)] // B=[-1/L(i) 0 1/L(i) 0]
 
 // Modulo de elasticidade
 E=Prop(i,1)
 // Area
 A=Prop(i,2)
 // Matriz constitutiva C
 C=E*A
 
 // Integrando
 f=B'*C*B*J // avaliado nos pontos de integracao. nesse caso e constante!
 // Matriz de rigidez
 k_l=k_l+f*w(p)
 end
 
 // Matriz de rigidez do elemento de barra (ja na forma direta)
 // k_l=(Prop(i,1)*Prop(i,2)/L(i,1))*[1 0 -1 0; 0 0 0 0; -1 0 1 0; 0 0 0 0] // Local
// Mostrar matriz de rigidez do elemento no sistema de coord local.
disp("Matriz de rigidez do elemento (local):"); disp(k_l)
// ------------------------------------------------------------------
 k_g=T'*k_l*T // Global. Este "'" diz respeito a Transposta
// Mostrar matriz de rigidez do elemento no sistema de coord global.
disp("Matriz de rigidez do elemento (global):"); disp(k_g)
 
 // Montagem da matriz de rigidez
 // Contribuicao do elemento i na matriz de rigidez global
 KK(Edof(i,:),Edof(i,:)) = k_g + KK(Edof(i,:),Edof(i,:))
end
// Mostrar matriz
disp("Matriz de rigidez global:"); disp(KK)
K1 = KK // armazenando matriz K para pos-processamento
// ------------------------------------------------------------------
// Aplicando as condicoes de contorno na matriz Rigidez
// Método da penalização
 Pen=max(KK)+10e8 // Valor do termo com penalização
 Num_BC=size(dp,"c") // Contagem do tamanho do vetor
 
 for i=1:Num_BC 
 KK(dp(i),dp(i))=Pen // Substituir termos pela penalização
 end
 
// ------------------------------------------------------------------
// Resultado da analise: campo de deslocamentos nodais
U=(KK^-1)*F // Solucao do sistema de equacoes
// Mostrar resultados
disp("Deslocamentos nodais [mm]:"); disp(U)
// ------------------------------------------------------------------
// Pos-processamento
// Forca de reacao nos nos com gl restrito
RF=K1*U
disp("Forcas nos nos [N]:"); disp(RF)
// Calculo de tensao no elemento, conforme Avelino, pag. 93. Equacao 3.12
Force=zeros(Num_el,1)
for i=1:Num_el
 // Considerando que o vetor theta ja foi definido
 c=cosd(theta(i))
 s=sind(theta(i))
 t=[c s] // Matriz de transformacao
 deltaU(1,i)= U(Edof(i,3))-U(Edof(i,1)) // U2-U1 direcao x
 deltaU(2,i)= U(Edof(i,4))-U(Edof(i,2)) // V2-V1 direcao y 
 // Calculo da forca interna no elemento
 Force(i)=Prop(i,1)*Prop(i,2)/L(i,1)*t*deltaU(:,i)
 // Calculo de tensao normal
 Sigma(i)=Force(i)/Prop(i,2) // Forca dividido pela area da secao
 // Considerando que é uma barra quadrada com área A, o momento de inércia é:
 I(i)=(Prop(i,2)^2)/12
 
 // Carga critica de flambagem
 P_cr(i)=-(%pi^2)*Prop(i,1)*I(i)/(L(i)^2)
 if(Force(i))<0 // Elementos com esforços compressivos
 if(Force(i))<P_cr(i) // comparativo com carga crítica
 disp("Flambagem no elemento:"); disp(i)
 end
 end
end
 disp("Tensao em cada elemento [MPa]:"); disp(Sigma) 
 disp("Forca interna em cada elemento [N]:"); disp(Force)
 disp("Carga critica de flambagem [N]:"); disp(P_cr)