Capítulo 5 - Formulação Corrotacional do Método dos

Professor Luiz Antonio Farani de Souza

1 Capítulo 5 - Formulação Corrotacional do Método dos Elementos Finitos para treliça 2D

Capítulo 5 - Formulação Corrotacional do Método dos Elementos

Finitos para treliça 2D

Conteúdo Capítulo 5 - Formulação Corrotacional do Método dos Elementos Finitos para treliça 2D ......... 1

5.1 Introdução ........................................................................................................................... 2

5.2 Formulação Corrotacional de Elementos Finitos ................................................................ 2

5.3 Problema numérico resolvido ............................................................................................. 3

5.3.1 Solução do problema ................................................................................................... 4

5.4 Programa computacional .................................................................................................... 5

5.4.1 Programa principal ....................................................................................................... 6

5.4.2 Sub-rotinas (funções) ................................................................................................... 7

5.5 Exercícios propostos .......................................................................................................... 10

5.5.1 Exercício proposto 1 ................................................................................................... 10

5.5.2 Exercício proposto 2 ................................................................................................... 11

Referências .............................................................................................................................. 13



5.1 Introdução

Neste capítulo, a formulação e o código computacional são apresentados para análise não

linear geométrica de treliças planas, com o programa livre Scilab, versão 6.1.1 (SCILAB, 2021).

O programa principal e as sub-rotinas (funções) são apresentados.

Na análise de treliça plana, consideram-se a não linearidade geométrica e a medida de

deformação de Green-Lagrange. A estrutura é discretizada por meio do Método Corrotacional de

Elementos Finitos (YAW, 2009).

O sistema de equações não lineares é resolvido pelo método incremental e iterativo de

Newton-Raphson Padrão associado à técnica de continuação Comprimento de Arco Linear

(hiperplano fixo).

O programa é composto por três partes:

a) Entrada de dados (pré-processamento);

b) Análise do problema (processamento); e

c) Saída de dados (pós-processamento).

5.2 Formulação Corrotacional de Elementos Finitos

A estrutura da treliça pode sofrer grandes deslocamentos e rotações no nível global, mas

pequenas deformações no nível local. Considere o elemento de treliça nas configurações inicial e

corrente como mostrado na Figura 1.

Figura 1: Configurações inicial e corrente do elemento de treliça.

Fonte: Adaptada de Yaw (2009).

Os comprimentos inicial L0 e corrente L são determinados conforme as equações,

respectivamente:

𝐿0 = √(𝑋2 − 𝑋1)2 + (𝑌2 − 𝑌1)

2, (5.1)

𝐿 = √(𝑥2 − 𝑥1)2 + (𝑦2 − 𝑦1)

2. (5.2)

As coordenadas globais permanecem fixas na formulação corrotacional. O ângulo

corrente do elemento corrotacionado com relação ao sistema de coordenadas globais é denotado

por . Na formulação bidimensional, os valores do seno e cosseno desse ângulo são determinados

por, respectivamente (YAW, 2009; CRISFIELD, 1991):

𝑐𝑜𝑠(𝜃) =(𝑋2 + 𝑢2) − (𝑋1 + 𝑢1)

𝐿, (5.3)

𝑠𝑒𝑛(𝜃) =(𝑌2 + 𝑣2) − (𝑌1 + 𝑣1)

𝐿. (5.4)

nas quais ui são os deslocamentos horizontais (na direção do eixo cartesiano global X) e vi são os

deslocamentos verticais (na direção do eixo cartesiano global Y), com i = 1, 2.



Considere 𝑐 = 𝑐𝑜𝑠(𝜃) e 𝑠 = 𝑠𝑒𝑛(𝜃). O vetor de força interna elementar Fel é avaliado

pela seguinte expressão:

𝑭𝒆𝒍 =𝑁𝐿

𝐿0𝒓, (5.5)

na qual o vetor 𝒓 = [−𝑐 −𝑠 𝑐 𝑠]𝑇 e N é a força axial interna local dada por:

