Nonlinear analysis of plane frames using a corotational formulation and plasticity by layers

Prévia do material em texto

NONLINEAR ANALYSIS OF PLANE FRAMES USING A CO-
ROTATIONAL FORMULATION AND PLASTICITY BY LAYERS 
 
Sebastião S. da Silva, William T. M. Silva 
PECC, Postgraduate Program in Structure and Civil Construction Engineering, Department 
of Civil and Environmental Engineering, University of Brasília -UnB, Campus Darcy Ribeiro, 
70910900, Brasília-DF, Brazil, http://www.estruturas.unb.br/ 
Keywords: Nonlinear analysis of structures; Co-rotational formulation; Plasticity in layers; 
Finite element method. 
Abstract. The purpose of this work is to perform a nonlinear analysis of plane frame 
structures using a co-rotational formulation and a layered plastic modeling. The plane frame is 
discretized with 2D beam elements. The main objective is to compare three different local 
elastoplastic elements. Plasticity is introduced by rate-independent one-dimensional model 
with isotropic hardening. Numerical integration over the cross-section is performed for obtain 
the internal force vector and tangent stiffness matrix of these elements. At each integration 
point, the return-mapping algorithm is used for integration of the constitutive equations. A 
some examples are used in order to assess the performances of the elements and the path 
following procedures. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1 INTRODUÇÃO 
A análise linear apresenta dificuldades em refletir o real comportamento de estruturas 
menos comuns, sob condições de carregamento não ordinárias ou ainda próximas ao colapso. 
Nesse contexto, a análise não linear física e geométrica tem grande aplicabilidade na 
engenharia estrutural, já que possibilita o conhecimento da capacidade resistente 
remanescente, depois da análise ultrapassar os pontos críticos e/ou atingir deformações 
inelásticas finitas. 
Para Menin (2006), na análise não-linear geométrica utilizando o método dos elementos 
finitos, três diferentes tipos de descrições cinemáticas têm sido amplamente utilizadas: 
descrição lagrangiana total, descrição lagrangiana atualizada e descrição co-rotacional. Esta 
última de acordo com Reddy (2004) tem origem no teorema da decomposição polar o qual 
estabelece que a deformação total de uma superfície contínua pode ser decomposta em 
movimento de corpo rígido e deformação relativa. O conceito da descrição cinemática co-
rotacional foi introduzido em um contexto do MEF com os trabalhos de Argyris (1965) e 
segundo Menin (2006), sendo a mais recente das formulações utilizadas na análise não linear, 
conta com uma grande variedade de assuntos a ser pesquisada. 
A não-linearidade física, por sua vez, está diretamente relacionada ao comportamento 
mecânico dos materiais constituintes da estrutura. Quando a resposta do sólido que constitui a 
estrutura é elástica, após cessada a carga o corpo não apresenta deformações. Já quando a 
resposta é plástica, após cessado o carregamento a deformação não se desfaz, sendo portanto 
irreversível ou permanente - Owen e Hinton (1980) e, Lemaitre e Chaboche (1994). A teoria 
da plasticidade fornece leis e modelos capazes de descrever o comportamento constitutivo de 
materiais que apresentam resposta elastoplástica. 
2 METODOLOGIA EMPREGADA NO TRABALHO 
Utiliza-se dois tipos de elementos de viga 2D na discretização dos pórticos - o modelo 
matemático de viga Euler-Bernoulli (C1) sem o acoplamento dos esforços axiais e de flexão, a 
partir de agora denominado elemento de viga Bernoulli e; o modelo anterior com acoplamento 
entre os esforços axiais e de flexão aqui denominado elemento de viga Bernoulli 2. 
Para realizar as simulações numéricas com a formulação proposta implementa-se no 
programa de elementos finitos 2D_Beam_f90 os elementos de viga descritos anteriormente 
utilizando a formulação co-rotacional. Este programa foi desenvolvido no âmbito do 
Programa de Pós Graduação em Estruturas e Construção Civil (PECC/UnB). 
No tratamento da plasticidade por camadas foi usado um modelo de plasticidade 
unidimensional bilinear com encruamento isotrópico. Integra-se este modelo utilizando um 
algoritmo implícito denominado na literatura de “return mapping”, Simo e Hughes (1998). 
Para obter os esforços seccionais, integram-se as tensões normais, utilizando sete ou quinze 
pontos de Gauss ao longo da altura da seção transversal. 
Os vetores de forças internas e a matriz de rigidez tangente dos três elementos acima 
mencionados foram determinados através do princípio dos trabalhos virtuais. Os coeficientes 
destes vetores e matrizes foram adquiridos por integração numérica. Para obter a trajetória de 
equilíbrio não-linear adota-se uma análise incremental-iterativa baseada no método de 
Newton-Raphson e na técnica do comprimento de arco (arc-length). Neste trabalho 
implementou-se seções transversais retangulares e trabalhou-se com outros tipos de seções 
transformando-as em uma seção retangular equivalente (de mesma área e inércia). 
3 A FORMULAÇÃO CO-ROTACIONAL DE PÓRTICOS PLANOS 
Sendo a energia de deformação elástica um polinômio de segunda ordem, para se obter a 
solução exata da integração, usou-se dois pontos de Gauss ao longo do eixo dos elementos. 
Suas posições e pesos são dados respectivamente por: 1 2(1 1 3)x L= − ; 2L e 
2 2(1 1 3)x L= + ; 2L . Os números de pontos de Gauss empregados ao longo da altura da 
seção transversal, foram sete e quinze. A formulação descrita aqui é a apresentada em Battini 
(2002), a qual se diferencia ligeiramente da exibida em Crisfield (1991). 
3.1 Descrição cinemática de um elemento de viga qualquer 
Na Fig. 1 apresenta-se as coordenadas dos nós 1 e 2 no sistema de coordenadas 
global( ),x z , ( )1 1,x z , ( )2 2,x z . O vetor de deslocamento global é definido por: 
 gu { }1 1 1 2 2 2θ θ=
