Introdução ao Método dos Elementos Finitos: Elasticidade ...luis/ae2/AE2_Elasticidade2D_v4.pdf · Análise de Estruturas II: Elasticidade Plana e Tridimensional 3 2. Hipóteses

Análise de Estruturas II: Elasticidade Plana e Tridimensional

1

Introdução ao Método dos Elementos Finitos:

Elasticidade Plana e Tridimensional

1. Introdução

Este texto resume a aplicação do elemento de deslocamento do Método dos Elementos Finitos na

solução de problemas de elasticidade plana e tridimensional.

Começa-se por definir as hipóteses da análise física e geometricamente linear, e recordar as

situações em que se pode admitir que corpos elásticos têm comportamento plano. É este o tipo de

comportamento utilizado na apresentação e interpretação do Método dos Elementos Finitos, o

qual é posteriormente generalizado para a análise de sólidos elásticos.

Definem-se as variáveis necessárias para caracterizar o problema e as equações que simulam

o comportamento da estrutura, recorrendo-se a um exemplo simples para ilustrar a identificação

das variáveis e a função de cada uma das equações, as quais são organizadas de forma idêntica às

da modelação do comportamento de peças lineares: as condições de equilíbrio e de

compatibilidade, no domínio e na fronteira da peça, e as relações de elasticidade.

O exemplo de aplicação é também usado para mostrar que a solução exacta dessas equações

não tem expressão analítica, na maioria das aplicações. Introduz-se, por isso, a ideia básica da

formulação do Método dos Elementos Finitos aqui adoptada, de aproximar o campo de

deslocamentos. Descreve-se depois, sumariamente, a generalização para o caso plano da

formulação anteriormente desenvolvida para a análise de peças lineares.

Essa generalização é ilustrada com um exemplo de introdução, o qual é utilizado para

formular e analisar os critérios de aproximação do campo de deslocamento e a consequente

aproximação do campo de deformações, de modo a assegurar a condição central de obter

soluções aproximadas cinematicamente admissíveis. A aproximação que se obtém para o campo

de tensões aplicando as relações de elasticidade é utilizada para mostrar que não é possível, em


2

geral, satisfazer localmente as condições de equilíbrio, tanto no domínio como na fronteira do

elemento. O conceito de força nodal equivalente é depois generalizado e utilizado para formular

a equação resolvente do Método dos Elementos Finitos, cuja interpretação é análoga à obtida

para peças lineares: define o equilíbrio nodal dos elementos da malha, combinando as forças

nodais devidas aos deslocamentos nodais (os termos da matriz de rigidez do elemento) e

resultante das forças nodais equivalentes às forças aplicadas, designadamente as forças de massa

e as forças de fronteira do elemento.

Estes conceitos são depois gradualmente generalizados. Começa-se por definir os elementos

mais simples usados na análise de estados planos e usam-se esses resultados para generalizar a

formulação do Método dos Elementos Finitos. As dificuldades de implementação do exemplo

usado para ilustrar essa generalização são depois usadas para justificar a formulação dos

elementos isoparamétricos, correntemente utilizados nos programas comerciais.

A grande vantagem destes elementos é permitirem generalizar, de uma maneira

computacionalmente muito eficaz, a representação da geometria da peça em análise e a

formulação de elementos com diferentes graus de aproximação do campo de deslocamentos.

Para além de demonstrar as vantagens alcançadas com o desenvolvimento de elementos

isoparamétricos, o objectivo principal do texto é o de realçar e fundamentar os cuidados a ter na

sua utilização. Dá-se particular atenção à análise da adequabilidade das malhas de discretização

que suportam a utilização desses elementos na análise de problemas bidimensionais.

Concluída a apresentação e a ilustração dos principais conceitos que suportam a aplicação do

Método dos Elementos Finitos na análise de problemas de elasticidade plana, sistematiza-se a

generalização da formulação para a análise de problemas tridimensionais. Como essa

generalização não exige a introdução de novos conceitos, o texto incide fundamentalmente na

identificação das variáveis e das equações que as relacionam, e na técnica usada para discretizar

a geometria do sólido em análise e para aproximar os campos de deslocamento, deformação e

tensão. Generaliza-se o conceito de elemento isoparamétrico e comenta-se brevemente o

processo utilizado na montagem da equação resolvente do Método dos Elementos Finitos.

Como a sistematização da aplicação do método é análoga à anteriormente descrita na sua

aplicação de estruturas articuladas e reticuladas, opta-se aqui por analisar as diferentes fases que

caracterizam a utilização de um programa de elementos finitos, dando-se especial atenção às

fases que dependem da intervenção directa do utilizador, ou seja a definição dos dados, a

avaliação da qualidade da malha de discretização gerada automaticamente e, principalmente, à

análise crítica dos resultados obtidos.


3

2. Hipóteses

Na extensão da aplicação do Método dos Elementos Finitos para análise de problemas de

elasticidade plana e tridimensional continua-se a admitir que o material é elástico linear

(linearidade física) e que os deslocamentos e as deformações são infinitesimais (linearidade

geométrica). Continua-se também a admitir que são desprezáveis as forças de inércia e de

amortecimento que se possam desenvolver durante o carregamento da estrutura (comportamento

quase-estático).

Relativamente à geometria da peça e a forma como é solicitada, consideram-se dois tipos de

problemas que podem ser adequadamente modelados recorrendo à formulação bidimensional da

Teoria da Elasticidade, designadamente, o comportamento de peças sujeitas a estados planos de

tensão ou a estados planos de deformação.

Uma placa é uma peça laminar plana solicitada no próprio plano, com uma espessura

suficientemente pequena em relação às restantes dimensões características para justificar a

hipótese de serem desprezáveis as componentes correspondentes do tensor das tensões. Esta

hipótese é frequentemente utilizada na análise de paredes resistentes.

Se, pelo contrário, a peça é prismática com uma espessura muito maior que as dimensões da

sua secção transversal, em cujo plano actuam as cargas, torna-se legítimo assumir que as secções

se deformam no próprio plano, sendo desprezáveis as três componentes do tensor das

deformações relativas à dimensão transversal. Esta é a hipótese adoptada, por exemplo, na

análise de muros de suporte e de secções de túneis.

Figura 1: Problemas de Elasticidade Plana e Tridimensional

Como se ilustra na Figura 1, o elemento estrutural é representado no plano em que actuam as

cargas, o qual corresponde ao plano médio da placa, num estado plano de tensão (EPT), ou a

z

y

x

Lx

Ly

Lz

x

y

Ω

Γ

:

:z x y

z x y

L L LEPT

L L LEPD

≈≈

≪

≫

y

z

x

V

: x y zE3D L L L≈ ≈

Γ


4

secção transversal da peça prismática, num estado plano de deformação (EPD). A área ocupada

pela peça define o domínio, Ω , e a linha que a limita define a fronteira, Γ .

Quando as dimensões características da peça são da mesma ordem, ou quando o tipo de acção

torna inadequadas as hipóteses simplificativas sobre o estado de tensão ou de deformação, a peça

é analisada recorrendo à formulação tridimensional da Teoria da Elasticidade. O volume

ocupado pela peça define o domínio, Ω , e a superfície que a limita define a fronteira, Γ .

Para simplificar a apresentação, a aplicação do Método dos Elementos Finitos a problemas de

elasticidade é feita tratando primeiro a análise de estados planos, estabelecendo-se depois a

generalização da formulação para a análise de problemas tridimensionais.

3. Variáveis

Em consequência das hipóteses que caracterizam os estados planos de tensão e de deformação, e

de acordo com a notação definida nas Figuras 1 e 2, o estado de tensão em cada ponto da peça é

definido por três componentes do tensor das tensões (simétrico), as quais são organizadas

matricialmente na forma,

xx

yy

xy

σσσ

=

s (1)

sendo nulas as restantes componentes nos estados planos de tensão: 0xz yz zzσ σ σ= = = .

Figura 2: Convenção adoptada na medição de variáveis

O estado de deformação é definido pelas três componentes correspondentes do tensor das

deformações (simétrico), as quais são organizadas matricialmente de forma análoga,

xx

yy

xy

εεγ

=

e (2)

x

y

,yy yyσ ε

,xy xyσ γ ,xx xxσ ε

x

y

,x xf u

,y yf u

,y yt u

,x xt u


5

em que 2xy xyγ ε= representa a distorção total. As restantes componentes são nulas nos estados

planos de deformação: 0xz yz zzγ γ ε= = = .

Relativamente às forças aplicadas, distinguem-se as forças de domínio ou de massa

(equivalentes às cargas de vão das barras),

x

y

f

f

=

f (3)

e as forças de fronteira (equivalentes às forças de extremidade):

x

y

t

t

=

t

As forças de domínio podem ser distribuídas, na área ou em linha, ou concentradas, e as

forças de fronteira distribuídas ou concentradas. Os deslocamentos correspondentes são

organizados de modo análogo:

x

y

u

u

=

u (4)

Variáveis estáticas Variáveis cinemáticas

Tensões, ( , )ij x yσ Deformações, ( , )ij x yε

Forças, ( , ), ( , )i if x y t x y Deslocamentos, ( , )iu x y

Quadro 1: Variáveis correspondentes

Admite-se, em geral, que não existem momentos aplicados no plano do carregamento, 0zm =

no domínio e na fronteira. A rotação correspondente é uma variável dependente, podendo ser

determinada a partir do campo de deslocamentos:

( )12z y x x yu uω = − ∂ − ∂ (5)

4. Balanço Energético

De acordo com a notação anteriormente definida, são as seguintes as definições do trabalho

realizado pelas forças interiores e exteriores:

T ( + + )i xx xx yy yy xy xyW d dΩ Ω

Ω ε σ ε σ γ σ Ω= =∫ ∫e s

T T ( ) ( )e x x y y x x y yW d d u f u f d u t u t dΩ Γ Ω Γ

Ω Γ Ω Γ= + = + + +∫ ∫ ∫ ∫u f u t (6)

A condição de balanço energético toma agora a seguinte forma:

T T Td d dΩ Ω Γ

Ω Ω Γ= +∫ ∫ ∫e s u f u t (7)


6

Figura 3: Forças descontínuas e deslocamentos correspondentes

A definição (6) mantém-se válida quando se pretende calcular o trabalho realizado por forças

de domínio distribuídas ao longo de linhas e de forças de domínio e de fronteira concentradas,

representando essas forças através de funções de Dirac. Na definição explícita que se obtém,

T T T TeW d d d

Ω Γ ΓΩ Γ Γ= + + +∫ ∫ ∫u f u t u f u F

ℓ

ℓ ℓ

fℓ é a força distribuída aplicada sobre a linha Γ

ℓ e F o vector que lista as forças concentradas.

O vector u define, no primeiro caso, as componentes do deslocamento ao longo da linha Γℓ e,

no segundo, as componentes do deslocamento correspondentes às forças concentradas, como se

ilustra na Figura 3.

5. Equações Básicas

As equações que caracterizam um problema de Elasticidade Plana estão resumidas no Quadro 2,

analisando-se em seguida o significado das condições de domínio (8) a (10) e das condições de

fronteira (11) e (12).

Equilíbrio Elasticidade Compatibilidade

Domínio (8) Domínio (9) Domínio (10)

T emΩ+ =A s f 0 emΩ=s D e emΩ=e Au

TtemΓ=N s t

uemΓ=u u

Fronteira (11) Fronteira (12)

Quadro 2: Equações básicas

Estas equações são a seguir escritas explicitamente, sublinhando-se o significado das

condições físicas que impõem. Dá-se particular atenção às condições de fronteira, pois são

x x

y y

Γℓ

,f uℓ ℓ

,x xF u

,x xF u

,y yF u

,y yF u


7

aquelas que mais directamente prestam informação sobre a qualidade das soluções aproximadas

que são obtidas utilizando o Método dos Elementos Finitos.

5.1 Condições de Domínio

A condição de equilíbrio no domínio assegura que a variação do campo de tensões equilibra, em

cada ponto, as forças de massa,

0

0x xx y xy x

y yy x xy y

fem

f

σ σΩ

σ σ∂ + ∂ + =

∂ + ∂ + = (13)

enquanto a condição de compatibilidade no domínio define, em cada ponto, as medidas de

deformação em função da variação do campo de deslocamentos:

xx x x

yy y y

xy y x x y

u

u em

u u

εε Ωγ

= ∂ = ∂ = ∂ + ∂

(14)

Quando se admite que o material elástico linear é homogéneo e isotrópico, obtém-se a

seguinte definição para as relações constitutivas, isto é, as condições que definem, em cada

ponto, o estado de tensão provocado por um determinado estado de deformação,

2

2

( )1

( )1

2(1 )

xx xx yy

yy yy xx

xy xy

E

Eem

E

σ ε ν εν

σ ε ν ε Ων

σ γν

= + −

= + −

= +

(15)

para estados planos de tensão, sendo E o módulo de elasticidade e ν o coeficiente de Poisson.

A relação constitutiva para estados planos de deformação é a seguinte:

1( )

1 2 1

1( )

1 2 1

2(1 )

xx xx yy

yy yy xx

xy xy

E

Eem

E

νσ ε ν εν ν

νσ ε ν ε Ων ν

σ γν

−= + − + −= + − +

= +

(16)

As condições de equilíbrio, compatibilidade e de elasticidade (13) a (15) ou (16) estão

expressas matricialmente no Quadro 2, de acordo com as definições (1) e (2) para os vectores de

tensão e deformação e as definições (3) e (4) para os vectores das forças de domínio e dos

deslocamentos, respectivamente.


8

As expressões que se encontram para os operadores diferenciais de compatibilidade e de

equilíbrio e para a matriz de rigidez são as seguintes:

0

0x

y

y x

∂ = ∂ ∂ ∂

A

T0

0x y

y x

∂ ∂ = ∂ ∂

A (17)

1 3 0

3 1 01

0 0 1

Gκ κ

κ κκ

κ

+ − = − + − −

D

em que / 2(1 )G E ν= + é o módulo de distorção, com (3 ) /(1 )κ ν ν= − + para estados planos de

tensão e 3 4κ ν= − para estados planos de deformação, respectivamente.

A relação de elasticidade (9) pode ser generalizada para incluir o efeito de campos residuais:

( )r r emΩ= − +s D e e s (18)

As tensões residuais, rs , resultam normalmente de erros de fabrico. Um exemplo típico da

utilização do campo de deformações residuais, re , é a modelação de variações térmicas, sendo,

T 0r T Tα ∆ α ∆=e

para materiais homogéneos e isótropos, em que α é o coeficiente de dilatação térmica e

( , )T x y∆ a variação de temperatura.

Exercício 1: Considere a consola triangular representada na Figura 4. Admitindo ser o

seguinte o campo de deslocamentos,

( , ) 0

( , ) ( / )x

y

u x y

u x y x L d

= =

(19)

com 2,6 / ( )d E m= − , verifique que, para um estado plano de tensão, são as seguintes as

definições dos campos de deformação e de tensão (corte puro),

( , ) ( , ) 0; ( , ) 2,6 /xx yy xyx y x y x y Eε ε γ= = = − (20)

2( , ) ( , ) 0; ( , ) 1 ( )xx yy xyx y x y x y kNmσ σ σ −= = = − (21)

calculadas recorrendo às condições de compatibilidade (14) e de elasticidade (15). Confirme

ainda que a condição de equilíbrio no domínio (13) só é verificada para forças de massa nulas:

0x yf f= = (22)


9

a) Geometria e material b) Deformada

c) Estado de deformação d) Estado de tensão

e) Forças segundo x f) Forças segundo y

Figura 4: Consola triangular com campo de deslocamentos linear

5.2 Condições de Fronteira

As restantes equações resumidas no Quadro 2 definem as condições de fronteira estáticas (11) e

as condições de fronteira cinemáticas (12). Essas condições também são designadas por

condições de Neumann (ou condições de fronteira naturais) e por condições de Dirichlet (ou

condições de fronteira essenciais), respectivamente.

A fronteira da peça, Γ , é decomposta nas duas regiões em que essas condições são impostas,

a fronteira estática tΓ (ou de Neumann) e a fronteira cinemática uΓ (ou de Dirichlet). Essas

regiões são complementares, t uΓ Γ Γ= ∪ e t uΓ Γ∅ = ∩ , por não ser fisicamente possível

impor num ponto simultaneamente uma força e o deslocamento correspondente.

x

y

L

L

2

1

0,3

. ( )

L m

E const kNm

ν−

===

xu

yu

x

y

2,6 / ( )d E m=

( , ) 0

( , ) ( / )x

y

u x y

u x y x L d

==

x

y

2,6 /xy Eγ = − 0xx yyε ε= =

y

x 21( )xy kNmσ −= − 0xx yyσ σ= =

1xt = −

12xt = + 0xt =

0xf = 0yf =

0yt =

1yt = + 12yt = −


10

A condição de fronteira estática (11) define o estado de tensão na fronteira que equilibra as

forças aplicadas. A forma explícita dessa equação é a seguinte,

x xx y xy x

ty yy x xy y

n n tem

n n t

σ σΓ

σ σ+ =

+ = (23)

em que xn e yn representam as componentes da normal exterior unitária, como se ilustra na

Figura 5. Quando esta condição é expressa matricialmente, conclui-se que a matriz de equilíbrio

na fronteira é análoga à matriz de equilíbrio no domínio (17), sendo obtida substituindo a

componente da normal unitária exterior o operador diferencial correspondente:

T0

0x y

y x

n n

n n

=

N (24)

Figura 5: Normal exterior unitária Figura 6: Condições de fronteira

A fronteira estática da placa representada na Figura 6 é definida pelas linhas de contorno ABC

e DEF, devendo aí ser impostas as seguintes condições:

1 12 2

1 12 2

( ) ( ) 0

( ) ( ) 0

xx xy x

AByy xy y

tem

t

σ σΓ

σ σ+ + − = =

− + + = =

(0) ( 1) 0

( 1) (0) 0xx xy x

BCyy xy y

tem

t

σ σΓ

σ σ+ − = =

− + = =

(0) (1) 0

(1) (0)xx xy x

DEFyy xy y

tem

t p

σ σΓ

σ σ+ = =

+ = = −

A condição de compatibilidade na fronteira, a equação (12), assegura que os deslocamentos

medidos no domínio, junto ao contorno, são coerentes com os deslocamentos que aí estejam

impostos. A forma explícita dessa condição é a seguinte,

x

y n

n

xn

xn

yn

yn

( )2 2 1x yn n+ =

A

C

D F

x

y

B

E

p

xxσ

yyσ

xyσ

12xn = +

12yn = −


11

x x

uy y

u uem

u uΓ

= =

ou, para o exemplo da Figura 6:

0

0x x

AFy y

u uem

u uΓ

= = = =

Ilustra-se na mesma figura um terceiro tipo de condição de fronteira, a condição de fronteira

mista ao longo do bordo CDΓ , em que é imposta uma força e a componente complementar do

deslocamento correspondente:

0

(0) (1) 0x x

CDyy xy y

u uem

tΓ

σ σ= =

+ = =

Esta condição de fronteira resulta, normalmente, de simplificações de simetria. A expressão

complementar definiria a condição de fronteira decorrente de uma condição de anti-simetria:

(1) (0) 0

0xx xy x

y y

t

u u

σ σ+ = = = =

Exercício 2: Verifique serem as forças de fronteira representadas na Figura 4 que equilibram

o campo de tensões aí definido. Confirme que essas forças estão globalmente em equilíbrio,

isto é, têm resultantes nulas e produzem um momento nulo em torno do eixo ortogonal ao

plano.

Em vez de impor uma força ou o deslocamento correspondente, é possível impor uma relação

entre a força e o deslocamento. Este tipo de condição de fronteira, também designada por

condição de Robin, é utilizado para simular o comportamento de apoios elásticos, como se

ilustra na Figura 7. A condição é escrita na forma,

e eemΓ=t D u (25)

reunindo a matriz eD os coeficientes de rigidez que caracterizam a deformabilidade do meio de

apoio.

Figura 7: Apoio elástico

eΓ


12

5.3 Fronteiras Interiores

Interessa referir, ainda, as condições de fronteira que são artificialmente introduzidas quando a

peça é discretizada em elementos finitos, como seria, por exemplo, o segmento BE se a consola

representada na Figura 5 fosse discretizada em dois elementos, como se mostra na Figura 8.

Figura 8: Fronteira interior Figura 9: Elemento de junta

As condições nesta fronteira devem garantir a continuidade dos deslocamentos,

j kiemΓ=u u (26)

em que os índices j e k identificam dois elementos que partilham a fronteira interior iΓ , assim

como o equilíbrio das forças em qualquer ponto da fronteira,

j kiemΓ+ =t t t (27)

sendo t a força aplicada exteriormente sobre a fronteira (=t 0 no exemplo da Figura 8).

As forças na fronteira de um elemento são calculadas recorrendo à condição de equilíbrio na

fronteira (11), a qual é válida para qualquer secção do domínio da peça,

j jiem para o elemento jΓ=σσσσN t

em que N é a matriz de equilíbrio (24) escrita para a fronteira iΓ do elemento j e jt o vector

das forças que nessa fronteira equilibram o estado de tensão no elemento.

Ou seja, se se seccionar uma peça, na representação de corpo livre daí resultante as forças

(indeterminadas) que se aplicam na secção de corte são calculadas em função do campo de

tensão recorrendo à equação (11). O conceito é idêntico ao utilizado na definição de diagramas

de corpo livre de peças lineares, em que se aplica na secção de corte como força exterior o

esforço aí existente.

Como as normais exteriores de dois elementos que partilham a mesma fronteira interior são

simétricas, a condição de equilíbrio (27) pode ser expressa directamente em termos do estado de

tensão em cada elemento na forma:

E

y

B

x

1 1,x xt u 2 2,x xt u

2 2,y yt u 1 1,y yt u

B

E

y

x


13

( )j kiemΓ− =σ σσ σσ σσ σN t (28)

É importante notar que, em consequência desta equação, as três componentes do estado de

tensão em cada elemento são relacionadas por apenas duas condições, e não as três que seriam

necessárias para estabelecer o equilíbrio em termos de tensões. Ou seja, permite que o estado de

tensão entre dois elementos seja diferente, j k≠σ σσ σσ σσ σ na fronteira que partilham, iΓ , mesmo que

seja nula a força exterior aí aplicada, =t 0 . Portanto, este enfraquecimento da condição de

continuidade do campo de tensão na peça, j k=σ σσ σσ σσ σ se =t 0 , resulta directamente da

discretização da peça em dois elementos distintos.

Exercício 3: Verifique que a condição (28) escrita para a fronteira entre os dois elementos do

exemplo ilustrado na Figura 7 impõe as condições de continuidade 1 2xx xxσ σ= e 1 2

xy xyσ σ= mas

permite que 1 2yy yyσ σ≠ :

1 2

1 2

1 2

1 0 0 0

0 0 1 0

xx xx

yy yy

xy xy

em BE

σ σσ σσ σ

− − −

=

Exercício 4: Os elementos de junta (unidimensionais, isto é, sem espessura em problemas

planos, como se ilustra na Figura 9) são usados para modelar diferentes condições de contacto

entre elementos. Escreva a condição de fronteira entre dois elementos ligados entre si por um

elemento de junta com as seguintes propriedades: a) Impede o afastamento mas permite o

deslizamento entre elementos, sem atrito; b) Impede o afastamento mas permite o

deslizamento entre elementos, com atrito; c) Permite o movimento elástico entre elementos.

6. Soluções Exactas e Aproximadas

Mantêm-se as classificações anteriormente definidas para os diferentes tipos de solução do

sistema de equações (8) a (12), resumidas no Quadro 2:

• Um campo de tensões, s, que satisfaz as condições de equilíbrio no domínio (8) e na

fronteira (11), incluindo as fronteiras interiores, é, por definição, uma solução

estaticamente admissível;

• Um campo de deslocamentos, u , que é contínuo no domínio da peça e que satisfaz as

condições de fronteira (12), incluindo as fronteiras interiores, é, por definição, uma solução

cinematicamente admissível, sendo o campo de deformações compatível associado, e,

definido pela condição de compatibilidade no domínio (10);


14

• A solução exacta é a solução que para além de ser estática e cinematicamente admissível