𝑁 = 𝐴𝐸𝜀𝐺 , (5.6)

em que EA é a rigidez axial e 𝜀𝐺 é a deformação de Green-Lagrange dada por:

𝜀𝐺 =𝐿2 − 𝐿0

2

2𝐿02 , (5.7)

A força axial N é calculada no sistema local de coordenadas cartesianas (Xl, Yl). A matriz

de rigidez elementar Kel é composta por duas parcelas:

𝑲𝒆𝒍 = 𝑲𝑴 +𝑲𝑮, (5.8)

na qual KM é a matriz de rigidez material e KG é a matriz de rigidez geométrica. Essas matrizes

sendo determinadas por, respectivamente:

𝑲𝑴 =𝐴𝐸

2𝐿0(3

𝐿2

𝐿02 − 1)𝒓𝒓𝑇 , (5.9)

𝑲𝑮 =𝑁

𝐿

𝐿

𝐿0[𝑰𝟐 −𝑰𝟐−𝑰𝟐 𝑰𝟐

] =𝐴𝐸𝜀𝐺𝐿0

𝑪, (5.10)

Na qual a matriz C é dada por:

𝑪 = [𝑰𝟐 −𝑰𝟐−𝑰𝟐 𝑰𝟐

] . (5.11)

em que I2 é a matriz identidade de ordem 2.

5.3 Problema numérico resolvido

Considere a estrutura correspondente a uma barra biarticulada com um grau de liberdade,

conforme é ilustrada na Figura 2. A barra tem rigidez axial adimensional EA = 5,0 107. Este

problema foi estudado por Menin (2006) e Crisfield (1991).



Figura 2: Modelo estrutural da barra com um grau de liberdade.

5.3.1 Solução do problema

Dados de entrada para o método de solução (Newton-Raphson padrão associado à técnica

continuação comprimento de Arco Linear):

Barra de progresso (processamento):



Trajetória de equilíbrio (deslocamento vertical versus força):

Resultados numéricos no Console do Scilab (ktotal, kmédio e t):

5.4 Programa computacional

Para a criação do programa para a solução do problema de barra descrito na Seção 5.3,

deve-se seguir os passos:

1. Executar o programa Scilab 6.1.1;

2. Criar uma pasta com o nome do programa (por exemplo, "Treliça_2D_Corr");

3. Abrir o programa Editor de texto SciNotes;

4. Criar no SciNotes o arquivo do programa principal (Subseção 5.4.1) com o nome “principal” e

extensão .sce, dentro da pasta criada; e

5. Criar novos arquivos no SciNotes para cada uma das funções (function) (Subseção 5.4.2) com

extensão .sci, dentro da pasta criada;



5.4.1 Programa principal //Programa principal - treliça 2D //Formulação corrotacional

//Análise não linear geométrica

//Deformação de Green-Lagrange clear

clc

exec('DKG.sci',-1); exec('DFG.sci',-1);

exec('dkelem.sci',0);

exec('ensamkg.sci',-1); exec('contkg.sci',-1);

exec('dfelem.sci',-1);

exec('ensamfg.sci',-1); exec('contfg.sci',-1);

exec('apontador.sci',-1);

exec('result.sci',-1); exec('entrada_dados.sci',-1);

exec('result.sci',-1);

//_____________________________________

//entrada de dados

//tol - tolerância

//deltal - comprimento de arco inicial //Nd - número desejável de iterações

//kmax - número máximo de iterações por passo

//nmax - número máximo de passos de carga //P - incremento de carga