T
u w u w . (1) 
1
2
x
z
θ
θ
θ
θ1
1
2
2
u1
w1
u
w2
2
zl
xl
0
u
 
Figura 1: Formulação co-rotacional – cinemática da viga. 
 
Por outro lado, o vetor de deslocamentos locais é dado por: 
 lu { }1 2θ θ= Tu . (2) 
As componentes de lu podem ser calculadas através das seguintes equações, 
 0= −nu l l , 1 1θ θ α= − e 2 2θ θ α= − . (3) 
Nas equações acima 0l e nl denotam os comprimentos inicial e atual, respectivamente. 
Após se fazer simples relações trigonométricas podem-se reescrevê-las como, 
 ( ) ( )
1
2 2 2
0 2 1 2 1
 = − − −
 
l x x z z e ( ) ( )
1
2 2 2
2 2 1 1 2 2 1 1
 = + − − − + − −
 n
l x u x u z w z w (4) 
onde α denota a rotação de corpo rígido. Esta pode estar relacionada trigonometricamente por 
 0 0α = −sen c s s c e 0 0cosα = +c c s s, (5) 
as quais, após desenvolvidas, chega-se a: 
 ( )0 0 2 1
0
1
cosβ= = −c x x
l
, ( )0 0 2 1
0
1β= = −s sen z z
l
, 
 ( )2 2 1 1
1
cosβ= = + − −
n
c x u x u
l
 e ( )2 2 1 1
1β= = + − −
n
s sen z w z w
l
. (6) 
Assim, para um α π< , α é dado por: 
 ( )1α α−= sen sen se 0α ≥sen e cos 0α ≥ ; ( )1cos cosα α−= se 0α ≥sen e cos 0α < ; 
 ( )1α α−= sen sen se 0α <sen e cos 0α ≥ e ( )1cos cosα α−= − se 0α <sen e cos 0α < (7) 
3.1.1 Deslocamentos virtuais 
Diferenciando-se as Eq. (3) obtêm-se os deslocamentos virtuais locais, 
 ( ) ( ) { }2 1 2 1 0 0δ δ δ δ δ= − + − = − −u c u u s w w c s c s guδ , (8) 
 1 1 1δθ δθ δα δθ δβ= − = − ( )0α β β= − e (9) 
 2 2 2δθ δθ δα δθ δβ= − = − .(10) 
Por outro lado, δβ pode ser calculada pela diferenciação da Eq. (6d) 
 ( ) ( )2 1 2 2 1 12
1δβ δ δ δ=  − − + − −  n n
n
w w l z w z w l
cl
, (11) 
onde nlδ = uδ é dado pela Eq. (8). Usando a Eq. (6d), a expressão de δβ torna-se 
 ( ) ( ) ( )22 1 2 1 2 1
1δβ δ δ δ δ δ δ = − − − − − n
n
w w l sc u u s w w
cl
, (12) 
a qual, após simplificações, produz 
 { }1 0 0δβ = − −
n
c s c s
l g
uδ . (13) 
Desta forma, a matriz de transformação B é definida como, 
 luδ = B uδ g , (14) 
e dada por, 
 B =
0 0
1 0
0 1
− − 
 − − 
 − − 
n n n n
n n n n
c s c s
s l c l s l c l
s l c l s l c l
. (15) 
3.1.2 Vetor de forças internas 
A relação entre o vetor de forças internas local lf e global fg é obtida pela equação dos 
trabalhos virtuais nos sistemas local e global 
 u f u f u B fδ δ δ= = =T T T Tg g l l g lV (16) 
A Eq. (16) deve ser aplicada a um uδ g arbitrário, e assim, o vetor de forças internas global 
gf é dado por 
 f B f= Tg l (17) 
no qual o vetor de forças internas local { }f = Tl 1 2N M M depende da definição do elemento 
finito de viga específico empregado. 
3.1.3 Matriz de rigidez tangente 
A matriz de rigidez tangente global gK definida por 
 g g gf K uδ δ= (18) 
é obtida pela variação da Eq. (17). Assim, 
 Tg l 1 1 2 2 3N b M b M bf B fδ δ δ δ δ= + + + . (19) 
Na equação anterior, 2b é, por exemplo, a segunda coluna de 
TB . As seguintes notações 
são introduzidas, 
 { }0 0 Tc s c sr = − − e { }0 0 Ts c s cz = − − (21) 
as quais por diferenciação tornam-se, 
 r zδ δβ= e z rδ δβ= − . (22) 