satisfaz também a relação de elasticidade (9);

• A solução exacta existe e é única, mas pode não ter expressão analítica.

Para o exemplo da consola triangular definido na Figura 4, o campo de deslocamentos (19)

define uma solução cinematicamente admissível, pois é continuo no domínio da peça, satisfaz a

condição de fronteira cinemática ao longo do bordo encastrado,

( 0) : 0u x yx u uΓ = = = (29)

e está associado a um campo de deformações (20) compatível.

Para o mesmo problema, o campo de tensões (21) define uma solução estaticamente

admissível se se admitir que as forças de massa são nulas, de acordo com a condição (22), e,

ainda, se as forças exteriores aplicadas equilibrarem esse campo de tensões, ou seja, se as

condições de fronteira estáticas forem as seguintes:

1 12 2

( 1) : 1; 0

( ) : ;

t x y

t x y

y t t

y x t t

ΓΓ

= = − = = = + = −

(30)

Nestas condições, a solução é exacta e única, pois os campos de tensão e de deformação

satisfazem também localmente a relação de elasticidade. As reacções no bordo encastrado da

consola são determinadas calculando as forças que aí equilibram o campo de tensão (21):

( 0) : 0; 1u x yx t tΓ = = = (31)

Figura 10: Placa carregada transversalmente Figura 11: Placa carregada axialmente

É muito limitado o número de problemas com relevância prática que têm solução analítica,

mesmo quando a geometria da peça e as condições de fronteira são muito simples. Por exemplo,

a consola carregada transversalmente ilustrada na Figura 10 não tem solução analítica exacta, e o

caso mais simples da carga de tracção uniforme só tem solução analítica exacta quando se

desprezam as forças de massa e o efeito de Poisson.

x

y

a

b

p

x

y

a

b

p


15

Exercício 5: Considere a placa carregada axialmente representada na Figura 11 e verifique o

seguinte, admitindo um estado plano de tensão: a) A solução ( /a)xu x d= , 0yu = é uma

solução cinematicamente admissível; b) A solução xx pσ = , 0yy xyσ σ= = é uma solução

estaticamente admissível quando se desprezam as forças de massa, 0x yf f= = ; c) Essas

soluções definem a solução exacta, com /d a p E= , quando se despreza o efeito de Poisson.

Para melhor esclarecer as razões que, por regra, impedem a determinação de soluções exactas

para problemas de elasticidade plana (e tridimensional), é conveniente eliminar variáveis e

compactar as equações (8) a (12) que governam o problema, tal como resumidas no Quadro 2.

Tal como se fez para as peças lineares, a condição de compatibilidade no domínio (10) é

usada para eliminar as deformações na condição de elasticidade (9), recorrendo-se à expressão

resultante para o campo de tensões,

emΩs = D Au

para exprimir as condições de equilíbrio no domínio (8) e na fronteira (11) em função dos

deslocamentos, e juntando a única condição ainda omissa, a condição de fronteira (12):

T emΩ+ =A D Au f 0 (32)

TtemΓ=N D Au t (33)

uemΓ=u u (34)

Qualquer solução do sistema de equações diferenciais (32) satisfaz localmente (ou de maneira

forte) todas as condições de domínio do problema, designadamente as condições de equilíbrio

(8), de elasticidade (9) e de compatibilidade (10).

A técnica geralmente usada para obter essas soluções consiste em combinar as soluções

complementar e particular,

0c emΩ= +u u u

definindo a primeira o conjunto das soluções da equação homogénea,

Tc emΩ=A D Au 0

e a segunda uma qualquer solução que represente o efeito das forças de massa:

T0 emΩ+ =A D Au f 0

Estas soluções da equação de Navier (32) foram há muito estabelecidas para materiais

homogéneos e isotrópicos, recorrendo a potenciais de deslocamento (funções cujas derivadas

definem o campo de deslocamento, u) ou a potenciais de tensão (funções cujas derivadas


16

definem o campo de tensão, s). Existem também soluções para materiais homogéneos

anisótropos, sendo mais limitados os casos resolvidos para materiais heterogéneos.

Na maioria das aplicações, a dificuldade de obter soluções exactas para problemas de

elasticidade plana (e tridimensional) não está, portanto, na determinação das soluções da equação

de Navier (32). A dificuldade está em conseguir que essas soluções satisfaçam, também de

maneira forte, as condições de fronteira do problema, definidas pelas equações (33) e (34).

7. Método dos Elementos Finitos

O modelo de deslocamento do Método dos Elementos Finitos é a seguir utilizado para

determinar soluções aproximadas do problema de elasticidade plana definido pelas equações

resumidas no Quadro 2, ou pelo sistema equivalente definido pelas equações (32) a (34). A

aplicação do método é conceptualmente idêntica à adoptada na análise de estruturas reticuladas,

podendo ser resumida em quatro fases:

Primeira fase: Aproximações

• A estrutura é discretizada em elementos finitos;

• O campo de deslocamentos é aproximado em cada elemento usando funções contínuas, na

forma,

=ΨΨΨΨu d (35)

sendo ΨΨΨΨ a matriz que reúne as funções de aproximação e d o vector dos deslocamentos

nodais do elemento;

• A condição de compatibilidade (10) é imposta localmente para definir a aproximação do

campo de deformações em cada elemento:

=e B d (36)

= ΨΨΨΨB A (37)

• A relação de elasticidade (9) é imposta localmente para definir o campo de tensões em

cada elemento:

=s DB d (38)

Segunda fase: Forças Nodais Equivalentes

• A equação resolvente do elemento é estabelecida interpretando-a como uma condição de

equilíbrio de forças nodais equivalentes,

=K d F (39)

reunindo a matriz de rigidez do elemento as forças devidas aos deslocamentos nodais,


17

T dΩ

Ω= ∫K B D B (40)

e o vector,

f t n= + +F F F F (41)

as forças nodais equivalentes às forças de massa,

Tf d

ΩΩ= ∫ ΨΨΨΨF f (42)

as forças nodais equivalentes às forças aplicadas na fronteira,

Tt d

ΓΓ= ∫ ΨΨΨΨF t (43)

e as forças concentradas aplicadas nos nós do elemento, nF .

Terceira fase: Equação Resolvente

• A solução aproximada em cada elemento está sujeita à condição de ser cinematicamente

admissível, sendo construída de maneira a satisfazer localmente (ou de maneira forte) a

condição de fronteira cinemática (34) da malha de elementos finitos e, ainda, a condição

equivalente (26) entre elementos;

• Essa condição é imposta relacionando os deslocamentos nodais dos elementos, d, com os

deslocamentos nodais da malha de elementos finitos, q, através de uma condição de

incidência nodal:

= ℑℑℑℑd q (44)

• As forças nodais equivalentes na malha de elementos finitos definem as resultantes das

contribuições das forças nodais equivalentes geradas em cada elemento,

T= ℑℑℑℑQ F (45)

e são utilizadas para estabelecer a equação resolvente do problema:

* =K q Q (46)

• Este sistema de equações define as condições de equilíbrio das forças nodais equivalentes,

impondo aproximadamente (ou de maneira fraca) as condições de equilíbrio no domínio

(8) e na fronteira estática (11), ou as equações equivalentes (32) e (33), estendidas de modo

a incluírem as condições de equilíbrio (28) nas fronteiras entre elementos.

Quarta fase: Análise da Solução

• Resolvido o sistema (46) nas incógnitas do problema, os deslocamentos nodais da malha

de elementos finitos, q, os deslocamentos, as deformações e as tensões são calculadas em

cada elemento recorrendo às aproximações (35), (36) e (38), depois de determinar os

deslocamentos nodais em cada elemento através da relação de incidência (44).


18

Este procedimento é a seguir apresentado e interpretado recorrendo a um exemplo simples,

sendo posteriormente sistematizado. Discute-se a definição das funções que são utilizadas para

estabelecer a aproximação em que o método se baseia, a aproximação (35) do campo de

deslocamento em cada elemento. Essa generalização da aproximação é depois utilizada para

ilustrar as operações (44) e (45) de reunião de elementos de que decorre a formulação da

equação resolvente (46) para uma malha de elementos finitos.

Exercício 6: Deduza a equação resolvente elementar (39) e obtenha as definições (40) a (43)

impondo as aproximações (35), (36) (38) na condição de balanço energético (7). Recupere os

mesmos resultados estacionarizando a energia potencial do elemento:

T T T T12Min d d d

Ω Ω ΓΠ Ω Ω Γ= − − −∫ ∫ ∫e s u f u t u F

Exercício 7: Generalize a expressão da equação resolvente (39) para incluir a definição do

vector das forças nodais equivalentes a campos de deformação e de tensão residuais (18):

T ( )r r r dΩ

Ω= −∫F B s De (47)

Exercício 8: Mostre ser a seguinte a contribuição de uma fronteira elástica (25) para a matriz

de rigidez elementar:

T

ee e ed

ΓΓ= ∫ Ψ ΨΨ ΨΨ ΨΨ ΨK D (48)

8. Exemplo de Introdução

O exemplo de aplicação escolhido é o problema de estado plano de tensão representado na

Figura 4, em que se despreza o efeito das forças de massa ( 3kNm− ), como estabelece a condição

(22), e se considera o carregamento indicado na Figura 12a).

De acordo com a identificação das fronteiras definida na Figura 12b), a placa triangular está

encastrada no bordo 0x = , sujeita a uma carga uniforme transversal no bordo 1y = , definindo o

lado y x= um bordo livre. A condição de fronteira cinemática (29) mantém-se válida, sendo a

condição de fronteira estática (30) substituída pela seguinte:

( 1) : 0; 1

( ) : 0t x y

t x y

y t t

y x t t

ΓΓ

= = = − = = =

(49)

Este exemplo é resolvido admitindo a mais simples aproximação possível. A placa é

discretizada num único elemento em que se admite ser linear a variação do campo de

deslocamentos, descrita pelos deslocamentos do vértice livre da consola, os deslocamentos

nodais 1d e 2d identificados na Figura 12c).


19

a) Carregamento b) Fronteiras

c) Aproximação dos deslocamentos d) Aproximação das tensões

Figura 12: Consola triangular sujeita a uma carga transversal uniforme

Como se usa apenas um elemento, este exemplo não serve para ilustrar as operações de

reunião de elementos para formar a malha de discretização, definidas pelas condições de

incidência nodal (44) e (45), as quais são analisadas posteriormente. O que se pretende ilustrar

são os dois aspectos fundamentais da formulação de deslocamento do Método dos Elementos

Finitos aplicado à solução de problemas de elasticidade plana:

• A solução aproximada é construída de maneira a ser cinematicamente admissível, sendo

escrita em função dos deslocamentos nodais, os deslocamentos 1d e 2d identificados na

Figura 12c);

• A aproximação resultante para o campo de tensões, definida na Figura 12d), viola, em

geral, as condições de equilíbrio, sendo o equilíbrio imposto sobre forças nodais

equivalentes às forças de massa e de fronteira, as forças 1F e 2F indicadas na mesma

figura.

8.1 Aproximação dos Deslocamentos

A hipótese do campo de deslocamentos variar linearmente no domínio do elemento é formulada

da seguinte maneira ( 1 )L m= ,

x

y

1m

1m 0,3

.E const

ν ==

21 ( )p kNm−=

( 1)t yΓ =

( 0)u xΓ =

( )t y xΓ =

x

y

1m

1m

1d

2d

xu

yu

1

2

( )

( )x

y

u x d

u x d

= =

x

y

1m

1m

xxσ

yyσ

xyσ

2F

1F

1

1

2

(1,099 )

(0,330 )

(0,385 )

xx

yy

xy

E d

E d

E d

σσσ

= = =


20

0 1 2

0 1 2

( , ) ( / ) ( / )

( , ) ( / ) ( / )x

y

u x y a x L a y L a

u x y b x L b y L b

= + += + +

sendo importante notar que apenas as componentes de translação do movimento são

aproximadas independentemente. A componente de rotação é uma variável dependente. É

calculada pela equação (5), encontrando-se:

11 22 ( )z b aω = −

Esta aproximação permite representar os três movimentos de corpo rígido no plano,

designadamente duas translações, 0 0a ≠ com 1 2 0a a= = e 0 0b ≠ com 1 2 0b b= = , e uma

rotação, 2 1 0a b≠ ≠ , e, ainda, a mais simples aproximação para o campo de deformações, o

estado de deformação constante, que adiante se define.

Os pesos das funções de aproximação, ia e ib , são escolhidos usando critérios análogos aos

adoptados na formulação dos elementos finitos para peças lineares. São definidos impondo

localmente a condição de fronteira cinemática (29),

0 2 0 2 0a a b b= = = =

a qual elimina as componentes de corpo rígido, e escolhidos de modo a representar as

componentes do deslocamento num nó do elemento, neste caso o vértice livre da placa,

1 1 1

2 1 2

(1,1)

(1,1)x

y

d u a d

d u b d

= ⇒ == ⇒ =

ver Figura 12c), obtendo-se a seguinte expressão para a aproximação (35) do campo de

deslocamentos:

1

2

0

0x

y

u dx

u dx

=

(50)

Sendo o campo de rotação dependente das componentes de translação, as rotações dos nós

não são escolhidas como incógnitas do problema, sendo o vector dos deslocamentos nodais

presente na definição (35) sempre definido em termos das componentes de translação medidas

no referencial global da estrutura.

8.2 Aproximação das Deformações e das Tensões

Em consequência das condições anteriormente impostas sobre a aproximação dos deslocamentos

(50), para assegurar que essa aproximação define uma solução cinematicamente admissível no

sentido forte, isto é, em todos os pontos do domínio da placa e da sua fronteira cinemática, basta

definir a aproximação do campo de deformações impondo localmente a condição de

compatibilidade no domínio (14), encontrando-se o seguinte resultado para a definição (36):


21

1

2

1 0

0 0

0 1

xx

yy

xy

d

d

εεγ

=

(51)

Como o campo de deformações é constante, a aproximação (38) que resulta para o campo de

tensões impondo localmente as relações de elasticidade (15) é:

1

2

1,099 0

0,330 0

0 0,385

xx

yy

xy

Ed

Ed

E

σσσ

=

(52)

8.3 Condições de Equilíbrio no Domínio e na Fronteira

Como a aproximação linear do campo de deslocamentos está associada a uma aproximação

constante (52) para o campo de tensões e se admite que as forças de massa (22) são nulas,

conclui-se que (contra o que geralmente sucede) a condição de equilíbrio (13) é localmente

satisfeita neste exemplo de aplicação.

a) Forças segundo x b) Forças segundo y

Figura 13: Forças de massa e de fronteira equilibradas pela aproximação do campo de tensões

No entanto, quando se recorre à definição (23) para calcular as forças na fronteira que

equilibram o estado de tensão,

2

1

1

(0,385 )0 0 1( 1) :

(0,330 )0 1 0

xxx

yyy

xy y

t E dy

t E d

σΓ σ

σ=

= = =

1 12 2 1 2

1 11 22 2

0 (0,777 ) (0,272 )( ) :

(0,233 ) (0,272 )0

xxx

yyy

xy y x

t E d E dy x

t E d E d

σΓ σ

σ

+ −

− +

=

+ − = = = − +

2(0,385 )E d

1(1,099 )E d

1

2

(0,777 )

(0,272 )

E d

E d

+−

0xf =

1(0,330 )E d

2(0,385 )E d

1

2

(0,233 )

(0,272 )

E d

E d

−+

0yf =


22

1

2

0

(1,099 )1 0 0( 0) :

(0,385 )0 0 1

xxx

yyy

xy x

t E dx

t E d

σΓ σ

σ=

−− = = = −−

conclui-se que não é possível escolher valores para os graus de liberdade do problema, os

deslocamentos nodais 1d e 2d , que permitam impor localmente as condições de fronteira

estáticas (49), como se ilustra na Figura 13.

Para obter uma solução aproximada útil torna-se necessário impor as condições de equilíbrio

de uma maneira fraca, recorrendo-se novamente ao critério de determinar os deslocamentos

nodais equacionando o trabalho das forças interiores ao das forças exteriores.

8.4 Equação Resolvente

A equação (39) que se obtém aplicando esse critério (ver Exercício 6) tem a seguinte expressão:

1

2

0,5(1,099 ) 0 0

0 0,5(0,385 ) 0,5

dE

dE

= −

(53)

A matriz de rigidez é calculada aplicando a definição (40), explorando o facto de serem

constantes os campos de tensão e de deformação no domínio da consola triangular (com área

20,5mΩ = ):

1 1T T T

01 0,5

xd dy dx

ΩΩ= = =∫ ∫ ∫K B D B B D B B D B

11 12

21 22

0,5(1,099 ) 0

0 0,5(0,385 )

K K E

K K E

= =

K (54)

Como as forças de massa são nulas, =f 0 em Ω , e não existem forças aplicadas no nó livre

da consola, só a força transversal (49) contribui para o vector das forças nodais equivalentes (41)

obtendo-se o seguinte resultado aplicando as definições (42) e (43):

1

1

0

0,5

F

F

= −

F = (55)

T

1 1

0

00 0

00 0x

f xy

fxdy dx

fx

= = =

∫ ∫F =

T T

1 1

0 0

0 00 0 0

1 00 0 0,5x x

ty yy = 1 y = x

t tx xdx dx

t tx x

= = + = = − = −

∫ ∫F =


23

8.5 Interpretação da Equação Resolvente

Na interpretação dada ao sistema resolvente (53) no contexto da aplicação do Método dos

Deslocamentos na análise de estruturas reticuladas, cada equação define uma condição de

equilíbrio entre dois conjuntos de forças nodais, as forças nodais devidas aos deslocamentos

nodais iq e as forças nodais aplicadas, iF .

a) Acção de 1 1d = b) Acção de 2 1d = c) Acção do carregamento

Figura 14: Interpretação dos termos da equação resolvente

Como se ilustra na Figura 14, continua a ser possível interpretar o coeficiente ijK da matriz

de rigidez K como a força nodal iF devida ao deslocamento nodal 1jq = e, no termo

independente do sistema (53), os coeficientes do vector de forças nodais F como sendo as

forças nodais iF aplicadas. A diferença essencial é que essas forças nodais representam agora

forças concentradas equivalentes às que equilibram o campo de tensões no domínio e na

fronteira do domínio bidimensional agora em análise.

Nas Figura 15a) e 15c) definem-se os sistemas de forças que equilibram a aproximação das

tensões devidos aos deslocamentos unitários, de acordo com os resultados apresentados na

Figura 13, e nas Figuras 15b) e 15d) definem-se as forças nodais equivalentes.

Sendo nulas as forças de massa e linear a aproximação do campo de deslocamentos, essas

forças são calculadas, como adiante se mostra, determinando a resultante das forças de fronteira

e atribuindo metade dessa resultante a cada nó de extremidade.

1 1d =

2 1d =

1m

1m

1p =

11K

21K

12K

22K

1F

2F


24

a) Forças de domínio e de fronteira devidas a 1 1d =

b) Forças nodais equivalentes devidas a 1 1d =

c) Forças de domínio e de fronteira devidas a 2 1d =

d) Forças nodais equivalentes devidas a 2 1d =

Figura 15: Forças nodais equivalentes devidas à aproximação do deslocamento

0,777xt E= +

0xt =

0xf = 1,099xt E= −

0yf = 0yt =

0,330yt E= +

0,233yt E= −

0

12 ( 1,099 )E−

0

22

12

( 0,777 )

( 1,099 )

E

E

+

= +

0

0

12 ( 0,330 )E+

22

12

( 0,233 )

( 0,330 )

E

E

−

= −

0xf =

0,385xt E= +

0,272xt E= −

0xt = 0yf =

0yt =

0,272yt E= +

0,385yt E= −

0

0

12 ( 0,385 )E+

22

12

( 0,272 )

( 0,385 )

E

E

−

= −

0 0

12 ( 0,385 )E−

22

12

( 0,272 )

( 0,385 )

E

E

+

= +


25

Os resultados das Figuras 15b) e 15d) recuperam a definição (54) da matriz de rigidez,

mostrando que:

• As forças nodais equivalentes 1F e 2F devidas ao deslocamento 1 1d = são

211 20 (0,777 ) 0,5(1,099 )K E E= + = e 21

21 2 2( 0,330 ) ( 0,233 ) 0K E E= + + − = ;

• As mesmas forças devidas ao deslocamento 2 1d = são, respectivamente,

2112 2 2( 0,330 ) ( 0,233 ) 0K E E= − + + = e 2

22 20 (0,272 ) 0,5(0,385 )K E E= + = .

Quando se aplica o mesmo critério para determinar as forças nodais equivalentes ao

carregamento da consola obtém-se a definição (55), como se ilustra na Figura 16.

a) Carregamento b) Forças equivalentes

Figura 16: Forças nodais equivalentes ao carregamento

8.6 Análise da Solução

Depois de calcular a solução do sistema (53), 1 0d = e 2 2,6 / ( )d E m= − , aplicam-se as

definições (50) a (52) para determinar as aproximações dos deslocamentos, deformações e

tensões. A solução que assim se obtém está representada na Figura 4, concluindo-se que:

• A solução é cinematicamente admissível e produz uma aproximação coerente para a

deformada, ilustrada na Figura 4b);

• Como a aproximação dos deslocamentos é linear e o material homogéneo e isótropo,

obtém-se uma aproximação constante para o campo de tensões, definida na Figura 4d);

• Essa aproximação verifica (neste caso particular) a condição de equilíbrio no domínio (por

serem nulas as forças de massa), mas viola as condições de equilíbrio na fronteira;

• O carregamento dado, representado na Figura 12a) e definido pelas condições de fronteira

(49), não é recuperado pelo sistema de forças que equilibra na fronteira estática a

aproximação do campo de tensão, definido nas Figuras 4e) e 4f) e que corresponde ao

carregamento definido na Figura 17;

1m

21p kNm−=

10,5kNm−


26

• No entanto, esses dois sistemas de forças têm as mesmas forças nodais equivalentes, como

se ilustra na Figura 18, por ser apenas isso o que se impõe ao estabelecer a equação

resolvente (39) do Método dos Elementos Finitos.

Figura 17: Forças na fronteira em equilíbrio com a aproximação das tensões

Figura 18: Forças nodais equivalentes às forças dadas e às forças calculadas

8.7 Refinamento e Convergência

Tal como se fez na análise de estruturas reticuladas, as soluções obtidas para problemas de

elasticidade plana podem ser melhoradas mantendo o grau de aproximação e aumentando o

número de elementos de discretização da peça (refinamento-h), ou mantendo a malha de

discretização e aumentando o grau da aproximação (refinamento-p).

Na Figura 19 apresentam-se as soluções obtidas para o problema da consola quando se

subdivide a malha e se mantém o grau de aproximação (linear) em cada elemento. A solução

mais refinada aí apresentada, em que se nota a localização da concentração de tensões nos cantos

da secção de encastramento, está já muito próxima da solução exacta.

A análise desses resultados permite estabelecer conclusões análogas às obtidas na análise de

peças lineares:

• A convergência em deslocamentos é relativamente rápida, sendo muito mais lenta a

convergência em tensões (a convergência para a função que se aproxima directamente é

mais rápida do que a convergência para a sua derivada);

• A solução é localmente compatível, no domínio e na fronteira, mas mais rígida do que a

solução exacta (subestima o deslocamento máximo);

• A solução é localmente desequilibrada, no domínio e na fronteira, e não está do lado da

segurança (subestima a tensão máxima);

• A convergência para a solução exacta é ainda fraca sempre que se reconhece a malha de

elementos finitos na representação dos campos de tensão (desequilíbrio entre elementos).

1m

1m

21kNm−

21kNm−

12 0,5F kNm−= −

1 0F =


27

1 elemento 4 elementos 384 elementos 15000 elementos

Figura 19: Refinamento-h da solução da consola triangular (elementos lineares)

A segunda opção de refinamento (refinamento-p) está ilustrada nas Figuras 20 e 21 para o

exemplo da consola quadrada ( )a b= representada na Figura 10. O problema é resolvido com

elementos de 3 nós (aproximação linear dos deslocamentos) e 6 nós (aproximação quadrática),

os quais são adiante definidos. As duas malhas têm o mesmo número de nós livres, pelo que os