txt = ['tolerância:';'número máximo de iterações:';'número de passos de carga:';'comprimento de arco inicial:';'número de

iterações desejadas por passo de carga:';'incremento de carga:']; sig = x_mdialog('Parâmetros método de solução',txt,['10^-8';'150';'25';'110';'5';'-1'])

tol = evstr(sig(1));

kmax = evstr(sig(2)); nmax = evstr(sig(3));

deltal = evstr(sig(4));

Nd=evstr(sig(5)); P=evstr(sig(6));

//Entrada de dados (construção da malha de EF) [NTNOS, NTEL, NNOSCC, NTGL, coord, inci, dofno, Fr, NOCC, E0, A, itipo]=entrada_dados(P);

//_____________________________________ //Processamento

//inicialização

lambda=0; deltal0=deltal;

coord0=coord;

udesl=zeros(NTGL,1); //vetor de deslocamento deltau=zeros(NTGL,1);

DELTAU=zeros(NTGL,1);

vu(1,1)=0; vf(1,1)=0;

ktotal=0; ierro=0;

tic //inicia um cronômetro winH=waitbar('Processamento ...'); //inicia barra de progresso

realtimeinit(0);

for np=1:nmax //passos de carga

[KG]=DKG(udesl,NOCC,NNOSCC,NTGL,NTEL,dofno,inci,coord0,E0,A,itipo); //matriz de rigidez

deltaur=KG\Fr;

Dlambda=deltal/norm(deltaur); if DELTAU'*deltaur<0 //determina o sinal do incremento de carga

Dlambda=-Dlambda;

end DELTAU0=Dlambda*deltaur;

DELTAU=DELTAU0;

[FG,vdef]=DFG(udesl+DELTAU,NOCC,NNOSCC,NTGL,NTEL,dofno,inci,coord0,E0,A,itipo); //vetor de força interna g=(lambda+Dlambda)*Fr-FG; //vetor de forças desequilibradas

k=0;

realtime(np); while k<kmax //ciclo iterativo

k=k+1; //contador de iterações

[KG]=DKG(udesl+DELTAU,NOCC,NNOSCC,NTGL,NTEL,dofno,inci,coord0,E0,A,itipo); //matriz de rigidez deltaug=KG\g;

deltaur=KG\Fr;



dlambda=-(DELTAU0'*deltaug)/(DELTAU0'*deltaur); //subincremento de carga

deltau=deltaug+dlambda*deltaur; //subincremento de deslocamento

DELTAU=DELTAU+deltau; //incremento de deslocamento Dlambda=Dlambda+dlambda; //incremento de carga

[FG,vdef]=DFG(udesl+DELTAU,NOCC,NNOSCC,NTGL,NTEL,dofno,inci,coord0,E0,A,itipo); //vetor de força interna

g=(lambda+Dlambda)*Fr-FG; //vetor de forças desequilibradas if norm(g)<=norm(Fr)*tol //critério de convergência

break

end end

if k==kmax

messagebox('não convergiu!') ierro=1;

break end

udesl=udesl+DELTAU; //deslocamento total

lambda=lambda+Dlambda; //parâmetro de carga total deltal = deltal0*(Nd/k)^0.5;

vu(1+np,1)=-udesl(4,1);

vf(1+np,1)=-lambda*Fr(4,1);

ktotal=ktotal+k; //contador iterações acumuladas

waitbar(lambda,winH); //barra de progresso

end close(winH); //fecha a barra de progresso

kmedio=ktotal/nmax; //número médio de iterações por passo

t=toc() //lê o cronômetro //_____________________________________

//Saída de Dados (pós-processamento)

if ierro==0 //gráfico - trajetória de equilíbrio

[vu1,vf1]=result();

plot(vu,vf,'b-','marker','s','markerFaceColor','b','markerEdgeColor','b','markersize',6);

set(gca(),"auto_clear","off");

plot(vu1,vf1,'k','marker','x','markersize',8) set(gca(),"auto_clear","off");

gca().grid=[1 1 1]; //Linhas de grade xlabel('Deslocamento vertical','fontsize',3); //eixo x

ylabel('Força P','fontsize',3); //eixo y

legend('Programa EF','Menin (2006)',3);

//resultados numéricos (console)

disp('Resultados numéricos') disp('a) Número total de iterações (ktotal):',ktotal)

disp('b) Número médio de iterações por passo (kmédio):',kmedio)

disp('c) Tempo de processamento em segundos (t):',t) end

5.4.2 Sub-rotinas (funções)

Função apontador function [IPO, TAM]=apontador(m, itipo) //Elemeno de barra com 2 NÓS, 2GL/NÓ

if (itipo(m,2)==1) then IPO(1)=1;

IPO(2)=2;

IPO(3)=3; IPO(4)=4;

TAM=4;

end endfunction

Função contfg function [FG]=contfg(NOCC, NNOSCC, FG)

//Impõe as condições de contorno no vetor de força global FG for I=1:NNOSCC

FG(NOCC(1,I))=0;

end endfunction



Função contkg function [KG]=contkg(NOCC, NNOSCC, NTGL, KG) //Impõe as condições de contorno na matriz de rigidez global KG

for J=1:NNOSCC

for I=1:NTGL KG(NOCC(1,J),I)=0;

KG(I,NOCC(1,J))=0;

end KG(NOCC(1,J),NOCC(1,J))=1;

end

endfunction

Função entrada_dados function [NTNOS, NTEL, NNOSCC, NTGL, coord, inci, dofno, Fr, NOCC, E0, A, itipo]=entrada_dados(P)

//Informa os dados da malha de elementos finitos, vetor de força e propriedades materiais e geométricas das barras

//NTNOS -> NÚMERO TOTAL DE NÓS //NTEL -> NÚMERO TOTAL DE ELEMENTOS

//NTGL -> NÚMERO TOTAL DE GRAUS DE LIBERDADE

//NNOSCC -> NÚMERO DE GRAUS RESTRITOS (CONDIÇÕES DE CONTORNO) NTNOS=2; //informar

NTEL=1; //informar

NNOSCC=3; //informar NTGL=NTNOS*2;

//Informa as coordenadas dos nós (informar a matriz) //coord(i,1)= coordenada x

//coord(i,2)= coordenada y

coord=[0 0;2500 2500]; //informar

//Informa a incidência dos elementos (informar a matriz)

//inci(i,1) = elemento //inci(i,2) = nó i

//inci(i,3) = nó j

inci=[1 1 2]; //informar

//determina os graus de liberdade por nó de cada elemento

for I=1:NTEL //elemento de treliça 2D

dofno(I,1)=inci(I,2)*2-1; //NÓ i

dofno(I,2)=inci(I,2)*2; dofno(I,3)=inci(I,3)*2-1; //NÓ j

dofno(I,4)=inci(I,3)*2;

end

//elemento barra -> itipo(nel,2)==1

//elemento viga -> itipo(nel,2)==2 //material 1 -> itipo(nel,3)==1

//material 2 -> itipo(nel,3)==2

for i=1:NTEL itipo(i,1) = i;

itipo(i,2) = 1;

itipo(i,3) = 1; end

//propriedades materiais e geométricas das barras

for m=1:NTEL

if itipo(m,3)==1

E0(m)=5*10^7; //informar o módulo de elasticidade

A(m)=1; //informar a área da seção transversal da barra

end end

//vetor de força externa Fext Fr=zeros(NTGL,1);

//informar as forças (carregamento) nos graus de liberdade

Fr(2*2,1)=P; //nó * 2gl (força na direção do eixo Y)

//informar as condições de contorno (graus de liberdade restringidos)

NOCC=[1 2 3]; //informar os graus de liberdade restringidos endfunction



Função dfelem function [FELEM, eps]=dfelem(m, E0, A, inci, udesl, dofno, coord0) //Determina o vetor de força elementar

FELEM=zeros(4,1);

for i=1:4 u(i)=udesl(dofno(m,i),1);

end

X1=coord0(inci(m,2),1); X2=coord0(inci(m,3),1);

Y1=coord0(inci(m,2),2);

Y2=coord0(inci(m,3),2); L0(m)=sqrt((X2-X1)^2+(Y2-Y1)^2); //comprimento indeformado

L(m)=sqrt( (X2+u(3)-X1-u(1))^2 + (Y2+u(4)-Y1-u(2))^2 ); //comprimento deformado

c=(X2+u(3)-X1-u(1))/L(m); s=(Y2+u(4)-Y1-u(2))/L(m);

r=[-c;-s;c;s];

eps=(L(m)^2-L0(m)^2)/(2*L0(m)^2); //deformação Green Lagrange FELEM = E0(m)*A(m)*eps*(L(m)/L0(m))*r; //vetor de força elementar

endfunction

Função DFG function [FG, vdef]=DFG(udesl, NOCC, NNOSCC, NTGL, NTEL, dofno, inci, coord0, E0, A, itipo) //Determinação do vetor de força interna global (FG)

FG=zeros(NTGL,1);

for m=1:NTEL [FELEM,eps]=dfelem(m,E0,A,inci,udesl,dofno,coord0);

[FG]=ensamfg(m,FELEM,dofno,itipo,FG);

vdef(m,1)=eps; end

[FG]=contfg(NOCC,NNOSCC,FG);

endfunction

Função dkelem function [KELEM]=dkelem(m, E0, A, inci, udesl, dofno, coord0)

//Monta a matriz de rigidez elementar Kelem

KELEM=zeros(4,4);

for i=1:4

u(i)=udesl(dofno(m,i),1);

end X1=coord0(inci(m,2),1);

X2=coord0(inci(m,3),1);

Y1=coord0(inci(m,2),2); Y2=coord0(inci(m,3),2);

L0(m)=sqrt((X2-X1)^2+(Y2-Y1)^2); //comprimento indeformado

L(m)=sqrt( (X2+u(3)-X1-u(1))^2 + (Y2+u(4)-Y1-u(2))^2 ); //comprimento deformado c=(X2+u(3)-X1-u(1))/L(m);

s=(Y2+u(4)-Y1-u(2))/L(m);

r=[-c -s c s]; C = [1 0 -1 0;

0 1 0 -1;

-1 0 1 0; 0 -1 0 1];

const=L(m)/L0(m);

eps=(L(m)^2-L0(m)^2)/(2*L0(m)^2); //deformação Green Lagrange KM = E0(m)*A(m)/L0(m)*0.5*(3*const^2-1)*r*r';

Kg = E0(m)*A(m)*eps/L(m)*L(m)/L0(m)*C;

KELEM=KM+Kg; //matriz de rigidez

endfunction

Função DKG function [KG]=DKG(udesl, NOCC, NNOSCC, NTGL, NTEL, dofno, inci, coord0, E0, A, itipo) //Determinação da matriz de rigidez global (KG)

KG=zeros(NTGL,NTGL);

for m=1:NTEL [KELEM]=dkelem(m,E0,A,inci,udesl,dofno,coord0);

[KG]=ensamkg(m,KELEM,dofno,itipo,KG);

end [KG]=contkg(NOCC,NNOSCC,NTGL,KG);

endfunction



Função ensamfg function [FG]=ensamfg(m, FELEM, dofno, itipo, FG) //Monta o vetor de força global

[IPO,TAM]=apontador(m,itipo);

for I=1:TAM P=dofno(m,IPO(I));

if (P>0)

FG(P,1)=FG(P,1)+FELEM(I,1); end

end

endfunction

Função ensamkg function [KG]=ensamkg(m, KELEM, dofno, itipo, KG)

//Monta a matriz de rigidez global

[IPO,TAM]=apontador(m,itipo); for I=1:TAM

for J=1:TAM

P=dofno(m,IPO(I)); Q=dofno(m,IPO(J));

if (P>0 & Q>0)

KG(P,Q)=KG(P,Q)+KELEM(I,J); end

end

end endfunction

Função result function [vu1, vf1]=result()

//resultados deformação de Green-Lagrange (Menin, 2006): com o programa Pega Ponto vu1=[7.80E+02

1.26E+03

1.76E+03 3.27E+03

3.78E+03

4.27E+03];