As Eq. (8) e (13) podem ser reescritas como, 
 Tn gu l r uδ δ δ= = e 
1 T
g
nl
z uδβ δ= . (23) 
Introduz-se expressões auxiliares, 
 1b r= , { }2
1
0 0 1 0 0 0
T
nl
b z= − e { }3
1
0 0 0 0 0 1
nl
T
b z= − (24) 
as quais diferenciadas produzem, 
 1
T
g
nl
zz
b r pδ δ δ= = e ( )2 3 1 T Tn g2 2
n n n
l
l l l
zz
b b rz zr p
δδδ δ δ= = − + = + . (25) 
 
O primeiro termo na Eq. (19) é calculado pela introdução da matriz de rigidez tangente 
local lK , a qual depende da definição do elemento. 
 l l l l gf K u K B uδ δ δ= = (26) 
Finalmente, das Eq. (18), (19), (25a), (25b) e (26), a expressão da matriz de rigidez 
tangente global torna-se, 
 ( )( )21
n
T T T T
g l 1 2
n
N
M M
l l
K B K B zz rz zr= + + + + (27) 
As Eq. (14) e (27) fornecem a relação entre forças internas e matrizes de rigidez tangente 
locais e globais. Estas relações independem da definição local do elemento. O último é obtido 
adotando-se como hipóteses: deformação de Bernoulli; relação constitutiva linear elástica e; o 
princípio dos trabalhos virtuais (PTV) - 1 1 2 2vV dv N u M Mσδε δ δθ δθ= = + +∫ . 
3.2 Elemento de viga Bernoulli 
O campo de deslocamento deste elemento é dado por: 
 
x
u u
L
= e 
2 22
1 21 1
x x x
w x
L L L
θ θ   = − + −   
   
. (28) 
Assim, a curvatura 2 2k w x= ∂ ∂ e a deformação u x kzε = ∂ ∂ − são dadas por, 
 1 22 2
4 2
6 6θ θ   = − + + − +   
   
x x
k
L L L L
 e 1 22 2
4 2
6 6ε θ θ    = + − + −    
    
u x x
z
L L L L L
. (29) 
 
3.2.1 Obtenção numérica do vetor de forças internas 
 
Utilizando o PTV, chega-se as forças internas locais, as quais são dadas por: 
 
v
N dv
L
σ= ∫ , 1 2
4
6
v
x
M z dv
L L
σ  = − 
 
∫ e 2 2
2
6
v
x
M z dv
L L
σ  = − 
 
∫ (30) 
Integra-se as Eq. (30) ao longo do eixo-x e sobre a altura da seção transversal do elemento, 
obtendo-se os coeficientes do vetor de forças internas no formato implementado – Quadro 1. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Quadro 1: Algoritmo para obtenção do vetor de forças internas local – Elemento de viga Bernoulli. 
 
As tensões são calculadas no âmbito elástico (ˆ =E E ) e, caso estas violem determinadas 
condições de admissibilidade, um rotina que considera a elastoplasticidade (ˆ = tE E ) é 
chamada. ( )pg i é a coordenada do i-ésimo ponto de Gauss e ( )wg i seu respectivo peso. 
 
0N =
 
1 0M = 
2 0M = 
 
faça 1,i = número de pontos de Gauss
 
 
( )
2
h
z pg i= ×
 
 
( ) ( )1 1 21 3 3 1ε θ θ = + + + − 
u z
l l 
 
( ) ( )1 1 21 3 1 3ε θ θ = + − − + 
u z
l l 
 
caso elástico: 1 1Eσ ε= 
 caso elastoplástico: chamar algoritmo “return mapping” → 1σ 
 caso elástico: 2 2Eσ ε= 
 caso elastoplástico: chamar algoritmo “return mapping” → 2σ 
 
( ) ( )
( ) ( )
( ) ( )
1 2
1 1 1 2
2 2 1 2
1
2 2 2
1 3 1 3
2 2 2 2
3 1 1 3
2 2 2 2
h h
N N b wg i b wg i
h h
M M z b wg i z b wg i
h h
M M z b wg i z b wg i
σ σ
σ σ
σ σ
 = + × × × + × × × 
 
+ −   = + × × × × + × × × ×   
   
− +   = + × × × × − × × × ×   
   
 
fim faça
 
3.2.2 Obtenção numérica da matriz de rigidez tangente 
 
A matriz de rigidez tangente local analítica lK é obtida pela diferenciação das Eq. (30). 
Aplicando integração gaussiana na direção do eixo e ao longo da altura da seção transversal 
do elemento, se obtêm os coeficientes tal qual programado - (Quadro 2). 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Quadro 2: Algoritmo para obtenção numérica da matriz de rigidez tangente – Elemento de viga Bernoulli. 
 
3.3 Elemento de viga Bernoulli 2 
O elemento de viga Bernoulli 2 usa ( )21 [ 1 2 ]f Lkz L u x w x dx kzε ε= − = ∂ ∂ + ∂ ∂ −∫ como 
definição para deformação local. Entretanto, toma-se uma deformação axial média para evitar 
o travamento de membrana (membrane locking). Assim: 
 2 21 1 2 2 1 22 2
1 1 1 4 2
6 6
15 30 15
u x x
z
L L L L L
ε θ θ θ θ θ θ    = + − + + − + −    
    