sistemas resolventes (46) correspondentes envolvem o mesmo número de incógnitas (dois

deslocamentos por nó livre).


28

512 elementos lineares 128 elementos quadráticos

Figura 20: Solução da consola quadrada com elementos lineares e quadráticos (289 nós)

512 elementos lineares 128 elementos quadráticos

Figura 21: Solução da consola quadrada com elementos lineares e quadráticos (289 nós)


29

Ambas as soluções são cinematicamente admissíveis mas qualquer delas ainda está longe de

aproximar correctamente o campo de tensão, pois é ainda visível o desequilíbrio entre elementos.

No entanto estes resultados servem para ilustrar outra conclusão tirada na análise de estruturas

reticuladas pelo Método dos Elementos Finitos:

• Para o mesmo número de graus de liberdade, a qualidade da solução cresce com o grau da

aproximação do campo de deslocamentos.

9. Discretização e Aproximação

Na aplicação do Método dos Elementos Finitos à solução de problemas de Elasticidade Plana o

domínio é decomposto em elementos com geometria simples de modo a viabilizar e facilitar as

várias fases de aplicação do método resumidas na Secção 7.

A utilização de elementos simples, tipicamente elementos com três e quatro lados, assegura a

primeira condição de facilitar a representação, com o rigor desejado, da geometria e das

condições de fronteira do domínio em análise. Na Figura 22 ilustra-se a discretização de um

açude usando elementos triangulares de 3 nós. A modelação da galeria no centro do açude pode

ser melhorada diminuindo a dimensão dos elementos e/ou aumentando o número de nós do

elemento, como adiante se mostra.

Figura 22: Discretização de um açude com elementos triangulares de 3 nós

Uma segunda condição importante é ser fácil definir expressões gerais para as funções usadas

na aproximação (35) do campo de deslocamentos. A terceira condição, essencial à formulação

em termos de deslocamentos do Método dos Elementos Finitos, que aqui se adopta, é assegurar

que essas funções produzam soluções cinematicamente admissíveis.

A última condição, muito importante em termos de implementação numérica, é a de facilitar o

cálculo dos integrais que definem a matriz de rigidez (40) de cada elemento e os vectores das

forças nodais equivalentes às forças de massa (41) e de fronteira (43).

2

1

2m

1 1

0,5

1 1 1 1

0,3

0,4

0,2

1

bd

ad dd

y

x 2

cd


30

Como adiante se mostra, a maneira mais eficaz de satisfazer esses objectivos é recorrer ao

conceito de elemento isoparamétrico. No entanto, nesta fase da apresentação é suficiente analisar

como se pode aproximar uma função sobre os dois tipos de elementos acima referidos, usando a

mais simples aproximação que cada um deles permite, escrita na forma,

1

( , ) ( , )N

i ii

f x y x y fΨ=∑= (56)

designadamente o elemento triangular de 3 nós ( 3)N = e o elemento rectangular de 4 nós

( 4)N = representados nas Figura 23 e 24.

Figura 23: Elemento triangular de 3 nós Figura 24: Elemento rectangular de 4 nós

Mantém-se o princípio de basear a aproximação em funções polinomiais, definindo-as de

modo a terem valor unitário num nó e nulo nos restantes,

1

( , )0i j j

se j ix y

se j iΨ

= ≠

= (57)

para assegurar que os coeficientes da aproximação, if , representam o valor da função nos nós do

elemento. Para além disso, é necessário assegurar que a aproximação (56) permite representar a

função polinomial mais simples, a função unitária:

1

( , ) 1N

ii

x yΨ=

=∑ (58)

9.1 Aproximação Linear

Em consequência das condições acima impostas, as funções presentes na aproximação (56) para

o elemento de 3 nós definem os planos representados na Tabela 1, com a seguinte expressão

geral,

1

( , ) ( ) 1, 2, 32i i i ix y x y com i

AΨ α β γ= + + =

em que A representa a área do elemento triangular, dada por,

x

y

1 1( , )x y

2 2( , )x y

3 3( , )x y

x

y

1 1( , )x y 2 2( , )x y

3 3( , )x y

4 4( , )x y

a

b


31

1 2 32A α α α= + +

se a numeração dos nós for sequencial e no sentido inverso ao dos ponteiros do relógio, tendo os

termos constantes as seguintes expressões,

i j k k jx y x yα = −

i j ky yβ = −

i k jx xγ = −

em que i j k≠ ≠ e onde se permutam os índices pela sequência natural:

1 2 3 3 2 2 3 1 1 3 3 1 2 2 1; ;x y x y x y x y x y x yα α α= − = − = −

1 1 1 1

1 1

1 1

( ) /2

/2

/2x

y

x y A

A

A

Ψ α β γΨ βΨ γ

= + +∂ =∂ =

2 2 2 2

2 2

2 2

( ) /2

/2

/2x

y

x y A

A

A


= + +∂ =∂ =

3 3 3 3

3 3

3 3

( ) /2

/2

/2x

y

x y A

A

A


= + +∂ =∂ =

Tabela 1: Elemento triangular de 3 nós

Estes resultados mostram que a aproximação (56) para 3N = é linear e completa, isto é,

envolve todos os termos lineares, x e y , e todos os monómios de grau inferior, que agora se

reduz ao termo constante.

Exercício 9: Recupere a aproximação (50) do campo de deslocamentos na consola triangular

representada na Figura 12 usando as funções de aproximação definidas na Tabela 1.

9.2 Aproximação Bilinear

Quando o mesmo processo de construção das funções de aproximação é aplicado ao domínio

rectangular representado na Figura 24, trabalhando apenas com os nós colocados nos vértices,

obtêm-se as funções definidas e representadas na Tabela 2, em que se usa a notação 1'x x x= −

1'y y y= − , para simplificar a definição das funções de aproximação.

1 1( , )x y

2 2( , )x y

3 3( , )x y

1

1 1( , )x y

2 2( , )x y

3 3( , )x y

1

1 1( , )x y

2 2( , )x y

1

3 3( , )x y


32

A análise dessas funções mostra que a aproximação (56) para 4N = é quadrática mas

incompleta, isto é, envolve apenas o termo bilinear, xy (ou ' 'x y na notação usada), dos três

monómios quadráticos possíveis, designadamente 2x , xy e 2y .

1

1

1

(1 ' / ) (1 ' / )

(1 ' / ) /

(1 ' / ) /x

y

x a y b

y b a

x a b

ΨΨΨ

= − −∂ = − −∂ = − −

2

2

2

(1 ' / ) ' /

(1 ' / ) /

' /x

y

y b x a

y b a

x ab

ΨΨΨ

= −∂ = −∂ = −

3

3

3

' ' /

' /

' /x

y

x y ab

y ab

x ab

ΨΨΨ

=∂ =∂ =

4

4

4

(1 ' / ) ' /

' /

(1 ' / ) /x

y

x a y b

y ab

x a b

ΨΨΨ

= −∂ = −∂ = −

Tabela 2: Elemento rectangular de 4 nós 1 1( ' ; ' )x x x y y y= − = −

9.3 Grau da Aproximação no Domínio e na Fronteira

A aproximação (56) é adiante generalizada para obter uma maior precisão na representação da

função, tanto para elementos triangulares como para elementos quadrangulares.

Esse enriquecimento é obtido aumentando o número de termos da aproximação, N , ou seja, o

número de nós do elemento de acordo com a condição (57), garantindo que se incluem na

aproximação tantos termos quantos possíveis de grau inferior ao termo de maior grau, M , na

aproximação (56).

1 1( , )x y

2 2( , )x y

3 3( , )x y

4 4( , )x y

1

1 1( , )x y

2 2( , )x y

3 3( , )x y

4 4( , )x y

1

1 1( , )x y

2 2( , )x y

3 3( , )x y

4 4( , )x y

1

1 1( , )x y

2 2( , )x y

3 3( , )x y

4 4( , )x y

1


33

O problema que se põe ao escolher um elemento com N nós para realizar uma dada análise é

saber qual é o grau, M , que se está a usar na aproximação no domínio do elemento, e o grau da

aproximação na fronteira, fM .

A maneira mais expedita para obter esta informação é recorrer ao Triângulo de Pascal

representado na Figura 25, o qual é construído de maneira a agrupar em cada linha todos os

monómios possíveis de um dado grau. Os lados do triângulo com origem no termo constante

definem o grau da aproximação (56) nos lados do elemento. Os termos interiores identificam o

grau dos monómios da aproximação (56) presentes no domínio do elemento.

2 2

3 2 2 3

4 3 2 2 3 4

1

x y

x xy y

x x y x y y

x x y x y xy y

Constante

Linear

Quadrático

Cúbico

Quártico

Figura 25: Triângulo de Pascal

Se se admitir que a aproximação na fronteira do elemento é completa, o que se verifica tanto

para elementos triangulares como quadrangulares, conclui-se que o grau de aproximação num

lado com fN nós é 1f fM N= − . É o que se verifica para o elemento triangular de 3 nós e para

o elemento rectangular de 4 nós anteriormente definidos: a aproximação é linear nos lados de

ambos os elementos ( 2fN = , 1fM = ), como se ilustra nas Tabelas 1 e 2.

Nas Figuras 26 e 27 definem-se outros tipos de elementos triangulares (de 6 e 10 nós) e

quadrangulares (de 8 e 9 nós). A todos eles se aplica a definição anteriormente dada para o grau

da aproximação (56) nos lados dos elementos: quadrática nos lados com 3 nós e cúbica nos lados

com 4 nós.

Figura 26: Elementos triangulares com 3, 6 e 10 nós.

1

1

1 2

3

ξ

η

1

1

1 2

3 ξ

η

4

5

6

1

1

1 2 3

ξ

η

4

5

6

7

8

9 10


34

Figura 27: Elementos quadrangulares com 4, 8 e 9 nós.

Na literatura sobre o Método dos Elementos Finitos designam-se por serendipianos os

elementos sem nós interiores e Lagrangianos os elementos que os têm. Essa distinção decorre da

escolha das funções de aproximação, sabendo-se que, para a mesma ordem de aproximação, os

elementos Lagrangianos têm melhor comportamento.

Se se admitir que a aproximação (56) num elemento triangular com N nós é completa, para

determinar o grau no domínio do elemento, M , basta varrer o Triângulo de Pascal por linhas e

contabilizar tantos monómios quanto o número de nós do elemento.

Aplicando este processo, conclui-se que a aproximação no elemento de 3 nós é linear, 1M = ,

como anteriormente se verificou, sendo quadrática, 2M = , e cúbica, 3M = para os elementos

triangulares com 6 e 10 nós representados na Figura 25.

O processo para determinar o grau da aproximação (56) num elemento quadrangular com N

nós é semelhante mas dificultado pelo facto da aproximação ser incompleta, isto é, não envolver

todos os monómios de grau M , ou mesmo 1M − , sendo M o maior grau envolvido na

aproximação.

Já se verificou que o rectângulo de 4 nós é completo nos termos lineares mas só contém o

termo bilinear dos monómios quadráticos do Triângulo de Pascal. O elemento de 8 nós

representado na Figura 26 contém todos os termos quadráticos, mas apenas dois dos monómios

cúbicos, os termos 2x y e 2xy . Para manter a simetria da aproximação, isto é, para não dar maior

peso aos monómios em x ou em y , o nó adicional no elemento de 9 nós é usado para introduzir

o termo 2 2x y , sendo portanto a aproximação incompleta nos termos cúbicos e quárticos.

Estes elementos e estes critérios de definição das funções de aproximação serão adiante

utilizados para definir os elementos isoparamétricos usados na solução de problemas de

Elasticidade Plana.

1 2

3

ξ

η

4

1 1

1

1

1

2

3

ξ

η

4

5

6

7

8

1 1

1

1

1

2

3

ξ

η

4

5

6

7

8 9

1 1

1

1


35

10. Generalização da Aproximação

As funções de aproximação anteriormente definidas podem ser utilizadas para interpolar uma

qualquer função definida num domínio plano. Por exemplo, para representar a topografia de uma

região, pode-se triangularizar o terreno, definir as cotas nos nós da malha de discretização e

aproximar o relevo usando a aproximação (56) em cada um dos elementos, com 3N = ,

representando if a cota do ponto do terreno que coincide com esse nó.

A mesma ideia é aplicada à aproximação do campo de deslocamentos num elemento finito

plano, em que os coeficientes if da aproximação (56) passam a representar as componentes do

deslocamento nos nós de cada elemento em que a estrutura foi discretizada. Tomando como base

essa aproximação do campo de deslocamentos, pode-se determinar as aproximações

correspondentes para os campos de deformação e de tensão. É esta primeira fase da descrição

sumária do Método dos Elementos Finitos apresentada na Secção 7 que a seguir se analisa:

• A discretização da estrutura em elementos finitos;

• A aproximação do campo de deslocamentos em cada elemento na forma (35);

• A determinação dos campos de deformação (36) e de tensão (38) correspondentes;

• A verificação das condições de continuidade desses campos quando os elementos são

reunidos para constituir a estrutura.

Figura 28: Dados e discretização para o exemplo da placa trapezoidal

10.1 Aproximação do Campo de Deslocamento

Independentemente da geometria do elemento, as componentes (de translação) do campo de

deslocamento são sempre medidas no referencial global da malha em que o elemento se insere.

Para além disso, a aproximação (56) é aplicada independentemente a cada componente do

deslocamento, sendo também os deslocamentos nodais medidos no referencial global, de acordo

com a sequência de numeração dos nós, como se mostra na Figura 29:

1 1

2m

x

y

10

20 kPa

3

0,3

x

y

f 0

f 20 kNm

E 200 GPa

ν

−

=

= −

==


36

1

1

( , ) ( , )

( , ) ( , )

N

x i ixi

N

y i iyi

u x y x y d

u x y x y d

Ψ

Ψ

=

=

=

=

∑

∑ (59)

De acordo com esta aproximação, quando se provoca um deslocamento no nó j segundo a

direcção x, jxd , e se impedem todos os outros deslocamentos nodais, 0jyd = e 0ix iyd d= = para

i j≠ , o elemento deforma-se apenas com deslocamentos segundo x,

( , ) ( , ) ; ( , ) 0x j jx yu x y x y d u x yΨ= = (60)

obtendo-se um resultado análogo, mas segundo y, quando se impõe um deslocamento nodal

nessa direcção:

( , ) 0; ( , ) ( , )x y j jyu x y u x y x y dΨ= =

a) Elemento triangular de 3 nós b) Elemento rectangular de 4 nós

Figura 29: Convenção para a medição das componentes do campo de deslocamento

A expressão matricial da aproximação (59) é escrita na forma que for mais adequada para a

questão que se estiver a tratar. Pode ser conveniente definir separadamente cada componente,

uα α=ΨΨΨΨ d (61)

com xα = ou yα = , reunindo o vector-linha ΨΨΨΨ as funções de aproximação,

1 2 NΨ Ψ Ψ=ΨΨΨΨ ⋯ (62)

e o vector-coluna αd as componentes do deslocamentos nodais:

1

2

N

d

d

d

α

αα

α

=

d⋮

x

y 1

2

3

xu

yu

2xd

3xd

1xd

1yd

2yd

3yd

x

y 1 2

3

xu

yu

2xd

3xd

1xd

1yd

2yd

3yd

4

4yd

4xd


37

Alternativamente, pode ser conveniente escrever matricialmente a definição (61)

simultaneamente para ambas as componentes,

x x

y y

u

u

=

ΨΨΨΨΨΨΨΨ

d0d0

(63)

em que 0 é o vector-linha com N coeficientes nulos, ou recorrer à forma mais compacta (35) da

aproximação (59), em que:

x

y

u

u

=

u

=

ΨΨΨΨΨΨΨΨΨΨΨΨ0

0

x

y

=

dd

d

Figura 30: Deslocamentos e forças nodais na placa trapezoidal

Identificam-se na Figura 30 os deslocamentos nodais da placa trapezoidal (assim como as

forças correspondentes, para uso posterior). Como adiante se irá verificar, esta discretização é

inaceitável em termos práticos, pois produz uma aproximação muito fraca da solução para o

carregamento e para as condições de apoio definidas na Figura 28. Para além disso, não é boa

prática usar elementos com graus de aproximação diferentes na análise de uma estrutura.

As expressões que se obtêm para o vector (62) que reúne as funções de aproximação das

componentes do deslocamento são as seguintes para cada elemento, de acordo com as definições

dadas nas Tabelas 1 e 2:

(1) (1) (1) (1)1 2 3

1 12 21 x x y y

Ψ Ψ Ψ=

− −

ΨΨΨΨ (64)

1

3 3xd

2xd

1xd

1yd

2yd

3yd

4

4yd

4xd

1 1

3xd

3yd

2xd

1xd

1yd

2yd

1

3

2 2

2

1

3 4

1

3

2 2

1xF

1xF 2xF

2xF

3xF

3xF 4xF

4yF

3yF 3yF

2yF

2yF 1yF

1yF


38

(2) (2) (2) (2) (2)1 2 3 4

1 1 1 12 2 2 2(2 )(2 ) (1 )(2 ) (1 ) (2 )x y x y x y x y

Ψ Ψ Ψ Ψ=

− − − − − − − −

ΨΨΨΨ (65)

10.2 Continuidade do Campo de Deslocamento

A condição de continuidade dos deslocamentos no domínio dos elementos é implicitamente

assegurada pelas funções de aproximação. Para além disso, e de acordo com a análise feita na

Secção 9, essas funções variam linearmente nas fronteiras dos elementos do exemplo em análise,

pois ambos os elementos considerados têm 2 nós por lado, como se ilustra na Figura 31 para o

deslocamento (1) (2)3 4x xd d d= = .

a) Deformadas dos elementos b) Deformada da estrutura

Figura 31: Deformadas devidas ao deslocamento (1) (2)3 4x xd d d= =

Por exemplo, no lado definido pelos nós 2 e 3 do elemento 1, o segmento definido por 1x = e

0 2y≤ ≤ (ver Figura 28) tem-se, usando a expressão (60) com (1)3xd d= , e a definição (64):

(1) (1) 1

3 2(1)

( 1, ) ( 1, ) ( )( 1)

( 1, ) 0x

y

u x y x y d y dx

u x y

ΨΓ

= = = == = =

Se se repetir o processo para o lado definido pelos nós 1 e 4 do elemento 2, usando agora a

expressão (60) com (2)4xd d= e a definição (65) para as funções do elemento, obtém-se o mesmo

resultado final:

(2) (1) 1

4 2(2)

( 1, ) ( 1, ) ( )( 1)

( 1, ) 0x

y

u x y x y d y dx

u x y

ΨΓ

= = = == = =

Consequentemente, os dois elementos têm deslocamentos compatíveis ao longo do lado que

partilham, mostrando-se na Figura 31b) a deformada cinematicamente admissível que se obtém

para a placa trapezoidal quando se reúnem os elementos.

É importante notar que se simplifica nessa figura a representação das condições de apoio

recorrendo a apoios rígidos pontuais. Como, para a aproximação feita, os deslocamentos variam

1

3 4

1

3

2 2

d d d


39

linearmente ao longo dos lados, as condições de encastramento deslizante são exactamente

reproduzidas colocando apoios móveis nos nós afectados por essa condição de apoio. O nó fixo

modela a combinação do efeito das duas condições de encastramento deslizante.

10.3 Aproximação do Campo de Deformação

Para garantir que a aproximação é cinematicamente admissível no sentido forte, a deformação

causada pela aproximação (59) é determinada impondo a condição de compatibilidade no

domínio (14), ficando:

1

1

1 1

Ni

xx ixi

Ni

yy iyi

N Ni i

xy ix iyi i

dx

dy

d dy x

Ψε

Ψε

Ψ Ψγ

=

=

= =

∂=∂

∂=∂

∂ ∂= +∂ ∂

∑

∑

∑ ∑

(66)

De acordo com esta aproximação, quando se provoca um deslocamento no nó j do elemento

segundo a direcção x, jxd , e se impedem todos os outros deslocamentos nodais, 0jyd = e

0ix iyd d= = para i j≠ , o elemento não tem deformação axial segundo y,

; 0;j jxx jx yy xy jxd d

x y

Ψ Ψε ε γ

∂ ∂= = =

∂ ∂

obtendo-se um resultado análogo, mas segundo y, quando se impõe um deslocamento nodal

nessa direcção:

0; ;j jxx yy jy xy jyd d

y x

Ψ Ψε ε γ

∂ ∂= = =

∂ ∂

Na notação matricial (36), é a seguinte a forma explícita da matriz (37) que define as

deformações compatíveis com a aproximação (59) dos deslocamentos:

1 2

1 2

1 2 1 2

0 0 0

0 0 0x x x N

y y y N

y y y N x x x N

Ψ Ψ ΨΨ Ψ Ψ

Ψ Ψ Ψ Ψ Ψ Ψ

∂ ∂ ∂ = ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂

B

⋯ ⋯

⋯ ⋯

⋯ ⋯

(67)

A partição feita nesta definição distingue as deformações devidas aos deslocamentos nodais

segundo x e segundo y, xd e yd , respectivamente, por ser por vezes conveniente escrever a

equação (36) na forma,

x

x y x x y yy

= = +

de B B B d B d

d (68)


40

em que, de acordo com a definições (2) e (67):

1 2

1 2

0 0 0x x x N

x

y y y N

Ψ Ψ Ψ

Ψ Ψ Ψ

∂ ∂ ∂ = ∂ ∂ ∂

B

⋯

⋯

⋯

1 2

1 2

0 0 0

y y y y N

x x x N

Ψ Ψ ΨΨ Ψ Ψ

= ∂ ∂ ∂ ∂ ∂ ∂

B

⋯

⋯

⋯

As definições que se obtêm para os elementos da placa trapezoidal são as seguintes,

(1)

2 2 0 0 0 01

0 0 0 0 1 12

0 1 1 2 2 0

− = − − −

B (69)

(2)

2 2 0 0 0 01

0 0 0 0 2 1 1 22

2 1 1 2 2 2

y y y y

x x x x

x x x x y y y y

− + − − = − + − − + − − + − − + − − + − −

B (70)

mostrando que a aproximação do campo de deformação é constante em elementos triangulares

de 3 nós e linear em elementos rectangulares de 4 nós.

10.4 Continuidade do Campo de Deformação

Estes resultados permitem analisar o campo de deformações ao longo do lado partilhado pelos

dois elementos de discretização da placa trapezoidal quando se impõe a deformada representada

na Figura 31a).

O campo de deformações no elemento 1 é determinado multiplicando a terceira coluna da

matriz de compatibilidade (69) por (1)3xd d= ,

(1) (1) (1) 120; 0; ( )xx yy xy dε ε γ= = =

e o campo de deformações no elemento 2 obtém-se multiplicando a quarta coluna da matriz de

compatibilidade (70) por (2)4xd d= , para 1x = :

(2) (2) (2)1 12 2( ) ; 0; ( )xx yy xyy d dε ε γ= − = =

Conclui-se, portanto, que a deformação axial yyε e a distorção xyγ são contínuas entre

elementos (o que geralmente não sucede nas soluções obtidas com elementos finitos) mas que a

deformação axial xxε não o é. No entanto, a deformada definida na Figura 31 é compatível no

domínio e na fronteira da placa, pois a definição de admissibilidade cinemática não exige a

continuidade do campo de deformações.


41

Tal como para as peças lineares, e pelas mesmas razões, as deformações em estados planos de

tensão ou de deformação podem ser descontínuas, o que sucederá se houver descontinuidades

nas propriedades do material estrutural ou nas forças aplicadas no domínio da estrutura.

10.5 Aproximação do Campo de Tensão

A aproximação do campo de tensão é determinada impondo localmente (ou de maneira forte) a

relação de elasticidade (9) sobre a aproximação (66) do campo de deformação. No entanto, é

mais prático obter esse resultado aplicando a definição matricial (38), ou a definição (68) se se

desejar isolar a contribuição dos deslocamentos nodais:

x

x y x x y yy

= = +

ds DB DB DB d DB d

d (71)

Os resultados que se obtêm para o exemplo da placa trapezoidal (estado plano de tensão) são

os seguintes, de acordo com a organização (1) do vector das componentes de tensão, quando se

admite que o módulo de elasticidade é constante e o coeficiente de Poisson 0,3ν = :