vf1=[3.16E+06

3.28E+06

2.41E+06 -2.49E+06

-3.32E+06

-3.07E+06]; endfunction

5.5 Exercícios propostos

5.5.1 Exercício proposto 1

Adaptar o programa para as seguintes medidas de deformação:

- Deformação de Engenharia;

- Deformação de Green-Lagrange;

- Deformação de Almansi; e

- Deformação de Biot.

Resolver o problema da barra biarticulada com um grau de liberdade ilustrado na Figura

1, e comparar as trajetórias de equilíbrio para as diferentes medidas de deformação. A solução do

problema (trajetórias de equilíbrio) apresentada por Menin (2006) aparece na Figura 3.

Tabela 1: Matriz de rigidez e o vetor de força interna elementares para as diferentes medidas de

deformação.

Deformação Fel Kel = KM + KG

Engenharia

(MENIN, 2006) 𝑭𝐞𝐥 = 𝐸𝐴𝑒𝒓

𝑲𝑴 =𝐸𝐴

𝐿0𝒓𝒓𝑻

𝑲𝑮 = 𝐸𝐴𝑒𝑪



Green-Lagrange

(MENIN, 2006) 𝑭𝐞𝐥 = 𝐸𝐴𝐺

𝐿

𝐿0𝒓

𝑲𝑴 =𝐸𝐴