. (31) 
A fim de simplificar as expressões subsequentes serão usadas as seguintes variáveis 
auxiliares: 2 21 1 2 2
1 1 1
15 30 15
u
g
L
θ θ θ θ= + − + , 1 1 2
2 1
15 30
g θ θ= − , 2 2 1
2 1
15 30
g θ θ= − . 
 
3.3.1 Obtenção numérica do vetor de forças internas 
 
Fazendo uso do PTV se chega as equações componentes do vetor de forças internas: 
11 0k = 
22 0k = 
23 0k = 
33 0k = 
 
faça 1,i = número de pontos de Gauss
 
 
( )
2
h
z pg i= ×
 
 caso elástico: 1ˆ =E E 
 caso elastoplástico: 1ˆ = +
EH
E
E H 
 caso elástico: 2ˆ =E E 
 caso elastoplástico: 2ˆ = +
EH
E
E H 
 
( ) ( )11 11 1 2
1 ˆ ˆ
2 2 2
h h
k k E b wg i E b wg i
l
 = + × × × + × × × 
  
 
( ) ( ) () ( )
2 2
2 2
22 22 1 2
1 3 1 3
ˆ ˆ
2 2 2 2
h h
k k E z b wg i E z b wg i
l l
+ −   = + × × × × + × × × ×   
    
 
( ) ( )2 223 23 1 2
1 1ˆ ˆ
2 2
h h
k k E z b wg i E z b wg i
l l
   = + × × × × + × × × ×   
    
 
( ) ( ) ( ) ( )
2 2
2 2
22 22 1 2
3 1 1 3
ˆ ˆ
2 2 2 2
h h
k k E z b wg i E z b wg i
l l
− +   = + × × × × + × × × ×   
    
fim faça 
 
v
N dv
L
σ= ∫ , 1 1 2
4
6
v v
x
M g dv z dv
L L
σ σ  = + − 
 
∫ ∫ e 2 2 2
2
6
v v
x
M g dv z dv
L L
σ σ  = + − 
 
∫ ∫ . (32) 
Integrando numericamente as Eq. (32) ao longo do eixo-x e da altura seção transversal do 
elemento encontra-se as expressões como implementadas - (Quadro 3). 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Quadro 3: Algoritmo para obtenção numérica do vetor de forças internas local – Elemento de viga Bernoulli 2. 
 
 
3.3.2 Obtenção numérica da matriz de rigidez tangente 
 
A matriz de rigidez tangente local é calculada pela diferenciação das Eq. (32). Integrando 
esta ao longo do eixo-x e sobre a altura da seção se obtêm os coeficientes de rigidezes da 
matriz como programado - Quadro 4. 
 
 
 
 
 
 
 
 
0N =
 
1 0M = 
2 0M = 
 
faça 1,i = número de pontos de Gauss
 
 
( )
2
h
z pg i= ×
 
 
( ) ( )1 1 21 3 3 1ε θ θ = + + + − 
z
g
l 
 
( ) ( )2 1 21 3 1 3ε θ θ = + − − + 
z
g
l 
 caso elástico: 1 1ˆσ ε= E 
 caso elastoplástico: chamar algoritmo “return mapping” → 1σ 
 caso elástico: 2 2ˆσ ε= E 
 caso elastoplástico: chamar algoritmo “return mapping” → 2σ 
 
( ) ( )1 2
1
2 2 2
h h
N N b wg i b wg iσ σ = + × × × + × × × 
  
 
( ) ( )1 1 1 1 22 2 2
l h h
M M g b wg i b wg iσ σ = + × × × × + × × × + 
  
 
( ) ( )1 2
1 3 1 3
2 2 2 2
h h
z b wg i z b wg iσ σ+ −   + × × × × + × × × ×   
    
 
( ) ( )2 2 2 1 22 2 2
l h h
M M g b wg i b wg iσ σ = + × × × × + × × × + 
  
 
( ) ( )1 2
3 1 1 3
2 2 2 2
h h
z b wg i z b wg iσ σ− +   + × × × × − × × × ×   
    
fim faça 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Quadro 4: Algoritmo para obtenção numérica da matriz de rigidez tangente – Elemento de viga Bernoulli 2. 
 
1σ 
2σ 
11 0k = 
12 0k = 
13 0k = 
22 0k = 
23 0k = 
33 0k = 
 
faça 1,i = número de pontos de Gauss
 
 
( )
2
h
z pg i= ×
 
 caso elástico: 1ˆ =E E 
 caso elastoplástico: 1ˆ = +
EH
E
E H 
 caso elástico: 2ˆ =E E 
 caso elastoplástico: 2ˆ = +
EH
E
E H 
( ) ( )11 11 1 2
1 ˆ ˆ
2 2 2
h h
k k E b wg i E b wg i
l
 = + × × × + × × × 
  
( ) ( )12 12 1 1 2
1 ˆ ˆ
2 2 2
h h
k k g E b wg i E b wg i = + × × × + × × × 
  