(1)

40 40 0 0 6 65

12 12 0 0 20 20182

0 7 7 14 14 0

E− − = − − − −

DB (72)

(2)

20(2 ) 20(2 ) 20 205

6(2 ) 6(2 ) 6 6182

7(2 ) 7 (1 ) 7 (1 ) 7(2 )x

y y y yE

y y y y

x x x x

− − − − = − − − − − − − − − −

DB (73)

(2)

6(2 ) 6(1 ) 6(1 ) 6(2 )5

20(2 ) 20(1 ) 20(1 ) 20(2 )182

7(2 ) 7(2 ) 7 7y

x x x xE

x x x x

y y y y

− − − − − − = − − − − − − − − − −

DB (74)

Estes resultados mostram que a aproximação do campo de tensão é constante em elementos

triangulares de 3 nós e linear no elemento rectangular de 4 nós. Mostram, também, que as

componentes do campo de tensão devidas a um deslocamento nodal são em geral não nulas,

devido ao efeito de Poisson.

Os campos de tensão nos elementos 1 e 2 da placa trapezoidal devidos ao deslocamento nodal

d, ver Figura 31, são definidos multiplicando a terceira coluna da matriz (72) por (1)3xd d= ,

(1) (1) (1) 351820; 0; ( )E

xx yy xy dσ σ σ= = = (75)

e a quarta coluna da matriz (73) por (2)4xd d= :

(2) (2) (2)100 30 35182 182 182( ) ( ) ; ( ) ( ) ; ( ) (2 )E E E

xx yy xyy d y d x dσ σ σ= − = − = − (76)


42

Estes resultados mostram que a aproximação feita para o campo de deslocamentos está

associada a campos de tensão que são descontínuos ao longo do lado comum dos elementos em

que a placa foi discretizada.

10.6 Equilíbrio do Campo de Tensão no Domínio do Elemento

A aplicação das condições de equilíbrio no domínio, definidas pela equação (13) na forma,

ˆ

ˆx x xx y xy

y y yy x xy

fem

f

σ σΩ

σ σ

= −∂ − ∂

= −∂ − ∂

permite determinar as forças de massa que equilibram o campo de tensões. De acordo com as

definições (8) e (38), a expressão deste resultado é a seguinte,

Tˆ ( ) emΩ−f = A DB d (77)

ou, se se usar a expressão (71):

T Tˆ ( ) ( )x x y y emΩ− −f = A DB d A DB d

É importante notar que estas forças não representam as forças de massa do problema que se

estiver a analisar, f , por exemplo as definidas na Figura 28 para a placa trapezoidal, mas as

forças que equilibram a aproximação feita para o campo de tensões.

O resultado (72) mostra que um elemento de 3 nós só pode equilibrar forças de massa nulas.

Conclui-se, portanto, que no exemplo em análise a condição de equilíbrio no elemento 1 é

sempre violada, independentemente dos valores que se venham a obter para os deslocamentos

nodais do elemento.

Segundo os resultados (73) e (74), a aproximação do campo de tensões num elemento de 4

nós só pode equilibrar forças de massa nulas ou constantes. No entanto, e como adiante se

ilustra, mesmo quando as forças de massa são constantes, como no caso do peso próprio,

verifica-se que a solução obtida aplicando o Método dos Elementos Finitos não satisfaz, em

geral, a condição de equilíbrio no domínio dos elementos.

Exercício 10: Verifique serem as seguintes as forças de massa equilibradas pelo campo de

tensão no elemento 2 da placa trapezoidal devido ao deslocamento nodal d:

(2) (2) 514

ˆ ˆ0; ( )Ex yf f d= = (78)

10.7 Equilíbrio do Campo de Tensão na Fronteira do Elemento

A aplicação das condições de equilíbrio (23),

ˆ

ˆx x xx y xy

y y yy x xy

t n nem

t n n

σ σΓ

σ σ = + = +

(79)


43

permite determinar as forças que equilibram a aproximação do campo de tensões na fronteira Γ

de um dado elemento. Obtêm-se as seguintes expressões gerais para este resultado usando as

definições (11) e (38) ou (71):

Tˆ ( ) emΓt = N DB d (80)

T Tˆ ( ) ( )x x y y emΓ+t = N DB d N DB d (81)

Como os lados dos elementos de 3 e 4 nós são sempre definidos por segmentos de recta, as

componentes da normal exterior unitária, xn e yn , são constantes, concluindo-se que essas forças

são constantes para o elemento de 3 nós e têm uma variação linear para o elemento de 4 nós.

Retomando o efeito do deslocamento d na aproximação feita sobre a placa trapezoidal, de

acordo com os resultados (75) e (76) e com a definição (79), as forças que equilibram esses

estados de tensão ao longo do lado que os elementos partilham são as seguintes para cada

elemento,

(1) (1) (1)

(1) (1) (1) 526

ˆ ( 1) (0) 0( 1)

ˆ (0) ( 1) ( )x xx xy

Ey yy xy

tx

t d

σ σΓ

σ σ = + + == = + + =

(2) (2) (2) 100

182(2) (2) (2) 35

182

ˆ ( 1) (0) ( ) ( )( 1)

ˆ (0) ( 1) ( )

Ex xx xy

Ey yy xy

t y dx

t d

σ σΓ

σ σ = − + == = + − = −

mostrando que a deformada representada na Figura 31 está associada a campos de tensão que

violam a condição de equilíbrio de forças entre os elementos segundo a direcção x, qualquer que

seja o valor (infinitesimal) do deslocamento nodal, 0d ≠ ,

(1) (2)

(1) (2)

ˆ ˆ 0( 1)

ˆ ˆ 0x x

y y

t tx

t tΓ

+ ≠= + =

as quais deveriam ser iguais e de sinal contrário, como em qualquer diagrama de corpo livre, tal

como impõe a condição (27) e se verifica (neste caso particular) para o equilíbrio de forças

segundo a direcção y.

Verifica-se que, em geral, a solução obtida aplicando o Método dos Elementos Finitos não

satisfaz as condições de equilíbrio na fronteira entre elementos. O mesmo se passa na fronteira

estática da malha de elementos finitos, como já se ilustrou na análise da consola triangular (ver

Figuras 12 e 16).

As condições de equilíbrio não são, por isso, impostas de maneira forte (isto é, em todos os

pontos) mas de maneira fraca recorrendo ao conceito de força nodal equivalente, através da

condição de equilíbrio nodal (39). Ilustra-se na Secção 11 a determinação dos termos dessa

equação.


44

Exercício 11: Verifique serem as forças representadas na Figura 32 que equilibram os campos

de tensão nos elementos 1 e 2 da placa trapezoidal devido ao deslocamento nodal 1d = :

a) Forças segundo x b) Forças segundo y

Figura 32: Forças de fronteira em equilíbrio com os campos de tensão

Exercício 12: Verifique que é globalmente equilibrado o sistema de forças aplicadas nos

elementos 1 e 2 da placa trapezoidal devidas ao deslocamento 1d = . Nota: Para o elemento 2

é necessário considerar a contribuição das forças de massa (78) que equilibram o campo de

tensão no domínio do elemento.

11. Forças Nodais Equivalentes

Prosseguindo para a segunda fase da aplicação do Método dos Elementos Finitos, resumida na

Secção 7, ilustra-se e interpreta-se o cálculo dos termos da matriz de rigidez (40), do vector das

forças nodais equivalentes às forças de massa (42) e do vector das forças nodais equivalentes às

forças de fronteira (43) escritas a nível de elemento.

Usam-se como aplicação os elementos de 3 e 4 nós definidos nas Tabelas 1 e 2, identificando-

se na Figura 33 as forças correspondentes aos deslocamentos representados na Figura 29.

a) Elemento triangular de 3 nós b) Elemento rectangular de 4 nós

Figura 33: Convenção para a medição das componentes do campo de forças

1

3 4

1

3

2 2

(2) 35182

ˆ (2 )xt x E= −

(2) 35182

ˆ ( 2)xt x E= −

(2) 100182xt yE= −

xt

(2) 100182xt yE=

(1) 526xt E=

(1)ˆ 0xt = (1) 526xt E= −

1

3 4

1

3

2 2

yt

(2)ˆ 0yt =

(2)ˆ 0yt =

(2) 200182yt E= − (2) 35

182yt E= −

(1)ˆ 0yt =

(1) 513xt E= −

(1) 526yt E=

x

y 1

2

3

xf

yf

1xF

2xF

3xF

3yF

2yF

1yF

x

y 1 2

3 4

1xF

2xF

3xF

4xF

4yF

3yF

2yF

1yF

xf

yf


45

11.1 Vector das Forças Nodais Equivalentes às Forças de Massa

Como adiante se mostra, num elemento triangular de 3 nós as forças nodais equivalentes a forças

de massa constantes,

x x y yf f e f f emΩ= = (82)

são obtidas calculando as resultantes dessas forças e distribuindo-as uniformemente pelos nós do

elemento, ou seja,

1 13 3ix x iy yF A f e F A f= = (83)

em que A representa a área do elemento. O resultado é análogo para o elemento rectangular de 4

nós, como se ilustra nas Figuras 34 e 35:

1 14 4ix x iy yF A f e F A f= = (84)

a) Forças de massa constantes, fα ,

elemento com área A

b) Forças nodais equivalentes, 13iF F A fα α α= =

Figura 34: Forças nodais equivalentes às forças de massa no elemento de 3 nós

a) Forças de massa constantes, fα ,

elemento com área A

b) Forças nodais equivalentes, 14iF F A fα α α= =

Figura 35: Forças nodais equivalentes a forças de massa no elemento de 4 nós

x

y 1

2

3

xf

yf

x

y 1

2

3 xF

xF

yF

yF

xF

yF

x

y 1 2

3 4

xf

yf

x

y 1 2

3 4

xF

xF

xF

xF

yF

yF

yF

yF


46

Aplicando directamente os resultados (83) e (84) e adoptando a identificação das forças

nodais definida na Figura 30, obtêm-se as seguintes expressões para as forças nodais

equivalentes (em 1kNm− ), representadas na Figura 36:

1

2

3(1)

1

2

3

0

0

0

20 / 3

20 / 3

20 / 3

x

x

x

fy

y

y

F

F

F

F

F

F

= = − −

−

F

1

2

3

4(2)

1

2

3

4

0

0

0

0

10

10

10

10

x

x

x

x

fy

y

y

y

F

F

F

F

F

F

F

F

= = − −

− −

F (85)

Figura 36: Forças nodais equivalentes às forças de massa da placa trapezoidal

Os resultados (83) e (84) decorrem da aplicação da definição geral (42) do vector das forças

nodais equivalentes às forças de massa. Se se introduzir a partição correspondente às

componentes em x e em y, de acordo com a relação (63),

T

T

x x

fy y

fd

fΩΩ

= =

∫

ΨΨΨΨΨΨΨΨ

F 0F

F 0

obtêm-se as seguintes expressões,

T f dα αΩΩ= ∫ ΨΨΨΨF (86)

com xα = ou yα = , ou para cada componente, de acordo com as definições (59) e (62):

i iF f dα αΩΨ Ω= ∫ (87)

Exercício 13: Verifique os resultados (83) e (84) para os elementos de 3 e 4 nós usados na

discretização da placa trapezoidal representada na Figura 28 substituindo as definições (64) e

(65) para as funções de aproximação na expressão geral da definição (87), ficando para os

elementos 1 e 2, respectivamente:

1 1

0

20x

y

f

f

== −

0

20x

y

f

f

== −

2m

1

3 4

1

3

2 2

203

203

203

10 10

10 10


47

121 2 2

0 0 0 0

x y

i i iF f dy dx f dx dyα α αΨ Ψ= =∫ ∫ ∫ ∫

2 2 2 2

1 0 0 1i i iF f dy dx f dx dyα α αΨ Ψ= =∫ ∫ ∫ ∫

Os resultados (83) e (84) são apenas válidos para forças de massa constantes e para os

elementos caracterizados nas Tabelas 1 e 2. Os resultados serão diferentes se as forças de massa

não forem constantes e/ou as funções de aproximação não forem lineares em elementos

triangulares ou bilineares em elementos rectangulares.

Todavia, o que a definição (87) assegura, para qualquer tipo de elemento e para qualquer

expressão das forças de massa, é que as forças nodais equivalentes realizam sobre os

deslocamentos nodais correspondentes o mesmo trabalho que as forças de massa realizam sobre

a aproximação (59) do campo de deslocamentos no domínio do elemento, ou seja,

( )i i i id F d f d u f dα α α α α αΩ ΩΨ Ω Ω= =∫ ∫ (88)

para o campo devido ao deslocamento nodal idα (nó i, direcção xα = ou yα = ) quando todos

os restantes são nulos.

Exercício 14: Recorrendo à generalização (18), defina o vector das forças nodais equivalentes

a uma variação de temperatura uniforme nos elementos da placa trapezoidal.

11.2 Vector das Forças Nodais Equivalentes às Forças de Fronteira

Se na definição geral (43) do vector das forças nodais equivalentes às forças de fronteira se

introduzir a partição das componentes em x e em y obtêm-se expressões análogas às

determinadas para as forças de massa, com a diferença que são agora definidas sobre a fronteira,

Γ , onde as forças estão aplicadas:

T

T

x x

ty y

td

tΓΓ

= =

∫

ΨΨΨΨΨΨΨΨ

F 0F

F 0

Tt dα αΓΓ= ∫ ΨΨΨΨF (89)

i iF t dα αΓΨ Γ= ∫ (90)

Tal como para as forças de massa, de acordo com a definição (90) uma força de fronteira

aplicada segundo uma dada direcção, x ou y, é substituída por forças nodais equivalentes com a

mesma direcção, ixF ou iyF . A principal diferença entre as definições (86) e (89) é que as forças

nodais equivalentes às forças de massa são distribuídas por todos os nós do elemento enquanto as

forças nodais equivalentes às forças de fronteira são apenas distribuídas pelos nós contidos no


48

lado do elemento em que essas forças estão aplicadas. Esta característica resulta da condição (57)

imposta sobre as funções de aproximação, a qual tem como consequência que a função associada

ao nó, i, iΨ , tem valor nulo nos lados do elemento que não contêm esse nó.

a) Elemento de 3 nós b) Elemento de 4 nós

Figura 37: Definição das fronteiras dos elementos de 3 e 4 nós

Definem-se na Figura 37 as equações dos lados dos elementos de 3 e 4 nós caracterizados nas

Secções 9.1 e 9.2, respectivamente, para apoiar o cálculo das forças nodais equivalentes a forças

de fronteira. Por exemplo, se no lado definido pelos nós 1 e 4 do elemento de 4 nós estiver

aplicada uma força constante, t pα = , a expressão (90),

( )1

11

y b

i i x xyF p dyα Ψ

+

== ∫

produz as seguintes definições para as forças nodais equivalentes, de acordo com os dados da

Tabela 2, em que 1'x x x= − e 1'y y y= − :

( ) ( )( )1

11

11 1 2' 00

1 '/ 1 '/ 'y b b

x x xyF p dy x a y b p dy pbα Ψ

+

= == = − − = ∫ ∫

( ) ( )( )1

112 2 ' 00

1 '/ '/ ' 0y b b

x x xyF p dy y b x a p dyα Ψ

+

= == = − = ∫ ∫

( ) ( ) ( )1

113 3 ' 00

'/ '/ ' 0y b b

x x xyF p dy x a y b p dyα Ψ

+

= == = = ∫ ∫

( ) ( ) ( )1

11

14 4 2' 00

1 '/ '/ 'y b b

x x xyF p dy x a y b p dy pbα Ψ

+

= == = − = ∫ ∫

Ou seja, a resultante da força aplicada no lado com comprimento b é repartida igualmente

pelos dois nós de extremidade.

Exercício 15: Mostre que se obtém o mesmo resultado para os restantes lados do elemento de

4 nós representado na Figura 37, assim como para qualquer lado do elemento de 3 nós com a

geometria geral definida na mesma figura.

x

y

1 1( , )x y

2 2( , )x y

3 3( , )x y

2 2 2 2 0y xΨ γ β α= + + =

1 1 1 1 0y xΨ γ β α= + + =

3 3 3 3 0y xΨ γ β α= + + =

x

y

1 1( , )x y 2 2( , )x y

3 3( , )x y

4 4( , )x y

1x x=

1y y=

1y y b= +

1x x a= +


49

Exercício 16: Determine os vectores das forças nodais equivalentes às forças de fronteira para

as duas discretizações da consola quadrada representada na Figura 38.

Figura 38: Consola quadrada sujeita a uma carga uniforme (estado plano de tensão)

Figura 39: Forças nodais equivalentes a forças de fronteira lineares em lados com 2 nós

Como as funções de aproximação do campo de deslocamentos, iΨ , variam linearmente ao

longo dos lados dos elementos de 3 e 4 nós, a definição (90) produz o mesmo resultado para a

mesma variação das forças aplicadas, como se ilustra na Figura 39 para forças de fronteira

lineares aplicadas em lados definidos por dois nós. No entanto, essa definição produz diferentes

resultados para elementos baseados em diferentes funções de aproximação, mantendo-se a

condição de definir forças nodais que realizam sobre os deslocamentos nodais o mesmo trabalho

que as forças de fronteira realizam sobre a aproximação dos deslocamentos nessa fronteira:

( )i i i id F d t d u t dα α α α α αΓ ΓΨ Γ Γ= =∫ ∫

Figura 40: Forças nodais equivalentes às forças aplicadas na fronteira da placa trapezoidal

1m

x

y

1m

2. ( )

0,3

0x y

E const kNm

f f

ν

−=== =

2( )p kNm−

Malha 1

1

3 4

2

Malha 2

1 1

2

2 3 3

kL

i j

xp 2

3ix xF R= 13jx xF R=

12 kR p Lα α=

i j

13jy yF R= 2

3iy yF R= yp

1

3 4

2

10

20 kPa

1

3

2

1

1

3 4

2 1

3

2

4yF

3yF

1 23 3 3( 5) ( 10) 5yF = + + − = −

2 14 3 3( 5) ( 10) 0yF = + + − =


50

São as seguintes as expressões dos vectores das forças nodais equivalentes ( 1kNm− ) às forças

de fronteira para o exemplo da placa trapezoidal, as quais estão representadas na Figura 40:

1

2

3(1)

1

2

3

0

0

0

0

0

0

x

x

x

ty

y

y

F

F

F

F

F

F

= =

F

1

2

3

4(2)

1

2

3

4

0

0

0

0

0

0

5

0

x

x

x

x

ty

y

y

y

F

F

F

F

F

F

F

F

= =

−

F (91)

11.3 Matriz de Rigidez

Se na definição geral (40) da matriz de rigidez se introduzir a partição correspondente às

partições (68) e (71) das deformações e das tensões devidas às componentes em x e em y dos

deslocamentos nodais obtém-se a seguinte expressão,

T

T

xx xy xx y

yx yy y

dΩ

Ω

= =

∫Κ ΚΚ ΚΚ ΚΚ Κ ΒΒΒΒ

Κ Β ΒΚ Β ΒΚ Β ΒΚ Β ΒΚ ΚΚ ΚΚ ΚΚ Κ ΒΒΒΒD

verificando-se as seguintes condições de simetria:

T T dαα αα α αΩΩ= = ∫ Β ΒΒ ΒΒ ΒΒ ΒK K D

T Txy yx x ydΩ

Ω= = ∫ Β ΒΒ ΒΒ ΒΒ ΒK K D

O termo geral da matriz pode ser determinado pela expressão seguinte, devendo ter-se em

conta quais são os deslocamentos envolvidos, segundo x ou segundo y, na identificação das

colunas iΒΒΒΒ e jΒΒΒΒ da matriz que define as deformações compatíveis com esses deslocamentos:

Tij ji i jK K d

ΩΩ= = ∫ Β ΒΒ ΒΒ ΒΒ ΒD (92)

Por exemplo, o termo 53K da matriz de rigidez do elemento 1 da discretização da placa

trapezoidal é definido pelas quinta e terceira colunas das matrizes (69) e (72),

1 25

53 35 1820 0

6

0 0 1 20

14

xEK K dy dx

− = = − −

∫ ∫

e o termo 46K da matriz de rigidez do elemento 2 é definido pelas quarta e segunda colunas das

matrizes (70) e (74), se se usar a partição anteriormente descrita:


51

2 251 1

46 64 2 2 1821 0

6(1 )

0 1 20(1 )

7(2 )

E

x

K K y x x dy dx

y

− = = − − − −

∫ ∫

Exercício 17: Usando um programa de computação simbólica, aplique a definição (40) e os

dados resumidos nas Tabelas 1 e 2 para obter a expressão geral da matriz de rigidez do

elemento de 3 nós,

1

4 1ij xx i j i j

GK

A

κγ γ β βκ

+ = + − (93)

1

4 1ij yy i j i j

GK

A

κβ β γ γκ

+ = + − (94)

3

4 1ij xy ji yx i j i j

GK K

A

κγ β β γκ

− = = − − (95)

e do elemento rectangular de 4 nós, em que /b aβ = e 2( 1)/( 1)α β κ κ= + − :

2 2 1 2 1 2

1 2 2 2 2 1

1 2 2 2 1 26

2 1 1 2 2 2

xx yy

G

α α α αα α α αα α α αβα α α α

+ − − − − − − + − − − − = = − − − − + − − − − − − +

K K (96)

T

1 2 1 2

2 1 2 1

1 2 1 22( 1)

2 1 2 1

xy yx

G

κ κκ κ

κ κκκ κ

+ − − − − − − + = = − − + −− − + − −

K K (97)

As expressões gerais (93) a (97) para os coeficientes da matriz de rigidez, assim como as

anteriormente dadas para o vector das forças nodais equivalentes, são úteis na aplicação manual

do Método dos Elementos Finitos. Os programas de cálculo baseiam-se em processos mais

gerais e eficazes, sumariamente descritos na Secção 17.

Exercício 18: Particularize os resultados do exercício anterior para obter as expressões das

matrizes de rigidez elementares para o exemplo da consola trapezoidal, definido na Figura 28:

(1)

400 400 0 0 60 60

400 435 35 70 130 60

0 35 35 70 70 0

0 70 70 140 140 0364

60 130 70 140 240 100

60 60 0 0 100 100

E

− − − − − − −

= − − − − −

− −

K (98)


52

(2)

290 255 145 110 65 5 65 5

255 290 110 145 5 65 5 65

145 110 290 255 65 5 65 5

110 145 255 290 5 65 5 65

65 5 65 5 160 60 80 20364

5 65 5 65 60 160 20 80

65 5 65 5 80 20 160 60

5 65 5 65 20 80 60 160

E

− − − − − − − − − − − − − − − − = − − − − −

− − − − − − − − − −

− − − − −

K

(99)

Tal como o termo independente da equação resolvente (39) combina as forças equivalentes às

forças de massa e de fronteira, o termo geral da matriz de rigidez define as forças equivalentes às

forças de massa (77) e de fronteira (80) que equilibram os campos de tensão provocados pela

aproximação do campo de deslocamentos: ijK representa a força nodal equivalente iF devida à

acção do deslocamento nodal 1jd = quando são nulos todos os restantes deslocamentos nodais.

Para demonstrar esta interpretação, introduz-se na definição (40) da matriz de rigidez a

expressão da matriz (37) que define as deformações compatíveis,

T( ) dΩ

Ω= ∫ ΨΨΨΨK A D B

e integra-se por partes para obter as forças nodais (42) e (43) equivalentes às forças de massa

(77) e às forças de fronteira (80) provocadas pelos deslocamentos nodais:

T T T T( ) ( )d dΩ Γ

Ω Γ= − +∫ ∫Ψ ΨΨ ΨΨ ΨΨ ΨK A D B N D B (100)

Esta interpretação é ilustrada para a deformada definida na Figura 31. Estão representadas na

Figura 41a) as forças equivalentes às forças de massa que equilibram o campo de tensões, nulas

no elemento 1 e dadas pela equação (78) para o elemento 2. As forças equivalentes às forças que

equilibram esse campo de tensões na fronteira dos elementos (ver Figura 32) estão representadas

na Figura 41b). As resultantes dessas forças recuperam os coeficientes da terceira e quarta