2𝐿0(3

𝐿2

𝐿02 − 1)𝒓𝒓𝑻

𝑲𝑮 = 𝐸𝐴𝐺𝐿

𝐿0𝑪

Biot

(MENIN, 2006) 𝑭𝐞𝐥 = 𝐸𝐴𝐵

𝐿0𝐿𝒓

𝑲𝑴 =𝐸𝐴

𝐿(3

𝐿02

𝐿2− 2

𝐿0𝐿)𝒓𝒓𝑻

𝑲𝑮 = 𝐸𝐴𝐵𝐿0𝐿𝑪

Almansi

(MENIN, 2006) 𝑭𝐞𝐥 = 𝐸𝐴𝐴𝐿0

2

𝐿2𝒓

𝑲𝑴 =𝐸𝐴

2𝐿(5

𝐿04

𝐿4− 3

𝐿02

𝐿2)𝒓𝒓𝑻

𝑲𝑮 = 𝐸𝐴𝐴𝐿0

2

𝐿2𝑪

Figura 3: Trajetórias de equilíbrio.

Fonte: Menin (2006).

5.5.2 Exercício proposto 2

Seja a treliça plana com dez barras, cuja geometria e o carregamento são apresentados na

Figura 4. As barras têm módulo de elasticidade E = 5,0 × 1010 N/m2 e área da seção transversal A