( ) ( )13 13 2 1 2
1 ˆ ˆ
2 2 2
h h
k k g E b wg i E b wg i = + × × × + × × × 
  
 
( ) ( )222 22 1 1 2
1 ˆ ˆlg
2 2 2
h h
k k E b wg i E b wg i = + × × × + × × × + 
  
 
( ) ( ) ( ) ( )
2 2
2 2
1 2
1 3 1 3
ˆ ˆ
2 2 2 2
h h
E z b wg i E z b wg i
l l
+ −   + × × × × + × × × × +   
    
 
( ) ( )1 215 2 2
l h h
b wg i b wg iσ σ + × × × + × × × 
  
 
( ) ( )23 23 1 2 1 2
1 ˆ ˆlg g
2 2 2
h h
k k E b wg i E b wg i = + × × × + × × × + 
  
 
( ) ( )2 21 2
1 1ˆ ˆ
2 2
h h
E z b wg i E z b wg i
l l
   + × × × × + × × × × −   
    
 
( ) ( )1 260 2 2
l h h
b wg i b wg iσ σ − × × × + × × × 
  
 
( ) ( )33 33 1 2ˆ ˆ2 2
h h
k k E b wg i E b wg i = + × × × + × × × + 
  
 
( ) ( ) ( ) ( )
2 2
2 2
1 2
3 1 1 3
ˆ ˆ
2 2 2 2
h h
E z b wg i E z b wg i
l l
− +   + × × × × + × × × × +   
    
 
( ) ( )1 215 2 2
l h h
b wg i b wg iσ σ + × × × + × × × 
  
fim faça 
4 MODELO DE ELASTOPLASTICIDADE UNIDIMENSIONAL 
Na Figura 2 ilustra-se o encruamento, fenômeno observado experimentalmente em muitos 
metais. Este causa uma expansão do espaço de tensões admissíveis σS . 
 
Figura 2: Plasticidade com encruamento. Fonte: Simo e Hughes (1998), adaptado. 
 
Na expansão de σS assume-se que: o encruamento é isotrópico; o encruamento é linear no 
ramo do fluxo plástico e independe de ( )εɺ psign . A primeira condição leva a: 
 ( ), 0σ µ σ σ µ = − + ≤ y yf H , 0µ ≥ , (33) 
onde: H é o módulo plástico e; µ é a variável de encruamento. A equação de evolução mais 
simples para esta variável é µ ε= ɺɺ p . O mecanismo que governa sua evolução é definido por: 
 ( )ε ϕ σ=ɺ p sign . (34) 
Este é capturado por meio das condições de Kuhn-Tucker, 
 0ϕ ≥ , ( ), 0σ µ ≤yf e ( ), 0ϕ σ µ =yf , (35) 
onde 0ϕ ≥ é determinada pela condição de consistência 
 ( ), 0ϕ σ µ =ɺyf . (36) 
Esta permite que se resolva explicitamente ϕ e se relacione as taxas de tensão com as de 
deformação. Das Eq (33) e (34), juntamente com a relação tensão-deformação elástica, 
 ( ) [ ] 0σ ε ϕ= − + ≤ɺ ɺyf sign E E H (37) 
Observe que 0>ɺyf não pode manter-se. De (34) e (36) segue que 0ϕ ≠ apenas se 
 0= =ɺy yf f ⇒ 
( )σ εϕ =
+
ɺsign E
E H
. (38) 
Então a forma do intervalo da relação elástica juntamente com (38) produz, 
 σ =
0
0
ε ϕ
ε ϕ
=

 > +
ɺ
ɺ
E se
EH
se
E H
, (39) 
onde a quantidade +EH E H é chamada módulo tangente elastoplástico. 
Algoritmo de mapeamento e retorno 
O processo de integração do modelo elastoplástico ocorre localmente. Emprega-se o 
método backward-Euler que leva ao algoritmo de mapeamento e retorno. Este considera 
inicialmente o seguinte estado auxiliar elástico teste obtido pelo congelamento do fluxo 
plástico: 
 
( )1 1σ ε ε σ ε+ += − ≡ + ∆teste pn n n n nE E , 
 1ε ε+ =
p teste p
n n , (40) 
 1µ µ+ =
teste
n n e 
 ( 1) 1σ σ µ+ +  = − + 
teste teste
y n n y nf H . 
Observa-se que o estado teste é determinado apenas em termos das condições iniciais 
{ }, ,ε ε µpn n n e de ε∆ n . Calculado este estado, duas situações podem ocorrer – Tabela 1. 
 
Tabela 1: Situações possíveis para um estado de tensão teste. 
 Função teste Admissibilidade Condições de Kuhn-Tucker Processo 
S
itu
aç
ão
 (
i) 
. 
( 1) 0+ ≤
teste
y nf Estado teste admissível, pois: 
1ε ε+ =
p p
n n , 
1µ µ+ =n n , 
1 1σ σ+ +=
teste
n n .
 
 
 ( 1) ( 1)
0+ += ≤
teste
y n y nf f 
e 
0ϕ∆ ≡ 
 
Elástico 
 
Fig. 3a 
S
itu
aç
ão
 (
ii)
 
. 
( 1) 0+ >
teste
y nf Estado teste inadmissível, pois viola: 
( ), 0σ µ ≤yf 
Assim: 
0ϕ∆ > , 
1ε ε+ ≠
p p
n n , 
1 1σ σ+ +≠
teste
n n . 
 
 
 0ϕ∆ > 