colunas das matrizes de rigidez (98) e (99) do elemento 1, (1)3 1xd = , e do elemento 2, (2)

4 1xd = :

(1) (1)1 13(1) (1)

2 23(1) (1)

3 33(1) (1)

1 43(1) (1)

2 53(1) (1)

3 63

0

35

35

70364

70

0

x

x

x

y

y

y

F K

F K

F K E

F K

F K

F K

− = = −

(2) (2)1 14(2) (2)

2 24(2) (2)

3 34(2) (2)

4 44(2) (2)

1 54(2) (2)

2 64(2) (2)

3 74(2) (2)

4 84

110

145

255

290

5364

65

5

65

x

x

x

x

y

y

y

y

F K

F K

F K

F K E

F K

F K

F K

F K

− − = = − −

(101)


53

Forças segundo x Forças segundo y

a) Forças nodais equivalentes às forças de massa, Eq. (77)


b) Forças nodais equivalentes às forças de superfície, Eq. (80)


c) Coeficientes da matriz de rigidez, Eq. (101)

Figura 41: Forças nodais devidas ao deslocamento (1) (2)3 4x xd d d= = na placa trapezoidal

Exercício 19: Verifique e represente o significado físico de uma coluna das matrizes de

rigidez dos elementos de 3 nós da consola quadrada usando os resultados (93) a (95):

(1)

35 0 35 0 35 35

0 100 100 30 0 30

35 100 135 30 35 65

0 30 30 100 0 100182

35 0 35 0 35 35

35 30 65 100 35 135

E

− − − − − − −

= − − − −

− − −

K (102)

1

3 4

1

3

2 2

1

3 4

1

3

2 2

528E 5

28E

528E 5

28E

1

3 4

1

3

2 2

552E

552E

552E− 5

52E−

10156

E

10156

E− 5156

E−

5156

E

200273

E 200273

E−

100273

E− 100273

E

1

3 4

1

3

2 2

526E− 5

26E

526E−

526E

1591

E− 526E− 15

91E−

526E−

1

3 4

1

3

2 2

(1)13K (1)

23K

(1)33K

(2)14K (2)

24K

(2)34K (2)

44K

1

3 4

1

3

2 2

(1)43K (1)

53K

(1)63K

(2)54K (2)

64K

(2)74K (2)

84K


54

(2)

100 100 0 0 30 30

100 135 35 35 65 30

0 35 35 35 35 0

0 35 35 35 35 0182

30 65 35 35 135 100

30 30 0 0 100 100

E

− − − − − − −

= − − − − −

− −

K (103)

Exercício 20: Verifique e represente o significado físico de uma das colunas das matrizes de

rigidez do elemento de 4 nós da consola quadrada aplicando os resultados (96) e (97), e

calcule um coeficiente aplicando a definição (92):

(1)

180 110 90 20 65 5 65 5

110 180 20 90 5 65 5 65

90 20 180 110 65 5 65 5

20 90 110 180 5 65 5 65

65 5 65 5 180 20 90 110364

5 65 5 65 20 180 110 90

65 5 65 5 90 110 180 20

5 65 5 65 110 90 20 180

E

− − − − − − − − − − − − − − − − = − − − −

− − − − − − − −

− − − −

K (104)

Exercício 21: Utilize os resultados (77), (80) e (100) para mostrar que a matriz de rigidez, K ,

define forças nodais equivalentes, K d , que realizam sobre os deslocamentos nodais um

trabalho, Td K d , equivalente ao realizado pelas forças de massa e de fronteira sobre a

aproximação do campo de deslocamentos: T T Tˆ ˆ= +d K d d f d t .

12. Formulação da Equação Resolvente

Analisa-se a seguir a terceira fase do processo de aplicação do Método dos Elementos Finitos

resumido na Secção 7. Mostra-se que a equação resolvente (46) da malha de elementos finitos é

obtida combinando os vectores das forças nodais e as matrizes de rigidez elementares usando um

processo de reunião análogo ao adoptado na análise de estruturas reticuladas.

12.1 Deslocamentos e Forças Nodais da Malha de Elementos Finitos

A identificação dos deslocamentos nodais independentes de uma malha de elementos finitos, q ,

consiste apenas em seleccionar os movimentos (de translação) livres em cada nó. Esses

deslocamentos, reunidos no vector q e cuja dimensão define o grau de indeterminação

cinemática da malha, β , são medidos no referencial global, pois é nesse referencial que são

definidos os deslocamentos nodais dos elementos. As forças correspondentes definem o vector

das forças nodais, Q .


55

Os deslocamentos e as forças nodais para o exemplo da placa trapezoidal estão definidos na

Figura 42, sendo portanto 5β = o grau de indeterminação cinemática para esta discretização

muito grosseira:

T1 2 3 4 5q q q q q=q

T1 2 3 4 5Q Q Q Q Q=Q

Figura 42: Deslocamentos nodais independentes e forças nodais correspondentes

No vector das forças nodais distinguem-se as parcelas associadas a cada tipo de carregamento,

designadamente as forças nodais equivalentes às forças de massa, fQ , às forças de fronteira, tQ ,

às tensões ou deformações residuais que possam existir, rQ , para além das forças concentradas

que possam estar aplicadas nos nós da malha, nQ :

f t n= + +Q Q Q Q (105)

Admite-se aqui que as forças nodais equivalentes a tensões ou deformações residuais, como

por exemplo as deformações térmicas, são incluídas no vector das forças nodais equivalentes às

forças de massa. Essas forças são nulas para o exemplo da placa trapezoidal, assim como as

forças nodais aplicadas, para o carregamento definido na Figura 28.

12.2 Incidência Nodal dos Deslocamentos Indeterminados

Na análise de estruturas reticuladas, para garantir a continuidade dos deslocamentos entre

elementos unidimensionais é suficiente impor que os elementos que partilham um mesmo nó

tenham aí deslocamentos idênticos aos deslocamentos da estrutura (continuidade nodal).

O mesmo acontece na análise de estados planos. Como se mostrou na Secção 10.2, as funções

de aproximação dos elementos bidimensionais são construídas de modo a assegurar que:

• A condição de continuidade nodal é suficiente para assegurar a continuidade dos

deslocamentos entre dois elementos, desde que esses elementos tenham o mesmo número

de nós no lado que partilham.

1q

2q

3q

4q

5q

5Q

4Q

3Q

2Q

1Q


56

A segunda condição imposta no processo de reunião de elementos de estruturas reticuladas

consiste em igualar as forças aplicadas num nó da estrutura à resultante das forças nodais dos

elementos que partilham esse nó (equilíbrio nodal). É também esse o critério adoptado em

problemas bidimensionais, mas com uma consequência diferente e importante, já ilustrada na

Secção 10.7 e que adiante se confirma:

• O equilíbrio das forças nodais equivalentes não é suficiente para garantir o equilíbrio entre

os elementos que partilham um mesmo nó.

As operações de reunião de elementos são formalmente descritas pelas condições (44) e (45),

servindo a primeira para identificar os deslocamentos nodais dos elementos considerados

isoladamente com os deslocamentos nodais independentes da malha e a segunda para definir as

forças nodais da malha em função das forças nodais dos elementos. No entanto, é mais prático,

mesmo computacionalmente, realizar estas operações recorrendo a uma tabela que resume a

informação presente nas matrizes de incidência nodal, ℑℑℑℑ .

Essa tabela de incidência é a seguinte para o exemplo da placa trapezoidal, de acordo com a

numeração definida nas Figuras 30 e 42 (admitindo que se numeram sequencialmente as

componentes em x e em y dos deslocamentos e forças nodais dos elementos):

Estrutura ( ,i iq Q ) 1 2 3 4 5 Elemento 1 ( (1) (1),i id F ) 1 2 3 6 - Elemento 2 ( (2) (2),i id F ) - 1 4 8 7

Tabela 3: Incidência nodal (deslocamentos indeterminados e forças correspondentes)

Tal como para as estruturas reticuladas, esta tabela serve para estabelecer a condição de

continuidade dos deslocamentos nodais e definir a resultante das forças nodais quando se realiza

a ligação dos elementos, recorrendo à informação resumida nas Figuras 30 e 42:

(1)1 1(1) (2)2 1 2(1) (2)3 4 3(1) (2)6 8 4(2)7 5

d q

d d q

d d q

d d q

d q

= = = = = = = =

(106)

(1)1 1

(1) (2)2 2 1

(1) (2)3 3 4

(1) (2)4 6 8

(2)5 7

Q F

Q F F

Q F F

Q F F

Q F

= = + = + = + =

(107)


57

Exercício 22: Defina matricialmente a relação (106) na forma (44), para cada elemento, e

verifique que a relação (107) tem a expressão geral (45).

Exercício 23: Defina as tabelas de incidência nodal para a discretização da consola quadrada

representada na Figura 38 adoptando a numeração indicada na Figura 43 para os

deslocamentos e forças nodais.

Figura 43: Deslocamentos nodais independentes e forças nodais correspondentes


Sendo nulas as forças concentradas aplicadas nos nós da estrutura, para o carregamento definido

na Figura 28,

T

1 2 3 4 5 0 0 0 0 0n n n n nQ Q Q Q Q =

aplica-se a definição (107) para determinar as forças nodais equivalentes às forças de massa,

T

1 2 3 4 5 0 0 0 50 / 3 10f f f f fQ Q Q Q Q = − −

e às forças de fronteira,

T

1 2 3 4 5 0 0 0 0 5t t t t tQ Q Q Q Q = −

de acordo com os resultados (85) e (91), como se ilustra nas Figuras 36, 39, 44 e 45.

a) Nos elementos b) Na estrutura

Figura 44: Forças nodais equivalentes às forças de massa

Malha 1

1

3 4

2

Malha 2

1 1

2

2 3 3

1 1,q Q

2 2,q Q

3 3,q Q

4 4,q Q

1 1,q Q

2 2,q Q

3 3,q Q

4 4,q Q

1

3

2

(1)4 fF (1)

5 fF

(2)8 10fF = −

1

3 4

2

(2)5 fF (2)

6 fF

(2)7 10fF = − (1) 20

6 3fF = −

1 fQ

2 fQ

3 fQ

4 fQ

5 fQ


58

a) Nos elementos b) Na estrutura

Figura 45: Forças nodais equivalentes às forças de fronteira

É a seguinte a resultante (105) das forças nodais equivalentes (ver Figura 42):

T

1 2 3 4 5 0 0 0 50 / 3 15Q Q Q Q Q = − − (108)

Para obter os coeficientes ijK da matriz de rigidez da estrutura a partir dos coeficientes ( )eijK

das matrizes de rigidez de cada elemento é necessário recorrer às condições de incidência em

termos de deslocamentos (106) e de forças nodais (107). Essa atribuição é feita atendendo a que

ijK representa a força nodal iQ na estrutura devida ao deslocamento unitário jq e que ( )eijK

representa a força nodal ( )eiF no elemento e devida ao deslocamento unitário ( )e

jd , concluindo-se

o seguinte para o exemplo em análise (ver Tabela 3):

(1) (1) (1)11 11 12 12 13 13

(1) (1) (2) (1) (2)21 21 22 22 11 23 23 14

(1) (1) (2) (1) (2)31 31 32 32 41 33 33 44

(1) (1) (2) (1) (2)41 61 42 62 81 43 63 84

(2)51 52 71 530

K K K K K K

K K K K K K K K

K K K K K K K K

K K K K K K K K

K K K K K

= = = = = + = + = = + = + = = + = + = = =

(1)1514 16

(2)(1) (2)25 1724 26 18

(2)(1) (2)35 4734 36 48

(2)(1) (2)45 8744 66 88

(2)(2) (2)55 7774 54 78

0KK K

K KK K K

K KK K K

K KK K K

K KK K

= = == + == +

== +

= =

(109)

A equação resolvente (46) para o exemplo da consola trapezoidal é obtida usando os

resultados (108) e (109), para as definições (98) e (99) das matrizes de rigidez elementares:

1

2

3

4

5

1,099 1,099 0 0.165 0 0

1,099 1,992 0,206 0,179 0,179 0

0 0,206 0,893 0,179 0,014 0

0,165 0,179 0,179 0,715 0,165 16,667

0 0,179 0,014 0,165 0,440 5

q

q

E q

q

q

− − − −

=− − − − − − − −

(110)

Este processo de reunião dos elementos é fácil de automatizar, sendo executado elemento a

elemento, tal como foi descrito no texto sobre a análise de estruturas reticuladas. Apresenta-se na

Figura 46 a interpretação física da terceira coluna da matriz de rigidez da estrutura.

1

3

2

1

3 4

2

(2)8 0tF = (2)

7 5tF = −

1tQ

2tQ

3tQ

4tQ

5tQ


59

a) Deslocamento imposto na malha b) Forças nodais equivalentes

c) Deslocamento imposto nos elementos d) Forças nodais equivalentes

Figura 46: Contribuição do deslocamento nodal 3 1q = para a matriz de rigidez da estrutura

Exercício 24: Verifique as seguintes definições da equação resolvente (46) para o exemplo da

consola quadrada definido nas Figuras 38 e 43 resolvido com um elemento de 4 nós,

1

2

3

4

0,4945 0,0549 0,1786 0,0137 0

0,0549 0,4945 0,0137 0,1786 0

0,1786 0,0137 0,4945 0,3022 0

0,0137 0,1786 0,3022 0,4945 0,5

q

qE

q

q p

− − =

− − − − −

(111)

e com dois elementos de 3 nós:

1

2

3

4

0,7416 0,1923 0,3571 0,1648 0

0,1923 0,7416 0,1923 0 0

0,3571 0,1923 0,7416 0,5495 0

0,1648 0 0,5495 0,7416 0,5

q

qE

q

q p

− − − =

− − − −

(112)

13. Análise da Solução

Depois de resolvida a equação do Método dos Deslocamentos, o resultado obtido é utilizado para

traçar a deformada da estrutura, os campos de tensão e, eventualmente, os campos de

deformação e as reacções de apoio. Estes resultados, determinados na fase de pós-processamento

dos programas comerciais, são apresentados graficamente e, se se desejar, definidos

numericamente em pontos ou em secções da estrutura.

3 1q =

13K

23K

33K

43K

53K

1

3 4

1

3

2 2

(1)3 1d = (2)

4 1d =

1

3 4

1

3

2 2

(1)13K (1)

23K

(1)33K

(1)43K (1)

53K

(1)63K

(2)14K (2)

24K

(2)34K (2)

44K

(2)54K (2)

64K

(2)74K (2)

84K


60

Para ilustrar a fase de análise dos resultados, e retomando o exemplo da placa trapezoidal,

representa-se na Figura 47 a deformada calculada depois de resolver a generalização (110) da

equação do Método dos Deslocamentos (em 810 m− ),

T

1 2 3 4 5 7,140 4,324 2,360 18,77 25,79q q q q q = − − − − − (113)

e as forças nodais equivalentes (107) ao carregamento dado.

a) Deformada compatível b) Forças equivalentes c) Carregamento

Figura 47: Solução da placa trapezoidal discretizada em dois elementos

Ao analisar esta solução, é essencial distinguir os resultados que têm necessariamente de

satisfazer as condições do problema, mesmo que a solução seja aproximada, dos resultados que

podem violar algumas dessas condições, assim como estimar se esse erro é aceitável ou se exige

o refinamento da solução.

13.1 Incidência Nodal dos Deslocamentos Impostos

Na Figura 48 estão identificados os deslocamentos impostos, 0iq = no exemplo da placa

trapezoidal, e as forças correspondentes, as reacções nodais equivalentes iR .

Esta informação é utilizada para relacionar esses deslocamentos e reacções com os

deslocamentos e as forças nodais de cada elemento, definidos na Figura 30, obtendo-se a Tabela

4 e as seguintes condições de incidência, complementares das definidas pelas equações (106) e

(107):

(2)2 1(2)3 2(1)4 3(1) (2)5 5 4(2)6 5

0

0

0

0

0

d q

d q

d q

d d q

d q

= = = = = = = = = = =

(114)

2,360

25,79 18,77

4,324 7,140

503 15

1 1

2m

x

y

10

20 kPa

0xf = 320yf kNm−= −


61

(2)1 2

(2)2 3

(1)3 4

(1) (2)4 5 5

(2)5 6

R F

R F

R F

R F F

R F

= = = = + =

(115)

Estrutura ( ,i iq R ) 1 2 3 4 5 Elemento 1 ( (1) (1),i id F ) - - 4 5 - Elemento 2 ( (2) (2),i id F ) 2 3 - 5 6

Tabela 4: Incidência nodal (deslocamentos impostos e reacções correspondentes)

Figura 48: Deslocamentos nodais impostos e reacções correspondentes

13.2 Determinação dos Campos de Deslocamento e de Deformação

A informação relativa aos deslocamentos nodais de cada elemento é retirada das equações (106)

e (114) e do resultado (113), obtendo-se (em 810 m− ):

(1) (1)1 1(1) (1)2 2(1) (1)3 3(1) (1)1 4(1) (1)2 5(1) (1)3 6

7,140

4,324

2,360

0

0

18,77

x

x

x

y

y

y

d d

d d

d d

d d

d d

d d

− − − = =

−

(2) (2)1 1(2) (2)2 2(2) (2)3 3(2) (2)4 4(2) (2)1 5(2) (2)2 6(2) (2)3 7(2) (2)4 8

4,324

0

0

2,360

0

0

25,79

18,77

x

x

x

x

y

y

y

y

d d

d d

d d

d d

d d

d d

d d

d d

− − = =

− −

(116)

O campo de deslocamento em cada elemento, medido no referencial global, é determinado

substituindo esta informação na aproximação (35) ou (59) em que o elemento se baseia, tomando

a seguinte expressão quando se usa o resultado (64) para o elemento 1,

(1)

(1)

7,140 2,816 0,982

9,385x

y

u x y

u y

= − + + = −

(117)

1q

2q

3q

4q 5q

5R

4R 3R

1R

2R


62

e o resultado (65) para o elemento 2:

(2)

(2)

(4,324 0,982 ) (2 )

(5,875 3,510 )x

y

u y x

u x y

= − − − = − +

(118)

Para além de verificarem as condições de fronteira,

(1) (2)

(2)

0 0

0 2y y

x

u u para y

u para x

= = = = =

estes resultados confirmam também a continuidade dos deslocamentos entre elementos, de

acordo com a condição (26), com a seguinte variação linear (elementos com dois nós por lado):

(1) (2)

(1) (2)

4,324 0,982 1

9,385 1x x

y y

u u y para x

u u y para x

= = − + = = = − =

Estas verificações são suficientes para assegurar que a solução é cinematicamente admissível,

pois está associada a um campo de deformações calculado impondo a condição de

compatibilidade (14) sobre as soluções (117) e (118) ou, o que é equivalente, aplicando a

aproximação escrita na forma (36) ou (68), encontrando-se (à escala 810− ),

(1)

(1)

(1)

2,816

9,385

0,982

xx

yy

xy

εεγ

= + = − = +

(119)

quando se usa o resultado (69) para o elemento 1, e o resultado (70) para o elemento 2:

(2)

(2)

(2)

4,324 0,982

5,875 3,510

0,982(2 ) 3,510

xx

yy

xy

y

x

x y

εεγ

= + − = − − = + − −

(120)

Pode verificar-se facilmente que o campo de deformações não é contínuo entre elementos, ao

longo da fronteira 1x = . Deverá sê-lo na solução exacta pois não há descontinuidades no

material que forma a placa trapezoidal ( .E const= e .constν = ) nem das forças aplicadas no

domínio, ao longo da fronteira entre elementos. No entanto, a continuidade do campo de

deformações não é uma condição formal de admissibilidade cinemática, pelo que a solução

aproximada que foi determinada é classificada como cinematicamente admissível.

Exercício 25: Sabendo ser a seguinte a solução do sistema de equações (112),

T

1 2 3 4 ( / ) 0,2423 0,3150 1,4575 1,6999q q q q p E= − + − − (121)

trace a deformada da consola quadrada discretizada em dois elementos de 3 nós, verifique as

seguintes expressões para os campos de deslocamento e de deformação,


63

(1)

(1)

0,3150 ( / )

1,6999 ( / )x

y

u x p E

u x p E

= + = −

(2)

(2)

(0,2423 0,5773 ) ( / )

(1,4575 0,2423 ) ( / )x

y

u x y p E

u x y p E

= − − = − +

(1)

(1)

(1)

0,3150 ( / )

0

1,6999 ( / )

xx

yy

xy

p E

p E

εεγ

= + = = −

(2)

(2)

(2)

0,2423 ( / )

0,2424 ( / )

0,9002 ( / )

xx

yy

xy

p E

p E

p E

εεγ

= − = − = −

(122)

e confirme a admissibilidade cinemática da solução.

Comparando a solução definida no exercício anterior com a solução do sistema (111),

T

1 2 3 4 ( / ) 0,9107 1,1115 1,9766 2,6457q q q q p E= − + − − (123)

conclui-se que a solução obtida discretizando a consola quadrada num elemento de 4 nós

(aproximação bilinear do campo de deslocamentos) é, para o mesmo número de graus de

liberdade, 4β = , menos rígida (maiores deslocamentos nodais) do que a solução obtida com

dois elementos de 3 nós (aproximação linear). A solução é também cinematicamente

admissível, sendo as seguintes as definições que se encontram para os campos de

deslocamento e de deformação, o qual varia agora linearmente em todo o domínio da consola:

(1)

(1)

(0,9107 2.0222 ) ( / )

(1,9766 0,6691 ) ( / )x

y

u y x p E

u y x p E

= − − = − +

(1)

(1)

(1)

(0,9107 2,0222 ) ( / )

0,6691 ( / )

(1,9766 2,0222 0,6691 ) ( / )

xx

yy

xy

y p E

x p E

x y p E

εεγ

= − − = − = − − +

(124)

13.3 Determinação dos Campos de Tensão

O campo de esforços é determinado pela relação de elasticidade (15) ou, o que é equivalente,

pela aproximação escrita na forma (38) ou (71). É a seguinte a expressão que se obtém para o

exemplo da placa trapezoidal quando se recorre aos resultados (73) e (74) (em kPa):

(1)

(1)

(1)

0

18,77

0,76

xx

yy

xy

σσσ

= = − =

(2)

(2)

(2)

5,629 2,314 2,158

10,061 7,714 0,647

0,755(2 ) 2,697

xx

yy

xy

x y

x y

x y

σσσ

= + − − = − − − = + − −

(125)

Os resultados correspondentes para a consola quadrada são os seguintes para a discretização

em dois elementos de três nós (aproximação constante),

(1)

(1)

(1)

0,3642

0,1038

0,6538

xx

yy

xy

p

p

p

σσσ

= + = + = −

(2)

(2)

(2)

0,3462

0,3462

0,3462

xx

yy

xy

p

p

p

σσσ

= − = − = −

(126)


64

e para um elemento de 4 nós (aproximação linear):

(1)

(1)

(1)

(1,0001 0,2206 2,2222 )

(0,3002 0,7353 0,6667 )

(0,7602 0,7778 0,2574 )

xx

yy

xy

x y p

x y p

x y p

σσσ

= − + − = − + − = − − +

(127)

O processo de convergência ilustrado na Figura 49, obtido por refinamento da malha de

elementos de 3 nós (refinamento-h), mostra que qualquer uma das soluções acima definidas dá

uma aproximação muito grosseira para o campo de tensões na consola. As quatro soluções aí

apresentadas são obtidas com malhas de 289, 1089 4225 e 16641 elementos, envolvendo 512,

2048, 8192 e 32768 nós, respectivamente.

xxσ yyσ xyσ

Figura 49: Solução da consola quadrada

13.4 Determinação dos Campos de Forças no Domínio e na Fronteira

As forças de massa que equilibram o campo de tensões determinado pelo Método dos Elementos

Finitos, assim como as forças que o equilibram na fronteira dos elementos, não são

imediatamente disponibilizadas pelos programas comerciais.


65

Essas forças são a seguir calculadas para os exemplos em análise, aplicando as definições (77)

e (80), para mostrar que, por regra, os campos de tensão determinados pelo Método dos

Elementos Finitos não são estaticamente admissíveis no sentido forte. Violam tanto as condições

de equilíbrio (13) no domínio dos elementos como as condições de equilíbrio (23) na fronteira