= 1,0 × 10-4 m². Obter a trajetória de equilíbrio para a medida de deformação de Engenharia. A

solução (trajetória de equilíbrio com dois pontos limites de força) é apresentada na Figura 5.

Figura 4: Modelo estrutural da treliça abatida com 10 barras.

Fonte: Adaptada de Wriggers, Wagner e Miehe (1988).



Figura 5: Trajetória de equilíbrio.



Referências

CRISFIELD, M. A. Non-linear finite element analysis of solids and structures. Vol. 1. New

York: Wiley, 1991.

MENIN, R. C. G. Aplicação da descrição cinemática co-rotacional na análise não-linear

geométrica de estruturas discretizadas por elementos finitos de treliças, vigas e cascas. Tese

(Doutorado em estruturas e construção civil) - Departamento de Engenharia Civil e Ambiental,

Faculdade de Tecnologia, Universidade de Brasília, Brasília, 2006.

ROCHA, G. Estratégias de incremento de carga e de iteração para análise não linear de

estruturas. Dissertação (Mestrado) - Departamento de Engenharia Civil, Universidade Federal

de Ouro Preto, Ouro Preto, 2000.

SCILAB, versão 6.1.1. France: ESI Group, 2021.

WRIGGERS, P.; WAGNER, W.; MIEHE, C. A quadratically convergent procedure for the

calculation of stability points in finite element analysis. Computer methods in applied

mechanics and engineering, v. 70, n. 3, p. 329-347, 1988.

YAW, L. L. 2D Co-rotational Truss Formulation. Walla Walla University, 2009.

Documents

Capítulo 5 - Formulação Corrotacional do Método dos