e 
( 1) 0ϕ +∆ =y nf ⇒ ( 1) 0+ =y nf , 
 
Plástico 
 
 
Fig. 3b 
 
 
 
 (a) (b) 
Figura 3: Passos incrementais: (a) elástico. (b) plástico. Fonte: Simo e Hughes (1998). 
 
Seja o problema incrementalmenteplástico [situação (ii)] caracterizado pelas condições, 
 ( )( 1) 1 10 , 0σ µ+ + +> ⇔ =testey n n nf f e 0ϕ∆ > . (41) 
O objetivo é determinar a solução { }1 1 1, , ,ε µ σ ϕ+ + + ∆pn n n . Para tanto, se expressa a tensão final 
1σ +n em termos de 1σ +
teste
n e ϕ∆ como segue: 
 ( )1 1 1σ σ ϕ σ+ + += − ∆testen n nE sign . (42) 
Dessa forma, desde que 0ϕ∆ > e, tendo em vista a Eq. (42): 
 
( )
( )
1 1 1
1 1
1
( 1) 1 1
,
,
0.
σ σ ϕ σ
ε ε ϕ σ
µ µ ϕ
σ σ µ
+ + +
+ +
+
+ + +
 = − ∆
 = + ∆
 = + ∆
  = − + = 
teste
n n n
p p
n n n
n n
y n n y n
E sign
sign
e
f H
 (43) 
Agora o problema da Eq. (43) é resolvido explicitamente em termos do estado elástico 
teste pelo seguinte procedimento: 
 ( ) ( ) ( )1 1 1 1 1testen n n n nsign sign Esignσ σ σ σ ϕ σ+ + + + += − ∆ , (44) 
 ( ) ( )1 1 1 1σ ϕ σ σ σ+ + + + + ∆ =  teste testen n n nE sign sign . 
Desde que 0ϕ∆ > e 0>E , os termos dentro dos colchetes na Eq. (44)b são 
necessariamente positivos, ( ) ( )1 1σ σ+ += testen nsign sign e 1 1σ ϕ σ+ ++ ∆ = testen nE . Finalmente 0ϕ∆ > 
é determinado da condição da Eq. (43)4 considerando a Eq. (40). O critério de escoamento é 
( 1)+y nf é dado por: 
 ( 1) 1 1
σ ϕ σ µ+ + + = − ∆ − + 
teste
y n n y nf E H ou ( )( 1) ( 1) ϕ+ += − ∆ +testey n y nf f E H . (45) 
Assim, 
 ( 1) 0+ =y nf ⇒ 
( 1) 0ϕ +∆ = >
+
teste
y nf
E H
. (46) 
Substituindo a Eq. (46) na Eq. (43), 
 
1 1
1
1
ϕσ σ
σ+ ++
 ∆
 = −
  
teste
n nteste
n
E
. (47) 
O algoritmo de mapeamento e retorno é resumido no Quadro 5. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Quadro 5: Algoritmo de mapeamento e retorno. 
 
 i) Dados em { }: ,ε µ∈ pn nx B 
 ii) Dado um campo de deformação em ∈x B : 1ε ε ε+ = + ∆n n n 
 iii) Calcule uma tensão elástica teste e verifique para o carregamento plástico 
 ( )1 1σ ε ε+ += −teste pn n nE 
 ( 1) 1σ σ µ+ +  = − + 
teste teste
y n n y nf H 
 SE ( 1) 0+ ≤
teste
y nf , então :FAZER 
 Passo elástico: Fixar ( ) ( )1 1+ +• = •
teste
n n
→ !SAIR 
 SE NÃO: 
 Passo plástico: )IR PARA PASSO iv 
 FIM SE. 
 iv) Mapeamento e retorno 
 ( 1) 0ϕ +∆ = >
+
teste
y nf
E H
 1 1
1
1
ϕσ σ
σ+ ++
 ∆
 = −
  
teste
n nteste
n
E
 
 ( )1 1ε ε ϕ σ+ += + ∆p p testen n nsign 1µ µ ϕ+ = + ∆n n 
5 EXEMPLOS NUMÉRICOS 
Selecionou-se quatros exemplos numéricos com os objetivos de estudar eficiência da 
formulação dos elementos implementados em capturar grandes deslocamentos e rotações na 
presença da elastoplasticidade; realizar uma análise comparativa entre as duas formulações de 
elementos de viga propostas; verificar a influência do número de pontos de Gauss na seção 
transversal; investigar o efeito do número de elementos adotados na discretização das 
estruturas estudadas e; analisar a eficiência da adoção de seções retangulares equivalentes no 
estudo de outros tipos de seções transversais. As estruturas escolhidas – pórtico de Lee, viga 
em balanço, pórtico assimétrico e pórtico toggle – são clássicas da literatura e, dessa forma, 
tem servido como benchmark para diversos trabalhos. 
5.1 Pórtico de Lee 
As condições de carregamento, de contorno, de geometria e os parâmetros materiais da 
estrutura analisada estão mostrados na Fig. 4. Esta foi apresentada inicialmente por Lee 
(1968) para análise de grandes deslocamentos e estabilidade de pórticos elásticos. Battini 
(2002), por sua vez, utilizou esse pórtico para estudar elementos de viga co-rotacional em 
problemas de instabilidade. 
 As configurações das trajetórias de equilíbrio da estrutura são conhecidas na literatura 