estática da malha. E violam também as condições de equilíbrio (27) de forças entre elementos.

Dadas Calculadas

a) Forças segundo x

Dadas Calculadas

b) Forças segundo y

Figura 50: Solução da consola quadrada discretizada num elemento de 4 nós (1p = )

Para o exemplo da consola quadrada discretizada num elemento de 4 nós, definido na Figura

40, as forças de massa que equilibram o campo de tensões (127) são:

ˆ 0,4779 0

ˆ 1,4444 0

x x

y y

f p f

f p f

= + ≠ =

= − ≠ =

Não sendo conhecidas as forças na fronteira cinemática (as reacções no lado encastrado da

consola), apenas se pode verificar o erro da solução na fronteira estática, encontrando-se os

seguintes resultados:

0xf =

0xt =

0xt =

0xt = ?xt =

ˆ 0,48xf = +

ˆ 1,22xt = −

ˆ 1,00xt = +

1,02− 0,24−

1,00+

1,22−

0,76+ 0,02−

1yt = −

0yt =

0yt =

?yt = 0yf =

ˆ 0,76yt = +

0,24−

ˆ 1,44yf = −

ˆ 1,02yt = + 0,37+

0,37−

0,02+

1,04+ 0,30+


66

ˆ (1,2214 2,2222 ) 0

( 1)ˆ (0,0176 0,2574 ) 0x x

y y

t y p tx

t y p tΓ

= − − ≠ == = + − ≠ =

ˆ (1,0176 0,7778 ) 0

( 1)ˆ (0,3664 0,7353 )x x

y y

t x p ty

t x p t pΓ

= − − ≠ == = + − ≠ = −

ˆ (0,7602 0,7778 ) 0

( 0)ˆ (0,3002 0,7353 ) 0x x

y y

t x p ty

t x p tΓ

= + − ≠ == = + + ≠ =

Estes resultados estão apresentados na Figura 50, estando os resultados correspondentes

obtidos para a discretização em dois elementos de 3 nós resumidos na Figura 51:

ˆ 0,346 0

( 1)ˆ 0,346 0x x

y y

t p tx

t p tΓ

= − ≠ == = − ≠ = (128)

ˆ 0,654 0

( 1)ˆ 0,104x x

y y

t p ty

t p t pΓ

= − ≠ == = + ≠ = − (129)

ˆ 0,346 0

( 0)ˆ 0,346 0x x

y y

t p ty

t p tΓ

= + ≠ == = + ≠ = (130)

Verifica-se que, em ambos os casos, a aproximação é ainda muito fraca pois os sistemas de

forças que equilibram os campos de tensão estão ainda longe de recuperar o carregamento dado

(ver Figuras 38 e 50).

Pelo contrário, e apesar de não se mostrar as escalas de cores para cada componente do campo

de tensão na definição dos campos de tensão apresentada na Figura 49, a análise da solução mais

refinada mostra o seguinte, relativamente às condições de fronteira estáticas do problema,

definidas pelas equações (128) a (130):

• Na face direita (fronteira 1x = ), as componentes de tensão xxσ e xyσ são constantes e

(praticamente) nulas ( 0; 0)x xx y xyt tσ σ= ≈ = ≈ ;

• Na face superior (fronteira 1y = ), a componente de tensão axial yyσ equilibra a carga

aplicada ( )y yyt pσ= ≈ − e a tensão xyσ é nula ( 0)x xyt σ= ≈ ;

• Na face inferior (fronteira 0y = ), as componentes de tensão yyσ e xyσ são nulas

( 0; 0)x xy y yyt tσ σ= − ≈ = − ≈ .

Não é possível concluir sobre a qualidade da solução do campo de tensões na secção de

encastramento (fronteira cinemática, 0x = ). No entanto, regista-se já a esperada concentração de

tensões nos cantos dessa secção.


67

Recordando que o sistema resolvente (46) do Método dos Elementos Finitos não inclui

nenhuma equação que imponha directamente o equilíbrio dos campos de tensão entre elementos,

pois apenas impõe o equilíbrio das forças nodais equivalentes, em ambas as direcções e nos nós

livres da malha, a sequência de resultados apresentada na Figura 49 confirma o que

intuitivamente se espera: à medida que se aumenta o número de nós em que essas condições são

impostas, mais próximo se está de impor as condições de equilíbrio de maneira forte, ou seja, em

todos os pontos do domínio em análise (cerca de 32000 para a solução mais refinada).

Diagrama de corpo livre Reunião dos elementos




Figura 51: Solução da consola quadrada discretizada em dois elementos de 3 nós (1p = )

As descontinuidades que se observam na definição dos campos de tensão apresentada na

Figura 49 para as duas primeiras soluções correspondem, essencialmente, ao desequilíbrio de

forças entre elementos, como se ilustra na Figura 51.

ˆ 0xf =

ˆ 0xf =

0,364 0,707

0,346

ˆ 0,654xt = −

ˆ 0,346xt = + ˆ 0,346xt = +

ˆ 0xf =

ˆ 0xf =

0,364 0,346

0,707

ˆ 0,654xt = −

0,346

ˆ 0,104yt = +

ˆ 0,346yt = +

ˆ 0yf =

ˆ 0yf =

0,654

0,536

ˆ 0,104yt = +

ˆ 0,346yt = +

ˆ 0yf =

ˆ 0yf =

0,346

0,654

0,536


68

O mesmo se conclui quando se analisam os resultados relativos à para a placa trapezoidal,

apresentados na Figura 52. As soluções (125) e (126) mostram que o campo de tensões é

descontínuo entre elementos, para 1x = na placa trapezoidal e para y x= na placa quadrada (ver

Figuras 28 e 38, respectivamente) e não satisfazem a condição mais fraca (mas formalmente

suficiente) de equilíbrio das forças:

(1) (2)

(1) (2)

ˆ ˆ 0ˆ ˆ 0x x

iy y

t t

t tΓ

+ ≠ + ≠

(131)

Dadas Calculadas


Dadas Calculadas


Figura 52: Solução da placa trapezoidal discretizada em dois elementos

13.5 Determinação das Reacções de Apoio

As reacções de apoio, na realidade a aproximação que se obtém com uma dada malha de

elementos finitos, são determinadas calculando as forças que equilibram o campo de tensões ao

longo da fronteira cinemática da estrutura.

?xt =

0xt =

0xt =

0xt =

0xf =

ˆ 0xf = ˆ 5xf =

ˆ 0,76xt = −

0,36 1,00

1,00

3,32

3,32 0,76

5,39 4,64

10

20 kPa

0yt = 0yt =

?yt =

20yf = −

0,76

ˆ 18,77yt = +

9,07

0,76

4,64

19,07 26,78

5,04

17,76 25,47

ˆ 0yf = ˆ 1,4yf = +


69

São disso exemplo as forças na secção de encastramento da consola quadrada definidas nas

Figura 50 e 51 para as discretizações num elemento de 4 nós e em dois elementos de 3 nós.

Relativamente ao exemplo da placa trapezoidal, essas reacções correspondem à força yt ao longo

do bordo 0y = e à força xt ao longo do bordo 2x = , representadas na Figura 52.

Os programas comerciais oferecem usualmente a opção de listar as forças nodais na malha de

elementos finitos, designadamente as forças equivalentes nos nós restringidos (geralmente

designadas por reacções de apoio) recorrendo às condições de incidência nodal, como se ilustrou

com a relação (115) para o exemplo da placa trapezoidal. Essa informação pode ser usada para

confirmar o equilíbrio global das forças aplicadas.

O processo mais prático para calcular as forças nodais equivalentes consiste em recorrer à

equação elementar (39), dado que se conhece a matriz de rigidez de cada elemento e, após a

solução do problema, os deslocamentos nodais da malha e, portanto, de cada elemento.

Os resultados para o exemplo da consola trapezoidal, resumidos na Figura 53, são obtidos

recorrendo às definições (98) e (99) das matrizes de rigidez elementares e à solução (116) para

os deslocamentos nodais. Pode verificar-se facilmente que as forças aplicadas a cada elemento e

à malha satisfazem as condições de equilíbrio global: são nulas as resultantes das forças segundo

as direcções x e y, assim como o momento resultante na direcção ortogonal ao plano da estrutura.


Figura 53: Forças nodais equivalentes para o exemplo da placa trapezoidal

Exercício 26: Determine as forças nodais equivalentes para o exemplo da consola quadrada

discretizada em dois elementos de 3 nós e num elemento de 4 nós e verifique as condições de

equilíbrio global. Utilize os resultados (102) a (104) e as soluções (121) e (123), para além

das condições de incidência nodal.

Exercício 27: Identifique a razão que assegura que as soluções obtidas pelo Método dos

Elementos Finitos sejam globalmente equilibradas.

1

3 4

1

3

2 2

0,38 0,38

9,38 7,29

0,38 0,38 0

0,75 11,94 10,35

1,94

1,94

15,0

10,13

15 16,67

0 0

10,35

1,94

1,94 0

0,75 21,97


70

14. Elementos Isoparamétricos

Ao desenhar malhas de elementos finitos, é necessário utilizar elementos com geometria

irregular, por exemplo, para fazer a transição entre elementos com dimensões diferentes em

zonas com diferentes graus de refinamento. É também importante dispor de elementos com lados

curvos, para representar adequadamente a geometria da estrutura.

Surgem três tipos de dificuldades quando se tenta definir esse tipo de elemento no sistema de

coordenadas global da malha, ou mesmo quando se tenta defini-los em coordenadas locais,

obtidas por translação e/ou rotação do referencial global.

Duas dessas dificuldades têm consequências importantes em termos da eficácia

computacional do método. A primeira é a complexidade das expressões que se obtêm para as

funções de aproximação, sujeitas às condições (57) e (58), a qual cresce substancialmente

quando se aumenta o grau da aproximação e se generaliza a geometria do elemento. Por outro

lado, torna-se inviável definir de um uma maneira computacionalmente eficaz os limites dos

integrais (40), (42) e (43) que definem cada um dos coeficientes da matriz de rigidez e dos

vectores das forças nodais equivalentes presentes na equação resolvente.

No entanto, é a terceira dificuldade que justifica o abandono da ideia, usada nas secções

anteriores, de definir as funções de aproximação no referencial global da malha ou num

referencial local do elemento. Como se ilustra na Figura 54, quando se definem no referencial

global da malha ( , )x y funções de aproximação de grau superior e se impõem as condições (57)

e (58), essas condições deixam de ser suficientes para garantir que as funções são coincidentes

ao longo do lado partilhado pelos elementos. Deixa de estar assegurada a continuidade dos

deslocamentos entre elementos, uma condição que está na base da formulação do Método dos

Elementos Finitos aqui usada.

a) Elementos adjacentes b) Variação das funções no lado comum

Figura 54: Descontinuidade de funções de aproximação definidas no referencial global

Estas dificuldades são ultrapassadas, com alguns custos computacionais, definindo elementos

com geometria simples (elementos-mestre), sobre os quais se definem de uma maneira fácil e

x

y

0s =

s L=

is i

s

0s= s L= is

s

1

(1)( , )i x yΨ

(2)( , )i x yΨ


71

sistemática as funções de aproximação sujeitas às condições (57) e (58). Como a geometria é

simples, também o é a definição dos limites dos integrais exigidos pela aplicação do método de

cálculo. Nas aplicações bidimensionais são usados dois tipos de elemento-mestre, o triângulo

rectângulo de lado unitário e o elemento quadrado com semi-lado unitário, representados na

Figura 55.

Figura 55: Elementos-mestre triangular e quadrangular

Para tornar essas definições válidas para elementos com geometria arbitrária, o elemento-

mestre e as funções de aproximação são caracterizados num espaço independente do espaço em

que a malha é definida (espaço das coordenadas naturais). O elemento-mestre é depois mapeado

sobre cada elemento da malha, por exemplo nos representados na Figura 54a), recorrendo a uma

mudança de coordenadas, do sistema de coordenadas naturais ( , )ξ η para o sistema de

coordenadas da malha de elementos finitos ( , )x y . É esta a fonte dos custos computacionais

adicionais que resultam desta opção para resolver as dificuldades acima referidas.

Em princípio, é possível usar diferentes expressões para as funções de aproximação e de

mapeamento do elemento-mestre. Um elemento é definido como isoparamétrico quando se usam

as mesmas funções para os dois fins, sendo implicitamente garantida a continuidade dos

deslocamentos entre elementos com os mesmos nós nos lados que partilham.

14.1 Elemento-Mestre

É sobre estes elementos que são definidas as funções de aproximação. Essas funções tomam

expressões relativamente simples, as quais são deduzidas uma única vez para cada grau de

aproximação.

Mantêm-se válidos os critérios definidos na Secção 9, designadamente as condições (57) e

(58), assim como a interpretação aí dada sobre o número total de nós e o número de nós por

fronteira para definir o grau das aproximações no domínio e na fronteira dos elementos-mestre.

ξ

η

1

1

1ξ η+ = 0ξ =

0η =

ξ

η

1 1

1

1

1η = +

1η = −

1ξ = − 1ξ = +


72

Para ilustrar a simplicidade relativa das expressões dessas funções, resume-se na Tabela 5 as

funções de aproximação para elementos triangulares de 3 nós (aproximação linear) e de 6 nós

(aproximação quadrática) representados na Figura 26. As funções de aproximação para os

elementos quadrangulares de 4 nós (aproximação quadrática incompleta) e de 9 nós

(aproximação quártica incompleta) representados na Figura 27 estão definidas na Tabela 6.

Elemento 3 nós 6 nós

1Ψ 1 ξ η− − (1 ) (1 2 2 )ξ η ξ η− − − −

2Ψ ξ 4 (1 )ξ ξ η− −

3Ψ η (2 1)ξ ξ −

4Ψ 4ξ η

5Ψ (2 1)η η −

6Ψ 4 (1 )η ξ η− −

Tabela 5: Elementos triangulares


1Ψ 14 (1 ) (1 )ξ η− − 1

4 (1 ) (1 )ξ η ξ η− −

2Ψ 14 (1 ) (1 )ξ η+ − 21

2 (1 ) (1 )η ξ η− − −

3Ψ 14 (1 ) (1 )ξ η+ + 1

4 (1 ) (1 )ξ η ξ η+ −

4Ψ 14 (1 ) (1 )ξ η− + 21

2 (1 ) (1 )ξ ξ η+ −

5Ψ 14 (1 ) (1 )ξ η ξ η+ +

6Ψ 212 (1 ) (1 )η ξ η− +

7Ψ 14 (1 ) (1 )ξ η ξ η− − +

8Ψ 212 (1 ) (1 )ξ ξ η− − −

9Ψ 2 2(1 ) (1 )ξ η− −

Tabela 6: Elementos quadrangulares

a) Função 1Ψ b) Função 2Ψ c) Função 9Ψ

Figura 56: Funções de aproximação do elemento quadrangular de 9 nós

-1

-0.5

0

0.5

1 -1

-0.5

0

0.5

1

-0.1

0

0.1

-1

-0.5

0

0.5

-1

-0.5

0

0.5

1-1

-0.5

0

0.5

1

00.250.5

0.751

-1

-0.5

0

0.5

-1

-0.5

0

0.5

1-1

-0.5

0

0.5

1

00.25

0.5

0.75

1

-1

-0.5

0

0.5


73

Os gráficos representados nas Tabelas 1 e 2 mantêm-se válidos para as funções de

aproximação dos elementos de 3 e 4 nós, sendo agora definidos no sistema de coordenadas

naturais ( , )ξ η . Na Figura 56 representam-se três funções típicas do elemento de 9 nós. Pode

confirmar-se que as funções variam quadraticamente (três nós por lado) ao longo dos lados que

contêm o nó onde a função toma valor unitário.

Exercício 26: Trace as funções de aproximação do elemento triangular de 6 nós.

14.2 Mapeamento do Elemento-Mestre

As funções de aproximação formuladas para cada elemento-mestre são usadas para mapear o

elemento em qualquer elemento da malha de elementos finitos (daí a designação alternativa de

funções de forma). Essa mudança de coordenadas é escrita de forma análoga à usada para

aproximar os deslocamentos, definida pela equação (59),

1

1

( , )

( , )

N

i ixi

N

i iyi

x c

y c

Ψ ξ η

Ψ ξ η

=

=

=

=

∑

∑ (132)

em que ixc e iyc representam agora as coordenadas em x e em y, medidas no referencial global da

malha, do nó i do elemento com N nós.

Figura 57: Mapeamento do elemento de 4 nós num elemento trapezoidal

Exemplifica-se na Figura 57 o mapeamento do elemento-mestre de 4 nós num elemento com

o domínio da placa trapezoidal, Figura 28. Depois de determinar as coordenadas nodais,

1 2 3 4 0 2 2 1x x x xc c c c =

1 2 3 4 0 0 2 2y y y yc c c c =

e aplicar as expressões definidas na Tabela 6, obtém-se a definição da mudança de coordenadas:

1 2

3

ξ

η

4

1 1

1

1

1 1

2m

x 1 2

3 4

y ξ

η


74

41

1 2 3 4 41

4

1 2 3 41

( , ) (0) (2) (2) (1) 1 (1 ) (1 )

( , ) (0) (0) (2) (2) 1

i ixi

i iyi

x c

y c

Ψ ξ η Ψ Ψ Ψ Ψ ξ ξ η

Ψ ξ η Ψ Ψ Ψ Ψ η

=

=

= = ⋅ + ⋅ + ⋅ + ⋅ = + + − +

= = ⋅ + ⋅ + ⋅ + ⋅ = +

∑

∑ (133)

É importante frisar que a mudança de coordenadas (132) não é, em geral, invertível. Ou seja,

não é em geral possível determinar univocamente as coordenadas ( , )ξ η no elemento-mestre de

um dado ponto do elemento com coordenadas ( , )x y no referencial global da malha.

Como adiante se mostra, esta limitação tem consequências que complicam a implementação

dos elementos isoparamétricos, pois impede que as funções de aproximação sejam expressas no

sistema de coordenadas globais. No entanto, continua a ser computacionalmente mais eficaz

basear o cálculo em coordenadas naturais mesmo quando a mudança de coordenadas é invertível,

por exemplo nos elementos de 3 nós e nos elementos de 4 nós com pelo menos dois lados

paralelos após o mapeamento.

Por exemplo, a mudança de coordenadas (133) tem a inversa,

1

4 ; 14

xy

yξ η−= = −

− (134)

a qual em nada simplifica a definição das funções de aproximação do elemento de 4 nós

definidas na Tabela 6, nem a manipulação dessas funções, como adiante se mostra.

Uma outra característica importante, mas esta positiva, é que a mudança de coordenadas (132)

permite representar domínios curvos. Como as funções de aproximação são definidas de modo a

ter grau 1n− num lado com n nós, o mapeamento permite que um lado de um elemento-mestre

com esse número de nós se transforme numa linha com grau menor ou igual a 1n− no

referencial da malha.

Uma consequência directa desta conclusão é o caso particular dos elementos triangulares de 3

nós e quadrangulares de 4 nós (ambos com dois nós por lado), os quais só podem ser mapeados

em polígonos de três e quatro lados (rectos), respectivamente.

Exercício 28: Considere o mapeamento do elemento-mestre de 3 nós definido na Figura 26

no elemento triangular de 3 nós definido na Figura 23. Defina e mostre que é invertível a

mudança de coordenadas (132) usando os resultados resumidos na Tabela 5. Mostre que essa

mudança de coordenadas é invertível e use o resultado obtido para recuperar as definições

para as funções de aproximação no sistema de coordenadas globais dadas na Tabela 1.

Um exercício análogo, mas que exige menos esforço de cálculo, consiste em mapear o

elemento-mestre de 4 nós definido na Figura 27 no elemento rectangular de 4 nós definido na


75

Figura 24 e usar os resultados resumidos nas Tabelas 2 e 6. O exercício seguinte serve para

ilustrar a representação de fronteiras não lineares e a aproximação implícita na mudança de

coordenadas:

Exercício 29: Defina o mapeamento do elemento-mestre de 6 nós definido na Figura 26 no

quarto de círculo representado na Figura 58 e trace sobre o arco de círculo 2 2( 1)x y+ = a

aproximação da geometria obtida (mapeamento do lado 1ξ η+ = do elemento-mestre).

Figura 58: Mapeamento do elemento de 6 nós num quarto de círculo

14.3 Aproximação do Campo de Deslocamentos

A aproximação dos deslocamentos continua a ser expressa na forma (59), com a diferença que as

funções de aproximação estão agora definidas no sistema de coordenadas naturais:

1

1

( , ) ( , )

( , ) ( , )

N

x i ixi

N

y i iyi

u x y d

u x y d

Ψ ξ η

Ψ ξ η

=

=

=

=

∑

∑ (135)

Para além de se continuar a aproximar independentemente as duas componentes (de

translação) do movimento, mantém-se a hipótese que tanto essas componentes como os

deslocamentos nodais são medidos no referencial da malha de elementos finitos, como se mostra

na Figura 29 para os elementos de 3 e 4 nós. Consequentemente, mantém-se também a

identificação das forças nodais correspondentes (ver Figura 33 para os mesmos exemplos).

Retoma-se na Figura 59 o exemplo da placa trapezoidal admitindo agora que é discretizada

num único elemento isoparamétrico de 4 nós, de acordo com o mapeamento (133), ilustrado na

Figura 56. Para a numeração dos deslocamentos nodais definida na Figura 58 e para as funções

definidas na Tabela 6, a aproximação (135) tem a seguinte expressão:

1 1 4 2

4 3 3 4

( , ) ( , ) ( , )

( , ) ( , ) ( , )x

y

u x y q q

u x y q q

Ψ ξ η Ψ ξ ηΨ ξ η Ψ ξ η

= += +

(136)

1

1

1

2 3

ξ

η

4

5

6

1

1

1

2 3

5

6

4

x

y


76

a) Deslocamentos e forças nodais b) Efeito do deslocamento nodal 2 1q =

Figura 59: Discretização da malha trapezoidal num elemento de 4 nós

Se a transformação de coordenadas não for invertível (e é invertível neste exemplo), a

aproximação (136) mostra que não é possível determinar o deslocamento num ponto de

coordenadas ( , )x y do elemento inserido na malha, mesmo depois de determinar os

deslocamentos nodais iq .

Só é possível escolher um ponto do elemento mestre, por exemplo com coordenadas

( , ) (0,0)ξ η = , determinar o ponto correspondente no elemento inserido na malha aplicando a

mudança de coordenadas (132), a que corresponderia neste exemplo o ponto 14( , ) (1 ,1)x y = + , e

calcular aí o deslocamento através da aproximação (135), ficando para este exemplo:

1 1 1

1 1 4 2 1 24 4 4

1 1 14 3 3 4 3 44 4 4

(1 ,1) (0,0) (0,0)

(1 ,1) (0,0) (0,0)x

y

u q q q q

u q q q q

Ψ ΨΨ Ψ

+ = + = ++ = + = +

14.4 Aproximação dos Campos de Deformação e de Tensão

As deformações compatíveis com a aproximação (135) do campo de deslocamento são definidas

pela condição (10), ou (14) na forma explícita. Como as funções de aproximação deixam de estar

expressas em coordenadas globais, as opções possíveis são recorrer à inversão da transformação

de coordenadas ou recorrer à derivada da função implícita.

Mesmo quando a mudança de coordenadas tem inversa, a expressão que as funções de

aproximação tomam pode não facilitar o cálculo das derivadas. Para o exemplo do mapeamento

(133) do elemento-mestre de 4 nós num elemento trapezoidal, a inversão (134) da mudança de

coordenadas origina a seguinte expressão para a primeira das funções de aproximação:

1

1( , ) 1

4

xx y y

yΨ −= − −

Por isso, esta via é abandonada em favor do processo mais geral, e mais sistemático, de

calcular as derivadas através das expressões, como se mostra adiante:

1 1,q Q

2 2,q Q

3 3,q Q

4 4,q Q

1 2

3 4

2 1q =

x

y

4xu Ψ=

0yu =


77

( , )

( , )

i i ix i

i i iy i

x x x

y y y

Ψ Ψ Ψξ ηΨ ξ ηξ η

Ψ Ψ Ψξ ηΨ ξ ηξ η

∂ ∂ ∂∂ ∂∂ = = +∂ ∂ ∂ ∂ ∂

∂ ∂ ∂∂ ∂∂ = = +∂ ∂ ∂ ∂ ∂

(137)

A aproximação do campo de deformações continua a ser expressa nas formas (36) ou (68). A

correspondente aproximação do campo de tensões é determinada impondo a relação de

elasticidade (9), para estados planos de tensão (15) ou de deformação (16), e expressa nas formas

alternativas (38) ou (71).

14.5 Processamento da Mudança de Coordenadas

Para sistematizar os cálculos exigidos pela mudança de coordenadas (132), define-se a matriz

Jacobiana da transformação,

x y

x y

ξ ξ

η η

∂ ∂ ∂ ∂ =

∂ ∂ ∂ ∂

J (138)

determina-se o seu determinante, o Jacobiano da transformação,

| | 0x y x y

Jξ η η ξ

∂ ∂ ∂ ∂= = − ≠∂ ∂ ∂ ∂

J

e a inversa da matriz Jacobiana, para obter os termos presentes na definição (137):

1 1

y y

Jx x

η ξ

η ξ

− −

∂ ∂ − ∂ ∂ =

∂ ∂ − ∂ ∂

J (139)

Sendo a seguinte a matriz Jacobiana da transformação inversa,

1 x x

y y

ξ η

ξ η−

∂ ∂ ∂ ∂ =

∂ ∂ ∂ ∂

J (140)

obtém-se a informação necessária para calcular as derivadas parciais (137),

1 1

1 1

;

;

y yJ J

x x

x xJ J

y y

ξ ηη ξ

ξ ηη ξ

− −

− −

∂ ∂ ∂ ∂+ −∂ ∂ ∂ ∂∂ ∂ ∂ ∂− +∂ ∂ ∂ ∂

= =

= = (141)

em que, para a mudança de coordenadas (132):


78

1 1

1 1

;

;

N Ni i

ix ixi i

N Ni i

iy iyi i

x xc c

y yc c

Ψ Ψξ ξ η η

Ψ Ψξ ξ η η

= =

= =

∂ ∂∂ ∂= =∂ ∂ ∂ ∂

∂ ∂∂ ∂= =∂ ∂ ∂ ∂

∑ ∑

∑ ∑

O cálculo das derivadas em coordenadas globais (137) não é simples, mas é facilmente

programável depois de tabelar as derivadas das funções nas coordenadas naturais.

São os seguintes os resultados que se obtêm para o exemplo da Figura 57, ver Eq. (133):

3 01

1 14

ηξ

− = −

J

14 (3 )J η= − (142)

1 1 04

1 33 ξ ηη−

= − −− J (143)

Para este exemplo, os resultados (141) e (143) mostram serem as seguintes as expressões das

derivadas (137) das funções de aproximação:

1

1

(1) (0)

( 1) (3 )

i i i

i i i

Jx

Jy

Ψ Ψ Ψξ η

Ψ Ψ Ψξ ηξ η

−

−

∂ ∂ ∂= + ∂ ∂ ∂

∂ ∂ ∂= − + − ∂ ∂ ∂

(144)

Para assegurar a estabilidade numérica da mudança de coordenadas, é importante garantir que

o Jacobiano da transformação seja positivo em qualquer ponto do domínio do elemento:

0J > (145)

Esta condição é violada em duas situações típicas. A primeira, ilustrada na Figura 60,

acontece quando não se respeita a sequência da numeração dos nós ao definir o mapeamento

(132) do elemento-mestre. A segunda, ilustrada na Figura 61 para o elemento de 4 nós, sucede

quando o elemento-mestre é mapeado num elemento não convexo, isto é, quando existe um

segmento definido por dois pontos da fronteira que não é interior ao domínio do elemento.

No primeiro caso o Jacobiano é negativo em todos os pontos do domínio, pois tem-se:

; ;x a y b J abη ξ= = = −

No segundo caso obtém-se,

18 [ 2 ( ) ( ) ]J a b b a ξ η= + − +

concluindo-se que o Jacobiano é positivo em qualquer ponto do domínio para 12b a> (elemento

convexo), nulo no nó 3 se 12b a= e negativo ou nulo em todos os pontos com coordenadas,


79

2

( )b

a bξ η+ ≤

−

se 12b a< (elemento não convexo). O resultado (142), obtido para o mapeamento convexo

definido na Figura 57, mostra que o Jacobiano é positivo para qualquer ponto do domínio, ou

seja 1 1ξ− ≤ ≤ + e 1 1η− ≤ ≤ + em todos os elementos isoparamétricos quadrangulares.

Figura 60: Alteração da sequência de identificação dos nós

Figura 61: Convexidade do mapeamento

Um terceiro aspecto que interessa sublinhar é que o valor do Jacobiano, calculado em

qualquer ponto do domínio, não deve ter valores próximos de zero, pois instabiliza o cálculo da

inversa da matriz Jacobiana (139). Como se mostra a seguir, estas situações ocorrem quando o

mapeamento define elementos muito distorcidos, com áreas tendencialmente nulas.

14.6 Cálculo de Integrais

Para definir a equação resolvente (39) de cada elemento da malha é necessário calcular integrais

de domínio, da forma,

( , )I f dΩ Ωξ η Ω= ∫ (146)

por exemplo na determinação da matriz de rigidez (40) e do vector das forças nodais

equivalentes às forças de massa (42) ou a campos de tensão ou de deformação residuais (47).

É ainda necessário calcular integrais de fronteira, da forma,

1 2

3

ξ

η

1

2

3

x

y

a

b 1

1

1 2

3

ξ

η

4

1 1

1

1

x 1 2

3

4

y

12 a

12 a

12 a 1

2 a

( , )b b

12b a>

ξ η α+ <


80

( , )I f dΓ ΓΓξ η Γ= ∫ (147)

por exemplo no cálculo do vector das forças nodais equivalentes às forças de fronteira (43) ou da

contribuição (48) de apoios elásticos para a matriz de rigidez do elemento. Na equação (147), a

notação ( , )f Γξ η representa o valor da função na fronteira em que se realiza a integração.

Uma grande vantagem dos elementos isoparamétricos é que esses integrais não têm de ser

calculados nos elementos mapeados, em que, por norma, seria muito difícil (e

computacionalmente muito ineficaz) definir os limites de integração. São calculados no

elemento-mestre, em que é simples a definição desses limites.

Para calcular integrais de domínio (146) em elementos-mestre quadrangulares,

1 1

1 1( , ) ( , )I f J d dΩ ξ η ξ η ξ η

+ +

− −= ∫ ∫ (148)

impõe-se a relação entre elementos de área no elemento mapeado e no elemento-mestre, a qual é

definida pelo Jacobiano da transformação:

dx dy J d dξ η=

Para definir o Jacobiano da transformação de um lado do elemento mapeado,

2 2ds dx dy= +

x x

dx d dξ ηξ η

∂ ∂= +∂ ∂

y y

dy d dξ ηξ η

∂ ∂= +∂ ∂

explora-se o facto desse lado ter definições simples no elemento-mestre. Os lados de um

elemento quadrangular são definidos pelas equações 1ξ ξ= = ± , com 0dξ = , ou 1η η= = ± ,

com 0dη = , como se pode verificar na Figura 27. O integral de fronteira (147) toma as seguintes

expressões para os lados em que ξ ξ= e η η= , respectivamente:

1

1( , ) ( , )I f J dΓξ ξ η ξ η η

+

−= ∫ (149)

1

1( , ) ( , )I f J dη Γξ η ξ η ξ

+

−= ∫ (150)

2 2

( , )x y

JΓξ ξ

ξ ηη η

∂ ∂= + ∂ ∂

2 2

( , )x y

JΓη η

ξ ηξ ξ

∂ ∂= + ∂ ∂


81

Os integrais de domínio (148) e de fronteira (149) e (150) são calculados numericamente,

aplicando-se normalmente a regra de quadratura de Gauss-Legendre,

1 1

( , ) ( , )M M

i j i j i ji j

I f J w wΩ ξ η ξ η= =

=∑∑

1

( , ) ( , )M

i i ii

I f J wΓξ ξ η ξ η=

=∑

1

( , ) ( , )M

i i ii

I f J wη Γξ η ξ η=

=∑

em que iξ e jη definem os M pontos de integração e iw e jw os pesos correspondentes (os

quais estão definidos no texto sobre Análise de Estruturas Articuladas).

O número de pontos de integração é escolhido atendendo ao grau das funções integrandas,

admitindo que são polinomiais. Tal só sucede para o caso particular de mudanças de coordenadas

isomorfas, em que o Jacobiano é constante.

Quando isso não acontece, o caso mais frequente, as funções integrandas não são polinomiais,

sendo consequentemente aproximado o resultado dessas integrações. Esse facto é ilustrado pelo

resultado (144), lembrando-se que essas derivadas intervêm na definição da matriz (37) de

aproximação do campo de deformação, presente, por exemplo, no cálculo da matriz de rigidez

(40) do elemento.

No caso de elementos triangulares o integral de domínio (148) toma a seguinte expressão (ver

Figura 26),

1 1

0 0( , ) ( , )I f J d d

ξ

Ω ξ η ξ η η ξ−

= ∫ ∫ (151)

sendo as expressões (149) e (150) válidas para os lados do elemento-mestre triangular definidos

pelas condições 0ξ ξ= = e 0η η= = , com limites de integração [ 0,1]:

1

0(0, ) (0, )I f J dΓξ η η η= ∫ (152)

1

0( ,0) ( ,0)I f J dη Γξ ξ ξ= ∫ (153)

No lado 1η ξ= − esse integral toma a seguinte expressão:

1

0( ,1 ) ( ,1 )I f J dη Γξ ξ ξ ξ ξ

+= − −∫ (154)

2 2

( , ) 1x x y y

J paraΓη η

ξ η η ξξ η ξ η

∂ ∂ ∂ ∂= − + − = − ∂ ∂ ∂ ∂


82

Estes integrais são também calculados numericamente. Existem diferentes regras de

quadratura para domínios triangulares, aplicáveis ao cálculo dos integrais (151). Os integrais de

fronteira (152) a (154) são ainda calculados pela regra de quadratura de Gauss-Legendre.

15. Traçado de Malhas de Discretização

Os programas de elementos finitos incluem rotinas que permitem decompor domínios com

topologias e geometrias complexas em subdomínio mais simples, os quais são discretizados em

malhas de elementos (de facto, células) triangulares ou quadrangulares.

Essas rotinas baseiam-se em considerações topológicas e procuram gerar malhas regulares,

formadas por elementos com formas e dimensões semelhantes, sem transições abruptas em

forma e em dimensões. Todavia, esse objectivo nem sempre é adequadamente atingido, pois a

lógica em que se baseiam não é facilmente generalizável.

Para além disso, essas rotinas são geralmente desenvolvidas sem atender ao problema que vai

ser dessa maneira discretizado, nem à técnica de solução adoptada na sua solução. Ou seja, e no

presente contexto, as malhas são traçadas sem atender às zonas do domínio estrutural onde se

esperam grandes gradientes (concentrações) de tensões, onde é necessário refinar a discretização.

Para além disso, não reconhecem os potenciais problemas numéricos associados à mudança de

coordenadas quando se usam elementos isoparamétricos. É por isso importante analisar a malha

gerada automaticamente e verificar se essas condições estão adequadamente asseguradas.

Figura 62: Esquemas de transição no refinamento de malhas de elementos finitos

Os problemas associados à instabilização da mudança de coordenadas (Jacobianos negativos

ou próximos de zero) são evitados rejeitando (ou corrigindo, se o programa o permitir) malhas

com elementos não convexos ou excessivamente distorcidos. Alguns programas permitem

refinar localmente as malhas, ilustrando-se na Figura 62 esquemas de transição usados no

refinamento de subdomínios triangulares e quadrangulares em elementos triangulares e

quadrangulares.

Na Figura 63a) ilustra-se uma primeira decomposição do domínio na análise de um provete de

betão fendilhado, simplesmente apoiado e sujeito a uma força concentrada. Na Figura 63b)


83

ilustra-se um primeiro refinamento da malha, com o qual não é ainda possível recuperar os

campos de tensão ilustrados na Figura 64, os quais representam as concentrações de tensões que

se desenvolvem na vizinhança da ponta da fenda e das cargas pontuais. Ilustra-se também na

Figura 64b) o forte efeito de descentrar a posição da fenda.

a) Decomposição do domínio

b) Refinamento da malha

Figura 63: Discretização do modelo de provete de betão

Distribuição da tensão σxx

Distribuição da tensão σyy

Distribuição da tensão σxy

a) Fenda a meio-vão b) Fenda descentrada

Figura 64: Campos de tensão no provete de betão fendilhado

P

120 120

80mm

20

A


84

Para facilitar o refinamento local das malhas, alguns programas disponibilizam elementos de

transição, elementos que compatibilizam os graus de aproximação usados em dois ou mais

elementos com diferentes graus de aproximação, ou seja com um número diferente de nós por

lado ou face. Por exemplo, um elemento com dois nós por lado (variação linear) é

compatibilizado com um elemento com três nós no lado comum (variação quadrática)

restringindo os deslocamentos no nó intermédio de modo a ser linear a variação do deslocamento

nesse lado do elemento. Para o exemplo da Figura 65, esta técnica consiste em obrigar os

deslocamentos do nó j a tomar o valor médio dos deslocamentos dos nós de extremidade do lado,

os nós i e k.

Figura 65: Transição entre elementos com dois e três nós por lado

16. Elasticidade Tridimensional

A formulação do Método dos Elementos Finitos anteriormente apresentada para a análise de

problemas de Elasticidade Plana (Estados Planos de Tensão e Estados Planos de Deformação) é

directamente generalizável à solução de problemas tridimensionais. Os principais aspectos dessa

generalização para a análise de sólidos elásticos são a seguir resumidos.

16.1 Variáveis

As variáveis continuam a ser organizadas em pares duais (variáveis estáticas e variáveis

cinemáticas correspondentes), de modo a preservar as definições dadas na Secção 4 para o

trabalho realizado pelas forças interiores e exteriores. De acordo com a notação aí usada, Ω

representa o domínio (volume) do corpo tridimensional e Γ a fronteira (superfície).

Os campos de tensão e de deformação são agora definidos por seis componentes

independentes das nove componentes dos tensores (simétricos) de tensão e de deformação,

Txx yy zz xy xz yzσ σ σ σ σ σ=s

(155)

Txx yy zz xy xz yzε ε ε γ γ γ=e

(156)

sendo três as componentes necessárias para caracterizar as forças de massa e de fronteira, assim

como os deslocamentos correspondentes:

i

j

k


85

Tx y zf f f=f

(157)

Tx y zt t t=t

(158)

Tx y zu u u=u

(159)

Todas as variáveis são medidas no referencial global da malha de elementos finitos.

Figura 66: Convenção adoptada na medição de variáveis

16.2 Equações

A forma explícita das equações de equilíbrio (8) e de compatibilidade (10) é a seguinte:

0

0

0

x xx y xy z xz x

y yy x xy z yz y

z zz x xz y yz z

f

f em

f

σ σ σσ σ σ Ωσ σ σ

∂ + ∂ + ∂ + = ∂ + ∂ + ∂ + = ∂ + ∂ + ∂ + =

xx x x

yy y y

zz z z

xy y x x y

xz z x x z

yz z y y z

u

u

uem

u u

u u

u u

εεε

Ωγγγ

= ∂ = ∂ = ∂ = ∂ + ∂ = ∂ + ∂ = ∂ + ∂

(160)

Quando estas relações são escritas matricialmente, de acordo com a organização

anteriormente definida para os vectores de tensão (155) e de deformação (156), e para os

vectores de forças (157) e de deslocamentos (159), obtendo-se a seguinte generalização para os

operadores diferenciais de compatibilidade e de equilíbrio,

0 0

0 0

0 0

0

0

0

x

y

z

y x

z x

z y

∂ ∂ ∂

= ∂ ∂ ∂ ∂

∂ ∂

A

dx

dy

dz

xxσ

xyσ

xzσ

yyσ

yxσ

yzσ

zzσ

zxσ

zyσ

x y

z ,x xt u

,y yt u

,z zt u FronteiraΓ

DomínioΩ


86

T

0 0 0

0 0 0

0 0 0

x y z

y x z

z x y

∂ ∂ ∂ = ∂ ∂ ∂ ∂ ∂ ∂

A (161)

mantendo-se a relação que caracteriza a análise baseada na hipótese de linearidade geométrica.

A forma explícita das condições de equilíbrio (11) na fronteira estática e de compatibilidade

(12) na fronteira cinemática é a seguinte,

x xx y xy z xz x

y yy x xy z yz y t

z zz x xz y yz z

n n n t

n n n t em

n n n t

σ σ σσ σ σ Γσ σ σ

+ + = + + = + + =

x x

y y u

z z

u u

u u em

u u

Γ = = =

em que xn , yn e yn continuam a representar as componentes da normal exterior unitária à

fronteira, no ponto da superfície em que estão aplicadas a força com componente xt , yt e zt ,

sendo xu , yu e zu as componentes impostas do vector de deslocamento.

Na formulação matricial da condição de equilíbrio (11) na fronteira, a matriz de equilíbrio

continua a ser análoga à matriz de equilíbrio no domínio, substituindo a componente da normal

unitária exterior o operador diferencial correspondente (161):

T

0 0 0

0 0 0

0 0 0

x y z

y x z

z x y

n n n

n n n

n n n

=

N

A matriz de rigidez local, presente na relação de elasticidade (9), toma a seguinte expressão

para materiais homogéneos e isotrópicos,

2 0 0 0

2 0 0 0

2 0 0 0

0 0 0 0 0

0 0 0 0 0

0 0 0 0 0

λ µ λ λλ λ µ λλ λ λ µ

µµ

µ

+ + +

=

D (162)

de acordo com a ordenação dos vectores de tensão (155) e de deformação (156), em que λ e µ

são as constantes de Lamé:


87

(1 ) (1 2 )

Eνλν ν

=+ −

2(1 )

EGµ

ν= =

+

16.3 Funções de Aproximação

O cálculo continua a ser baseado em elementos isoparamétricos, sendo agora os elementos-

mestre o tetraedro recto, com aresta unitária, e o hexaedro cúbico, com semi-aresta unitária,

como se ilustra na Figura 67. As funções de aproximação continuam a ser polinómios

construídos de maneira a satisfazer as condições (57) e (58):

1

( , , )0i j j j

se j i

se j iΨ ξ η ζ

= ≠

=

1

( , , ) 1N

ii

Ψ ξ η ζ=

=∑

a) Elemento tetraédrico b) Elemento hexaédrico

Figura 67: Elementos-mestre tridimensionais

Figura 68: Pirâmide de Pascal

ξ η

ζ

1 1

1

ξ

ζ

η

2

2

2

ξη

1

ξ η

ζ 2ξ

3ξ

3ζ

2ζ 2η

3η

ηζ

ξζ

2ξ η

2ξη 2η ζ

2ηζ

2ξζ


88

O número de nós do elemento está relacionado com o grau de aproximação, sendo ainda

válidas as identificações apresentadas na Secção 9, recorrendo agora à pirâmide de Pascal

representada na Figura 68, definida em termos das coordenadas naturais ( , , )ξ η ζ .

As funções de aproximação do elemento tetraédrico de 4 nós (linear) e do elemento

hexaédrico de 8 nós (trilinear) representados na Figura 69 são as resumidas na Tabela 7.


1Ψ 1 ξ η ζ− − − 18 (1 ) (1 ) (1 )ξ η ζ− − −

2Ψ ξ 18 (1 ) (1 ) (1 )ξ η ζ+ − −

3Ψ η 18 (1 ) (1 ) (1 )ξ η ζ+ + −

4Ψ ζ 18 (1 ) (1 ) (1 )ξ η ζ− + −

5Ψ 18 (1 ) (1 ) (1 )ξ η ζ+ − +

6Ψ 18 (1 ) (1 ) (1 )ξ η ζ+ − +

7Ψ 18 (1 ) (1 ) (1 )ξ η ζ+ + +

8Ψ 18 (1 ) (1 ) (1 )ξ η ζ− + +

Tabela 7: Elementos tetraédricos e hexaédricos

a) Elemento tetraédrico b) Elemento hexaédrico

Figura 69: Elementos isoparamétricos de 4 e 8 nós

Tal como para os elementos planos, a aproximação é, em geral, completa nos elementos

tetraédricos e incompleta nos elementos hexaédricos. Por exemplo, e de acordo com a Pirâmide

de Pascal, a aproximação é quadrática completa no elemento tetraédrico com 3 nós por aresta

(tetraedro de 10 nós), sendo cúbica incompleta (apenas contém o termo trilinear ξηζ ) no

elemento hexaédrico de 8 nós (a componente quadrática é também incompleta, pois não contém

os termos 2ξ , 2η e 2ζ , sendo apenas completa a componente linear, com a presença dos quatro

primeiros termos da Pirâmide de Pascal). Interessa também notar que as definições das funções

de aproximação nas faces dos elementos tetraédricos e hexaédricos contêm as funções de

aproximação dos elementos planos correspondentes.

ξ η

ζ

1

2 3

4

ξ

ζ

η 1

2 3

4

5

6 7

8


89

16.4 Aproximações

Os elementos isoparamétricos tridimensionais continuam a ser caracterizados por usarem as

mesmas funções de aproximação para mapear o elemento-mestre no elemento inserido na malha

de discretização,

1

1

1

( , , )

( , , )

( , , )

N

i ixi

N

i iyi

N

i izi

x c

y c

z c

Ψ ξ η ζ

Ψ ξ η ζ

Ψ ξ η ζ

=

=

=

=

=

=

∑

∑

∑

(163)

e para aproximar independentemente cada componente do campo de deslocamento:

1

1

1

( , , ) ( , , )

( , , ) ( , , )

( , , ) ( , , )

N

x i ixi

N

y i iyi

N

z i izi

u x y z d

u x y z d

u x y z d

Ψ ξ η ζ

Ψ ξ η ζ

Ψ ξ η ζ

=

=

=

=

=

=

∑

∑

∑

(164)

Os termos ( , , )ix iy izc c c na definição (163) caracterizam as coordenadas do nó i do elemento,

medidas no referencial global da malha. Na definição (164) os termos xu , yu e zu identificam as

componentes de deslocamento medidas no mesmo referencial e expressas em função das

componentes dos deslocamentos nodais, ( , , )ix iy izd d d .

É a seguinte a expressão geral da matriz das funções de aproximação na definição matricial

(35), de acordo com a definição (159) do vector dos deslocamentos:

1

1

1

0 0 0 0

0 0 0 0

0 0 0 0

N

N

N

Ψ ΨΨ Ψ

Ψ Ψ

=

ΨΨΨΨ⋯ ⋯ ⋯

⋯ ⋯ ⋯

⋯ ⋯ ⋯

A definição (37) mantém-se válida para a definição (36) das deformações compatíveis com a

aproximação do campo de deslocamentos, ficando,

1

1

1

1 1

1 1

1 1

0 0 0 0

0 0 0 0

0 0 0 0

0 0

0 0

0 0

x x N

y y N

z z N

y y N x x N

z z N x x N

z z N y y N

Ψ ΨΨ Ψ

Ψ ΨΨ Ψ Ψ ΨΨ Ψ Ψ Ψ

Ψ Ψ Ψ Ψ

∂ ∂ ∂ ∂ ∂ ∂

= ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂

∂ ∂ ∂ ∂

B

⋯ ⋯ ⋯

⋯ ⋯ ⋯

⋯ ⋯ ⋯

⋯ ⋯ ⋯

⋯ ⋯ ⋯

⋯ ⋯ ⋯

(165)


90

de acordo com a condição de compatibilidade (160) e a organização (156) do vector das

componentes de deformação. Também se mantém válida a definição (38) para a aproximação do

campo de tensão, usando agora a matriz de rigidez local (162).

Recorre-se à generalização da regra de diferenciação (137) para calcular os coeficientes da

matriz de compatibilidade (165),

( , , )

( , , )

( , , )

i i i ix i

i i i iy i

i i i iz i

x x x x

y y y y

z z z z

Ψ Ψ Ψ Ψξ η ζΨ ξ η ζξ η ζ