como snap-through no grau de liberdade u e v apresentando, portanto dois pontos limites no 
diagrama carga-deslocamento. 
 
P
3,0
2,0
120
24
12
0
v
u
 
Figura 4: Geometria, condições de contorno e parâmetros físicos utilizados no pórtico de Lee. 
 
Obteve-se a resposta elastoplástica apresentada pela estrutura fazendo-se uso de 20 
elementos de viga Bernoulli e procedendo-se a comparação com Battini (2002). Conforme 
pode-se visualizar na Fig. 5, há uma grande aproximação entre os resultados. Verifica-se 
ainda uma queda acentuada da capacidade portante da estrutura. Todavia, nota-se uma 
elevação da curva elastoplástica no segundo ponto limite como onsequência do 
encruamento do material. 
Realizou-se ainda uma análise da estrutura usando 7 e 15 pontos de Gauss ao longo da 
altura da viga. Conforme exibido nas curvas da Fig. 6, há uma boa aproximação de resultados 
entre as duas verificações, o que evidencia uma integração numérica eficiente das equações 
E = 720 
υ = 0,3 
Et = E/10 
σy = 10,44 
com a adoção de apenas 7 pontos. Deformadas da estrutura são exibidas na Fig. 7, onde se 
pode verificar que as acentuadas rotações e translações ocorridas são capturadas pela 
formulação implementada. 
 
 
 
Figura 5: Resposta elastoplástica do pórtico de Lee para o elemento de viga Bernoulli e uma malha com 20 
elementos. 
 
 
 
Figura 6: Resposta elastoplástica do pórtico de Lee para o elemento de viga Bernoulli 2; malha com 20 
elementos e; 7 e 15 pontos de Gauss. 
 
 
Figura 7: Perfil de deformadas para a resposta elastoplástica. 
 
 
5.2 Viga em balanço 
Este problema foi resolvido inicialmente por Euler em 1744. O fato de se conhecer a sua 
solução analítica e, pelas características de seu comportamento, leva o mesmo a ser muito 
utilizado para testes de formulações destinadas a análise não linear. Remo (2000) utiliza essa 
estrutura nos seus estudos de pórticos inelásticos que experimenta grandes deslocamentos. As 
dimensões reais e equivalentes, as condições de contorno e os parâmetros físicos do material 
usados neste exemplo estão mostrados na Fig. 8. Avaliou-se o efeito provocado por uma 
carga pontual aplicada na extremidade livre da estrutura. 
 
 
400
P
v
u
R r = 17,393
r R = 17,773
33
33
R - r2 2
R + r4 4
R + r2 2
 
Figura 8: Dimensões (real e equivalente), condições de contorno e propriedades do material usadas na viga em 
balanço. 
 
A trajetória de equilíbrio da resposta elastoplástica é mostrada na Fig. 9. Os resultados 
concordam com os obtidos por Remo (2000) o que mostra a eficácia do uso da estratégia de 
seções retangulares equivalentes. Observa-se ainda que, para a hipótese de material linear, os 
resultados obtidos são exatamente os mesmos dos obtidos pelo citado autor. 
E = 20 
 υ = 0,3 
Et = E/30 
 
Figura 9: Resposta elastoplástica para viga em balanço usando elemento de viga Bernoulli 2; malha com 10 
elementos e 15 pontos de Gauss. 
 
5.3 Pórtico Toggle 
As condições de carregamento, geométrica (real e equivalente), de contorno, e de material 
do pórtico toggle aqui estudado pode ser visto na Fig. 10. A necessidade de se verificar a 
eficiência de formulações numéricas no tratamento de trajetórias de equilíbrio com caráter 
"snap-through" justifica o fato de essa estrutura constar como exemplo em grande número de 
referências bibliográficas. Ela foi resolvida inicialmente de maneira analítica e experimental 
por Williams (1964). Remo (2000) também a estudouna investigação de grandes 
deslocamentos experimentados por pórticos inelásticos. 
 
Figura 10: Dimensões (real e equivalente), condições de contorno e propriedades do material usadas no pórtico 
toggle. 
As respostas elástica e elastoplástica do pórtico são mostradas na Fig. 11 onde se pode 
perceber o caráter “snap-through” das curvas e a existência de dois pontos limites. Observa-
se uma queda na capacidade portante da estrutura quando se contabiliza na sua trajetória de 
equilíbrio os efeitos da elastoplasticidade. Plota-se sobre as trajetórias de equilíbrio elásticas e 
elastoplásticas pontos extraídos de Remo (2000) e conclui-se que os resultados obtidos entre 
os dois trabalhos são praticamente coincidentes. Assim, pode-se depreender que o uso de uma 
seção retangular equivalente é eficiente. Para a curva com tensão de escoamento infinito, as 
curvas elásticas se sobrepõem, exatamente. 
E = 199714 
G = 76813 
 υ = 0,3 
Et = E/2 
σy = 3x103 
65,715 
 
Figura 11: Resposta elástica e elastoplástica para o pórtico toggle usando elemento de viga Bernoulli 2; malha 
com 8 elementos e 15 pontos de Gauss. 
 
5.4 Pórtico assimétrico 
O último exemplo estudado trata-se de um pórtico assimétrico composto de duas barras de 
mesmo comprimento, uma vertical e a outra horizontal, ligadas rigidamente entre si. As 
condições de contorno, de geometria e os parâmetros materiais são mostrados na Fig. 12. 
Battini (2002) utiliza essa estrutura com a finalidade de estudar elementos de viga formulados 
com a descrição co-rotacional em problemas de instabilidade. 
 
P
0,1
0,1
1,0
1,
0
d
 
Figura 12: Dimensões, condições de contorno e propriedades do material usadas do pórtico assimétrico. 
 
Estuda-se a resposta da estrutura submetida a uma carga pontual aplicada sobre ligação 
rígida na direção vertical. Como o programa 2D_beam_nl.f90 não possui estratégias 
computacionais para obtenção das trajetórias secundárias de equilíbrio, induziu-se a resposta 
pós flambagem nas duas direções possíveis, aplicando uma carga pontual de pequena 
 
E = 2x1011 
 υ = 0,3 
Et = E/2 
σy = 10
6 
magnitude na direção do grau de liberdade d. Na Fig. 13 exibem-se as trajetórias 
elastoplástica onde se observa a existência de um ponto de bifurcação do qual se iniciam as 
trajetórias assimétricas. 
 
 
Figura 13: Resposta elastoplástica para o pórtico assimétrico usando elemento de viga Bernoulli 2; malha com 
16 elementos e 15 pontos de Gauss. 
 
 
6 CONCLUSÕES 
 
Neste trabalho avaliou-se a eficiência de dois elementos de viga 2D aplicados em um 
programa de elementos finitos que leva em consideração a não-linearidade geométrica (e 
física. Nos exemplos estudados – pórtico de Lee, pórtico assimétrico, viga em balanço e 
pórtico toggle - verifica-se uma queda na capacidade portante da estrutura quando se 
considera o fenômeno da elastoplasticidade. Além disso, evidencia-se a capacidade da 
formulação co-rotacional em capturar grandes deslocamentos e rotações para todos os 
elementos implantados. 
Na análise elastoplástica, o uso de quinze pontos de Gauss na altura da seção transversal 
em vez de sete não produz uma melhora significativa dos resultados. Os resultados 
apresentados mostram que a estratégia utilizada para evitar enrijecimento por membrana no 
elemento Bernoulli 2 foi eficiente. No exemplo do pórtico assimétrico deve-se mensurar com 
cuidado o valor da carga indutora da resposta no sentido desejado - minorá-la pode fazer com 
que o programa não consiga obter a trajetória secundária; majorá-la pode simular a resposta 
da estrutura com imperfeição. Por fim, depreende-se que o uso da simples estratégia de seções 
retangulares equivalentes no estudo de outros tipos de seções transversais fornecem resultados 
satisfatórios. 
 
Agradecimentos 
 
Ao CNPq pelo apoio financeiro, essencial para viabilizar a pesquisa no nosso país. Ao 
Programa de Pós-Graduação em Estruturas e Construção Civil - UnB. 
 
 
REFERÊNCIAS 
 
Argyris, J.H., (1965), “Continua and discontinua”, Proceedings 1st Conference on Matrix 
Methods in Structural Mechanics, AFFDL-TR-66-80, Air Force Institute of Tecnology, 
Dayton, Ohio-USA. 
Battini, J.M., (2002), “Co-rotational beam elements in instability problems”, Ph.D Thesis, 
Royal Institute of Tecnology - Departament of Mechanics, Stockholm / Sweeden. 
Crisfield, M.A., (1991), “Non-linear finite element analysis of solids and structures”, Volume 
1: Essential, John Wiley & Sons, Chichester, UK. 
De Souza, R. M., (2000), "Force-Based Finite Element for Large Displacement Inelastic 
Analsysis of Frames", Ph.D. Dissertation, University of California at Berkeley, Berkeley, 
CA, USA. 
Felippa, C.A., (2001), “Non-linear finite element methods / NFEM”, Lecture notes for the 
course non-linear finite element methods, Center for Aerospace Structures, University of 
Colorado, Boulder/USA. 
J. N. Reddy., (2004), “An Introduction to Nonlinear Finite Element Analysis”, Oxford 
University Press : Oxford, U.K. 
Owen, D.R.J., Hinton, E., (1980), Finite Elements in Plasticity: Theory and Practice, 
Pineridge. Swansea, U.K. 
Lemaitre, J., Chaboche, J-L., (1994), Mechanics of Solids Materials, Cambrige University 
Press, Cambrige, , U.K. 
Menin, R. C. G., (2006), “Aplicação da descrição co-rotacional na análise não-linear 
geométrica de estruturas discretizadas por elementos finitos de treliças, vigas e cascas”, 
Tese de doutorado, E.TD-004A/06, Brasília : ENC/FT/UnB. 
Simo, J. C., (1998), Computational Inelastic, Interdisciplinary Apllied Mathematics, vol. 7, 
Springer, New York, USA. 
Williams, F. W. (1964), ‘An approach to the nonlinear behaviour of the members of a rigid 
jointed plane framework with finite deflection’, Quart. J. Mech. Appl. Maths. 17(4), 451–
469.