∂ ∂ ∂ ∂∂ ∂ ∂∂ = = + +∂ ∂ ∂ ∂ ∂ ∂ ∂

∂ ∂ ∂ ∂∂ ∂ ∂∂ = = + +∂ ∂ ∂ ∂ ∂ ∂ ∂

∂ ∂ ∂ ∂∂ ∂ ∂∂ = = + +∂ ∂ ∂ ∂ ∂ ∂ ∂

tomando a matriz Jacobiana da transformação (138) a seguinte expressão:

x y z

x y z

x y z

ξ ξ ξ

η η η

ζ ζ ζ

∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂= ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂

J

A matriz inversa é calculada numericamente, em cada ponto,

1

x x x

y y y

z z z

ξ η ζ

ξ η ζ

ξ η ζ

−

∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂= ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂

J

mantendo-se os cuidados a ter no mapeamento do elemento-mestre nos elementos da malha de

discretização, para garantir a estabilidade das operações decorrentes da transformação de

coordenadas.


A equação resolvente para um elemento tridimensional mantém a expressão geral (39),

continuando igualmente válidas as definições (40) para a matriz de rigidez e as definições (42) e

(43) para os vectores das forças nodais equivalentes às forças de massa e de fronteira. Também

continua válida a definição (47) para o vector das forças nodais equivalentes a campos de tensão

ou de deformação residuais, e a definição (48) para a contribuição de apoios elásticos para a

matriz de rigidez do elemento.


91

Todos estes integrais são calculados numericamente. Para elementos hexaédricos, os integrais

de domínio (146) e (148) são substituídos pelos seguintes,

1 1 1

1 1 1( , , ) ( , , ) ( , , )I f d f J d d dΩ Ωξ η ζ Ω ξ η ζ ξ η ζ ξ η ζ

+ + +

− − −= =∫ ∫ ∫ ∫

ou, numericamente:

1 1 1

( , , ) ( , , )M M M

i j k i j k i j ki j k

I f J w w wΩ ξ η ζ ξ η ζ= = =

=∑∑∑

É semelhante a generalização do cálculo de integrais de domínio para elementos tetraédricos.

Relativamente aos integrais de fronteira, são calculados nas faces desses elementos, tendo

portanto expressões análogas às dos integrais de domínio definidas para os elementos planos.

Após o cálculo dos termos da equação resolvente para cada elemento da malha de

discretização, aplica-se o processo de espalhamento anteriormente descrito para determinar os

termos correspondentes da equação resolvente (46) do modelo da estrutura em análise,

recorrendo a tabelas de incidência nodal análogas às anteriormente descritas para problemas

planos. A única diferença é que essa incidência é agora estabelecida para três componentes de

deslocamento e de força em cada nó da malha.

Em consequência da concepção em que se baseiam os elementos isoparamétricos, a condição

(44) de continuidade dos deslocamentos nodais assegura a continuidade de cada componente do

campo de deslocamentos nas interfaces dos elementos. Para além disso, o equilíbrio das

resultantes das forças nodais dos elementos que partilham um mesmo nó, definidas pela

condição (45), continua a caracterizar a função da equação resolvente (46).

Tal como para os problemas unidimensionais (estruturas articuladas e reticuladas) e para

problemas bidimensionais (estados planos de tensão e de deformação), a aproximação do campo

de deslocamentos em sólidos permite obter soluções aproximadas que são compatíveis mas que,

por regra, violam todas as condições de equilíbrio, tanto no domínio como na fronteira dos

elementos. A equação resolvente (46), ao estabelecer o equilíbrio das forças nodais equivalentes,

assegura que se obtém a melhor solução possível (em termos energéticos) com a aproximação

assumida para o campo de deslocamentos.

17. Utilização de Programas

As três grandes fases da análise de uma estrutura pelo Método dos Elementos Finitos foram

anteriormente caracterizadas e ilustradas na solução de exemplos de aplicação: a fase de pré-

processamento, em que se recolhem e armazenam os dados sobre a estrutura e o carregamento, e

se define a malha de discretização e o tipo de elemento; a fase de processamento, em que se


92

determina a equação resolvente de cada elemento da malha, se atribuem esses coeficientes ao

sistema resolvente da malha e se determina de solução do problema; e a fase de pós-

processamento, em que se caracteriza a resposta da estrutura ao carregamento dado, através da

representação da deformada e dos campos de tensão.

A execução dessas operações exige o recurso a técnicas específicas que visam a optimização

de cada fase do cálculo, para minimizar as exigências de armazenamento da informação, para

acelerar as operações e para assegurar a estabilidade numérica do processo de cálculo. Não são

essas técnicas, que quem deseje desenvolver ou adaptar um programa de elementos finitos deve

conhecer, que a seguir se descrevem.

O que se pretende é resumir a informação mais importante para apoiar a utilização de um

programa de elementos finitos, dando-se especial atenção às fases em que tem de existir, ou deve

existir, uma intervenção directa do utilizador, designadamente na definição dos dados e do

modelo de elementos finitos e, principalmente, de análise dos resultados.

Os programas comerciais partilham a mesma estrutura em termos de concepção mas podem

optar por diferentes estratégias na sequenciação de algumas operações. Todavia, o que mais os

distingue são os níveis de sofisticação que permitem na representação de geometrias e de acções,

na caracterização dos materiais estruturais e na selecção de modelos de análise. Incluem,

tipicamente, análises física e/ou geometricamente lineares, em regime estático ou dinâmico,

podendo ou não modelar a interacção entre a estrutura, a fundação e o carregamento.

Para além das potencialidades e das limitações do método e dos modelos de cálculo em que se

baseia, a boa utilização de um programa de elementos exige muita experimentação e o

conhecimento das consequências de cada uma das opções que o utilizador tem de tomar durante

o processo de análise. São essas opções que adiante se abordam, ainda que de uma maneira

muito limitada, de iniciação à utilização de programas de elementos finitos, pois uma informação

mais detalhada só pode ser apresentada no contexto de um dado programa e para um tipo

específico de análise estrutural.

17.1 Caracterização da Análise

Em consequência da generalidade conceptual do Método dos Elementos Finitos, que mais não é

do que um método de solução numérica de problemas às derivadas parciais, a maioria dos

programas oferece a possibilidade de resolver diferentes tipos de problemas definidos pelo

mesmo tipo de equações, designadamente da mecânica de estruturas, de fundações e de fluidos,

ou da térmica ou acústica, e, muitas vezes, o acoplamento desses problemas.


93

Por isso, a primeira informação a prestar é o tipo de problema (de análise estrutural) que se

pretende resolver. A segunda é sobre o modelo de análise a adoptar, em consequência das

simplificações que se podem alcançar atendendo à geometria das peças e às características do

carregamento. Em consequência da concepção do Método dos Elementos Finitos, esta segunda

informação é dada definindo (ainda em termos gerais) os graus de liberdade dos nós da malha,

isto é, as componentes do deslocamento e da rotação dos nós que podem ser não nulas.

Esta decisão do utilizador, ainda independente dos elementos a usar na análise, condiciona

todo o processo de cálculo, pelo que não pode ter dúvidas sobre os graus de liberdade de cada

tipo de peça estrutural que será utilizado e as acções a que estará sujeito. Um elemento

unidimensional poderá ter de dois a seis graus de liberdade por nó, as translações de um

elemento de treliça com comportamento plano e as três translações e as três rotações de um

elemento de pórtico com comportamento tridimensional, respectivamente. O mesmo sucede com

elementos planos, correspondendo o primeiro a um elemento para estados planos de tensão ou de

deformação e o segundo a um elemento de casca.

É também nesta fase de entrada de dados que se define o tipo de análise estrutural,

designadamente se é uma análise estática ou dinâmica, e se é geometricamente linear

(deformações e deslocamentos infinitesimais) ou não linear (grandes deformações e/ou

deslocamentos).

Dependendo das potencialidades e da flexibilidade do programa, pode ser necessário definir

ainda nesta fase outros aspectos da análise, podendo uns ser do âmbito da concepção do modelo

estrutural e outros do âmbito do processamento numérico. São exemplos a modelação de tensões

ou de deformações residuais, os métodos de armazenamento da matriz do sistema resolvente e de

solução desse sistema, assim como opções específicas sobre execução de determinados tipos de

análise, designadamente análises dinâmicas ou que envolvam a modelação da interacção da

estrutura com a fundação ou com a acção, por exemplo em problemas de acústica estrutural e de

aerodinâmica estrutural.

17.2 Caracterização da Geometria

Apesar da caracterização da geometria estar muito facilitada nos programas comerciais, esta fase

continua a ser a que consome mais tempo na intervenção do utilizador. A estratégia mais

recomendada é a de decompor o domínio em subdomínios com geometrias simples e, se

possível, que ocorram com alguma frequência, tanto no problema que se pretende analisar como

nos tipos de problemas frequentemente analisados.


94

É importante conceber cuidadosamente a definição de cada subdomínio, pois uma

sistematização adequada facilita substancialmente as fases seguintes de caracterização do

material estrutural e das condições de fronteira do problema. Para além disso, permite criar

ficheiros de dados que podem ser usados em análises posteriores.

A geometria desses subdomínios é caracterizada definindo as coordenadas dos nós

estritamente necessários para definir os elementos que a caracterizam, designadamente

segmentos, superfícies e volumes, consoante o tipo de aplicação. É em geral possível

sistematizar a caracterização de elementos geométricos através de expressões analíticas simples.

Esses dados são definidos separadamente, para apoiar a reunião dos subdomínios,

designadamente os segmentos ou as superfícies comuns em geometrias planas ou

tridimensionais.

17.3 Caracterização do Material Estrutural

Os programas permitem a utilização de uma grande variedade de modelos para as relações

constitutivas, designadamente modelos elásticos, elastoplásticos e viscoelásticos. Podem incluir

diferentes tipos de anisotropia e de não linearidade, dependendo do nível de deformação e/ou da

evolução das propriedades mecânicas no tempo ou, por exemplo, com a temperatura.

Alguns programas permitem que seja o utilizador a criar (ou seja, a programar) um

determinado modelo constitutivo. No entanto, a tendência é para cada programa se especializar

num determinado tipo de material estrutural, tipicamente para estruturas metálicas ou para

estruturas de betão armado pré-esforçado, caso em que disponibilizam os modelos definidos nos

regulamentos que lhes são específicos e nas condições aí definidas.

17.4 Caracterização das Condições de Fronteira

É na caracterização do carregamento e, principalmente, nas condições de apoio que ocorrem com

mais frequência os erros de utilização de programas de elementos finitos.

Para além da representação gráfica desses dados nem sempre ser suficientemente clara, é fácil

confundir o referencial em que os dados são definidos, no referencial global da estrutura ou nos

referenciais locais dos elementos geométricos a que as cargas ou as condições de apoio serão

atribuídos.

A definição das forças aplicadas e das condições de apoio também devem ser tipificadas por

grupos, para depois facilitar a sua atribuição aos elementos usados na caracterização da

geometria, designadamente, pontos, segmentos e superfícies.

Para além de forças (momentos) e de deslocamentos (rotações) impostos em pontos, os

programas definem um conjunto de opções fixas que traduzem os carregamentos e as condições


95

de apoio usadas com mais frequência. Está também generalizada a oferta da opção de fixar a

posição relativa de pontos da estrutura ligando-os por bielas (ou planos) rígidos.

No entanto, é muito variável a flexibilidade que os programas oferecem para a caracterização

de forças e de deslocamentos variáveis no espaço (e no tempo, em análises dinâmicas), o mesmo

sucedendo em relação à modelação de apoios deformáveis (por exemplo apoios elásticos) ou

unidireccionais (que apenas resistem a forças de contacto), com ou sem atrito.

Tal como para a caracterização dos materiais estruturais, os programas que se especializam

em determinado tipo de solução estrutural incluem os modelos regulamentares que lhes são

específicos, tanto em termos de acções como de condições de apoio.

17.5 Caracterização da Malha de Discretização

A definição de subdomínios com geometrias simples facilita a geração de malhas regulares, com

uma transição progressiva na forma e nas dimensões das células, que posteriormente serão

identificadas com os elementos finitos. Para além disso, o cuidado posto na definição dos

elementos geométricos de interface (os segmentos ou superfícies que são partilhados por

subdomínios) facilita a regularização e a compatibilização das malhas definidas para cada

subdomínio.

Os programas oferecem diferentes opções para gerar as malhas de discretização dos diferentes

subdomínios definidos na análise. Todas elas visam assegurar um determinado nível de

densidade da discretização, podendo escolher-se diferentes critérios de segmentação (tipicamente

leis de progressão aritmética ou geométrica) e escolher dimensões mínimas e máximas. É

também possível escolher o tipo de célula, triangular (tetraédrica) ou quadrangular (hexaédrica)

e, com alguma frequência, o algoritmo de geração da malha.

É essencial analisar criticamente as malhas obtidas antes de avançar para a escolha do tipo de

elemento a utilizar no cálculo. Devem ser rejeitadas as malhas que apresentem distorções fortes

ou incompatibilidades, e deve-se assegurar que é suficiente o refinamento nas zonas do modelo

em que se esperam grandes concentrações de tensões ou esforços.

Só depois se deve escolher o tipo de elemento finito a adoptar na análise. Tem de ser do tipo

especificamente desenvolvido para o modelo estrutural em análise, deve ser coerente com as

células da malha de discretização e deve ter um grau de aproximação adequado à densidade da

malha que foi criada.

Por coerência com as células da malha entende-se que devem ser usados elementos com a

mesma topologia, triangulares ou quadrangulares. Quando se usa este último tipo de elemento

numa malha com células triangulares, ou se opta pelo inverso, está-se implicitamente a introduzir


96

uma distorção sistemática dos elementos finitos. O mesmo se aplica, naturalmente, a malhas

tridimensionais definidas por células tetraédricas/hexaédricas, as quais devem ser implementadas

com elementos do mesmo tipo.

É óbvia a necessidade de usar o tipo de elemento finito adequado ao modelo estrutural da

peça, designadamente elementos triangulares ou quadrangulares com dois graus por nó para

estados planos de tensão ou de deformação e com três graus de liberdade por nó para estados

tridimensionais. É também intuitivo que não devem ser usadas malhas relativamente grosseiras

com elementos com poucos nós ou que pode ser excessivo usar malhas muito densas com

elementos com muitos nós.

No entanto, e antes de se começar a conhecer razoavelmente um dado programa, é frequente a

dúvida sobre os elementos a utilizar na modelação, especialmente na análise de lajes e cascas,

finas ou espessas em ambos os casos. Há uma gama de elementos disponíveis bastante vasta,

podendo serem distintos os tipos de elementos disponibilizados por diferentes programas. É

indispensável consultar o manual do programa e seguir as recomendações aí dadas sobre a

adequabilidade de cada tipo de elemento para os diferentes tipos de análise.

Concluída a caracterização da malha de discretização de cada subdomínio e a escolha dos

elementos finitos a adoptar, é importante representar a malha e activar a opção que permite

sobrepor os elementos utilizados para verificar a compatibilidade da malha em termos de

incidência nodal. Apesar da solução obtida com uma malha incompatível não ser

necessariamente pior que a obtida com uma malha compatível, deixam de estar garantidas as

condições de convergência da formulação, as quais é importante garantir para poder ajuizar,

posteriormente, se será necessário refinar essa solução.

Para potenciar a sistematização da definição de geometrias e malhas, os programas permitem

importar e combinar ficheiros em que esses dados estão caracterizados, podendo pôr-se o

problema de compatibilizar malhas e elementos. Como se referiu anteriormente, os programas

podem disponibilizar elementos de transição para facilitar essa compatibilização.

17.6 Análise da Solução

Concluída a caracterização do modelo de elementos finitos, processa-se a montagem do sistema

resolvente e o cálculo da solução (ou das soluções em cada incremento em análises dinâmicas

e/ou não lineares). Essa informação é registada em ficheiros, assim como a informação adicional

que o utilizador possa solicitar, tipicamente a energia de deformação, os deslocamentos em

determinados nós e as tensões ou as deformações em determinados pontos (tipicamente os

pontos de Gauss, onde o erro é minimizado) de determinados elementos.


97

É indispensável analisar cuidadosamente os resultados obtidos para responder a duas questões

essenciais: se o problema resolvido corresponde ao que se pretendia analisar e, em caso

afirmativo, se o nível de precisão obtido é adequado para os fins dessa análise.

Essa decisão é fundamentada combinando dois conjuntos de informação: as condições de

fronteira do problema, pois o carregamento e as condições de apoio são dados do problema; as

condições de aproximação em que o método de cálculo se baseia, pois sabe-se quais as

condições que devem ser satisfeitas localmente e quais as que são impostas de maneira

aproximada. Os comentários seguintes pressupõem uma análise elástica linear, pois podem ser

diferentes as condições impostas de maneira exacta em problemas baseados noutras hipóteses.

No contexto da análise elástica linear, a solução obtida tem de ser estritamente compatível. É

desnecessário verificar se as deformações são compatíveis com a aproximação dos

deslocamentos pois essa condição é implícita à formulação do método de cálculo. Todavia é

necessário verificar se o campo de deslocamentos é contínuo e se satisfaz as condições de apoio

do problema dado.

Essa verificação faz-se facilmente analisando a deformada da estrutura. Terá havido erro do

utilizador na definição da malha de discretização ou na escolha dos elementos finitos se se

verificarem vazios ou sobreposições entre elementos da malha. E terá havido erro do utilizador

na definição das condições de apoio se estiverem livres movimentos que deveriam estar

restringidos, ou se se verificar o inverso.

É também desnecessário verificar se as relações de elasticidade são impostas localmente, pois

está também implícito na formulação que as tensões são calculadas a partir da aproximação do

campo de deformações impondo essas relações. O que deve ser cuidadosamente verificado é o

erro nas condições de equilíbrio no domínio e na fronteira dos elementos, as quais se sabe serem

impostas apenas de maneira aproximada, através do equilíbrio das forças nodais equivalentes. A

facilidade com que essas condições podem ser avaliadas depende muito da sofisticação do

programa de elementos finitos.

Tal como se fez nos exemplos de aplicação, é possível determinar e representar graficamente

o campo das forças de massa que equilibram as tensões calculadas e compará-las com as forças

de massa dadas para avaliar a qualidade do equilíbrio. Também como se fez nesses exemplos, é

possível calcular e representar graficamente as forças de fronteira que equilibram esse campo de

tensões e verificar o desequilíbrio de forças entre elementos e em relação às forças de fronteira

definidas no carregamento (as condições de fronteira estáticas do problema).


98

No entanto, o que todos os programas disponibilizam é suficiente para avaliar o equilíbrio da

solução e os eventuais erros na definição do carregamento. Todos os programas disponibilizam

representações gráficas dos campos de tensão, estando sempre presentes duas opções: a

representação das tensões calculadas e a representação da regularização dessa solução (sendo

esta opção designada por smoothing, averaging ou regularization). Os programas disponibilizam

também, em modo numérico ou gráfico, as reacções de apoio (de facto, as forças nodais

equivalentes correspondentes aos deslocamentos restringidos), as quais podem ser usadas para

avaliar se o carregamento foi correctamente definido.

A representação do campo de tensões não regularizado serve sempre para decidir sobre a

inadequabilidade da solução. O nível de convergência é necessariamente insatisfatório se as

descontinuidades no campo de tensão são tais que reflectem a malha de elementos finitos usada

no cálculo. Esse resultado deve servir, no entanto, para decidir como melhorar a solução,

refinando as zonas onde os maiores desequilíbrios se verificam. Pode também servir para expor

zonas de concentração de tensões que não haviam sido consideradas na fase de definição da

malha de discretização.

Devem ser feitas duas verificações quando se obtêm campos praticamente contínuos na

representação não regularizada das componentes de tensão. A primeira é confirmar se as tensões

equilibram as forças de fronteira do problema. Se tal não suceder, terá havido erro na definição

do carregamento. A segunda verificação a fazer é sobre as concentrações de tensão reveladas

nessa solução. Se essas concentrações não podem ser justificadas pelo carregamento (cargas

pontuais ou descontínuas), pelas descontinuidades das condições de apoio ou da geometria

(tipicamente cantos reentrantes), ou por restrições impostas sobre o movimento relativo de nós

da malha (tipicamente bielas rígidas), é muito provável que tenha havido erro nas restrições

postas sobre os graus de liberdade dos nós da malha.

É aceitável uma solução que satisfaça os critérios de análise acima resumidos. No entanto, e

como a maioria dos programas não disponibiliza estimadores de erro, se a aplicação exige um

nível de segurança mais apertado, a única opção possível é realizar estudos de convergência, os

quais podem ser feitos sobre a mesma malha e melhorando a qualidade dos elementos, refinando

local ou globalmente a malha ou mesmo gerando uma malha diferente. Para cada um desses

casos traça-se a variação da energia de deformação (disponibilizada pela maioria dos

programas), de deslocamentos nodais representativos e de tensões em elementos mais sensíveis.

Ao analisar esses resultados é necessário notar que a energia de deformação tende a convergir

mais rapidamente que os deslocamentos e estes mais rapidamente que as tensões.


99

18. Conclusão

Baseando-se na mesma formulação do Método dos Elementos Finitos, os conceitos fundamentais

que a seguir se resumem repõem os já identificados na análise de estruturas articuladas e

reticuladas. De facto, o único conceito fundamentalmente diferente que foi introduzido nesta

generalização do método para a solução de problemas planos e tridimensionais foi o conceito de

elemento isoparamétrico, o qual também é utilizado no contexto das peças lineares:

• Definem-se elementos com geometria simples num espaço independente daquele em que

se realiza a análise da estrutura, os elementos-mestre triangulares/tetraédricos ou

quadrangulares/hexaédricos, definidos no espaço das coordenadas naturais;

• A estrutura é decomposta em elementos finitos isoparamétricos, os quais se caracterizam

por usar o mesmo conjunto de funções para mapear o elemento-mestre nos elementos da

malha de discretização da estrutura e para aproximar a variação do campo de

deslocamentos nesses elementos;

• Essas funções, cujo grau de variação é determinado pelo número de nós do elemento, são

contínuas e definidas de maneira a facilitar a imposição das condições de fronteira

cinemáticas e de continuidade dos deslocamentos entre elementos;

• A aproximação das componentes de deformação é calculada impondo localmente a

condição de compatibilidade no domínio do elemento e a aproximação das componentes de

tensão é calculada impondo localmente a condição de elasticidade, transferindo-se o erro

da aproximação para as condições de equilíbrio no domínio e na fronteira do elemento;

• Essas condições são impostas exigindo apenas o equilíbrio dos nós do elemento,

recorrendo-se à definição de forças nodais estaticamente equivalentes para obter a equação

resolvente elementar, cujos coeficientes são calculados por integração numérica no espaço

do elemento-mestre;

• As equações resolventes elementares são combinadas para obter a equação resolvente do

modelo de elementos finitos, sendo essa combinação determinada por critérios de

incidência nodal que identificam os deslocamentos nodais dos elementos com os

deslocamentos nodais da estrutura e definem, implicitamente, as forças nodais da estrutura

como as resultantes das forças nodais dos elementos;

• A solução que se obtém é (necessariamente) cinematicamente admissível e recupera

(necessariamente) a solução exacta do problema sempre que a aproximação em cada

elemento permita obter uma solução que é aí estaticamente admissível;


100

• Quando tal não sucede, as condições de equilíbrio no domínio e na fronteira desses

elementos são violadas localmente mas impostas de maneira fraca, através de uma

condição de equilíbrio imposta sobre as forças nodais equivalentes, as quais satisfazem

(necessariamente) as condições de equilíbrio global da estrutura, mesmo que a solução não

seja exacta;

• Uma solução aproximada pode ser sempre melhorada, convergindo para a solução exacta

(a qual geralmente não tem expressão analítica) quando se subdivide os elementos

(refinamento-h), se aumenta o grau da aproximação em cada elemento (refinamento-p), ou

se combinam os dois processos.

Documents

Introdução ao Método dos Elementos Finitos: Elasticidade ...luis/ae2/AE2_Elasticidade2D_v4.pdf · Análise de Estruturas II: Elasticidade Plana e Tridimensional 3 2. Hipóteses