APOSTILA - MEF - 2011 - Rade

Prévia do material em texto

UNIVERSIDADE FEDERAL DE UBERLÂNDIA 
 FACULDADE DE ENGENHARIA MECÂNICA 
Programa de Pós-Graduação em Engenharia Mecânica 
INCT de Estruturas Inteligentes em Engenharia 
Laboratório de Mecânica de Estruturas Prof. José Eduardo Tannús Reis 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Prof. Domingos Alves Rade 
 
 
 
 
 
 
 
2011 
 
MÉTODO DOS ELEMENTOS FINITOS
APLICADOS À ENGENHARIA 
MECÂNICA 
Introdução ao Método dos Elementos Finitos D.A. Rade 
 1
CAPÍTULO 1 
 
INTRODUÇÃO AO MÉTODO DOS ELEMENTOS FINITOS 
 
 
1.1 - Fundamentos do método dos elementos finitos 
 
O método dos elementos finitos (MEF) é uma técnica de análise numérica 
destinada à obtenção de soluções aproximadas de problemas regidos por equações 
diferenciais. Embora o método tenha sido originalmente desenvolvido para a análise 
estática de sistemas estruturais, ele tem sido utilizado no estudo de uma grande 
variedade de problemas de Engenharia, nos domínios da Mecânica dos Sólidos, 
Mecânica dos Fluidos, Transmissão de Calor e Eletromagnetismo. Devido à sua 
eficiência e flexibilidade, além de sua adequação à implementação em computadores 
digitais, o MEF tem hoje uma grande difusão tanto no meio acadêmico como no 
industrial, estando disponível em grande número de “pacotes” comerciais existentes 
no mercado (ANSYS®, NASTRAN®, ABAQUS®, SYSTUS®, COMSOL®, etc.). 
Contudo, deve ser lembrado que a utilização eficaz destes programas e a correta 
interpretação dos resultados requerem o amplo conhecimento, por parte do 
Engenheiro, dos fundamentos do MEF. 
 A principal motivação para o uso do MEF reside no fato que, devido à 
complexidade dos problemas práticos de Engenharia, soluções analíticas em forma 
fechada tornam-se inviáveis ou mesmo impossíveis. Assim, devemos recorrer a 
técnicas capazes de fornecer soluções numéricas aproximadas. A título de exemplo, 
consideremos os problemas de determinação da capacidade de carga de uma placa 
contendo enrijecedores e entalhes de formas complexas, ou de determinação da 
concentração de poluentes sob condições atmosféricas não uniformes, ou ainda de 
caracterização do perfil de velocidades em torno de um aerofólio. Para cada um 
destes problemas podemos obter, sem grande esforço, as equações governantes e as 
condições de contorno, utilizando princípios elementares da Física. Contudo, 
nenhuma solução analítica simples poderá ser obtida quando o problema exibir 
geometria e/ou condições de contorno complicadas, o que quase sempre ocorre em 
situações práticas. Para contornar esta dificuldade, uma estratégia possível é a 
simplificação do problema (em termos de sua geometria e/ou condições de contorno) 
de modo a viabilizar a construção de um modelo matemático cuja resolução analítica 
seja possível. Contudo, em grande número de casos (talvez na maioria das vezes), 
este procedimento tem como conseqüência graves imprecisões nas previsões do 
modelo. Uma segunda alternativa consiste em preservar a complexidade do modelo 
e empregar técnicas aproximadas de resolução. Esta segunda estratégia, na qual 
está inserido o MEF, tem sido cada vez mais viabilizada pela crescente capacidade 
de processamento dos computadores digitais. 
 Em todo problema formulado em domínios contínuos, as incógnitas do 
problema, denominadas variáveis de campo (que podem ser grandezas escalares, 
como temperaturas ou vetoriais, como deslocamentos) podem assumir valores 
independentes em cada ponto do domínio. Conseqüentemente, o problema tem 
número infinito de incógnitas, sendo caracterizado como um problema 
infinito-dimensional. Este tipo de problema é geralmente modelado por equações 
diferenciais parciais, cuja solução analítica é dada por funções que fornecem os 
Introdução ao Método dos Elementos Finitos D.A. Rade 
 2
valores das variáveis de campo em função das coordenadas espaciais para todos os 
pontos do domínio. 
O MEF é essencialmente um processo de discretização, que visa transformar 
um problema infinito-dimensional em um problema finito-dimensional, com número 
finito de incógnitas. O método consiste em dividir o domínio sobre o qual o problema 
é estudado em várias regiões interconectadas, denominadas elementos. Cada 
elemento dispõe de um certo número de pontos (interiores e/ou limítrofes), 
denominados nós ou pontos nodais. O conjunto de elementos utilizados na 
discretização é denominado malha. Um exemplo é apresentado na Figura 1.1, que 
mostra a seção transversal de uma palheta de turbina de geometria complexa, 
discretizada em elementos de forma triangular, tendo, em cada vértice, um nó. 
Neste exemplo, o problema em questão poderia ser a determinação da distribuição 
de temperaturas sobre a seção da palheta, conhecidos o fluxo de calor e as condições 
de contorno. 
Uma vez definidos os elementos e seus respectivos nós, no interior de cada 
elemento são admitidas soluções aproximadas para as variáveis de campo, 
expressas como funções arbitrárias dos valores que as incógnitas assumem nos nós 
(valores nodais). Estas funções são denominadas funções de interpolação ou funções 
de forma. São também impostas condições garantindo a continuidade da solução no 
nós compartilhados por vários elementos. As incógnitas do problema, denominadas 
graus de liberdade (g.d.l.), passam a ser os valores das variáveis de campo nos 
pontos nodais, sendo o número destas incógnitas (agora finito), denominado número 
de graus de liberdade do modelo. Dependendo da natureza do problema, após a 
discretização, o modelo matemático regente resulta representado por um número 
finito de equações diferenciais ordinárias ou de equações algébricas, cuja resolução 
numérica conduz aos valores das incógnitas nodais. Uma vez determinadas estas 
incógnitas, os valores das variáveis de campo no interior dos elementos podem ser 
avaliados empregando as funções de interpolação. 
Conforme será visto mais adiante, a precisão da solução obtida depende 
essencialmente do número de elementos e do tipo de funções de forma empregadas 
na discretização. Sendo satisfeitas algumas condições, admite-se que a solução do 
problema discretizado convirja para a solução exata do problema contínuo à medida 
que se aumenta o número de incógnitas nodais. 
 
 
 
Figura 1.1 – Ilustração da malha de um modelo de elementos finitos 
Introdução ao Método dos Elementos Finitos D.A. Rade 
 3
Em comparação com outras técnicas numéricas, as principais vantagens o 
método dos elementos finitos são as seguintes: 
 
• elementos de diferentes formas e tamanhos podem ser associados para 
discretizar domínios de geometria complexa. 
• a divisão do contínuo em regiões facilita a modelagem de problemas 
envolvendo domínios não homogêneos, onde as propriedades físicas variam em 
função das coordenadas espaciais. 
• o método pode ser todo formulado matricialmente, facilitando sua 
implementação computacional. 
 
 A implementação do MEF pode sempre ser efetuada em etapas sucessivas, de 
forma estruturada. As principais etapas são as seguintes: 
 
1ª) Discretização do domínio. O primeiro passo é a divisão do domínio em 
elementos. O tipo e número de elementos a serem utilizados devem ser escolhidos de 
modo a representar adequadamente a geometria do problema e caracterizar 
convenientemente as variações da solução ao longo do domínio. Alguns tipos de 
elementos freqüentemente empregados para a discretização de domínios 
unidimensionais, bidimensionais e tridimensionais são ilustrados na Figura 1.2. 
Neste aspecto, deve-se observar que problemas unidimensionais são aqueles 
definidos em domínios representados por apenas uma coordenada espacial (linhas),ao passo que problemas bidimensionais e tridimensionais são aqueles definidos em 
domínios representados por duas coordenadas espaciais (superfícies) e três 
coordenadas espaciais (volumes), respectivamente. Os elementos axissimétricos, 
mostrados na Figura 1.2, são elementos utilizados para a discretização de 
problemas tridimensionais caracterizados pela existência de simetria geométrica e 
de carregamento em relação a um dado eixo. Neste caso, o problema tridimensional 
pode ser formulado como um problema bidimensional. 
 
 
 
(a) 
 
elementos unidimensionais 
 
 
(b) 
 
elementos bidimensionais 
Introdução ao Método dos Elementos Finitos D.A. Rade 
 4
 
 
(c) 
 
elementos tridimensionais 
 
(d) 
 
elemento axissimétrico 
 
Figura 1.2 – Ilustração de diferentes tipos de elementos 
 
2ª) Escolha das funções de interpolação. Nesta etapa são escolhidas as funções 
de interpolação que representam as variáveis de campo no interior de cada 
elemento. Freqüentemente, mas nem sempre, funções polinomiais são escolhidas 
como funções de interpolação, devido à facilidade que oferecem para derivação e 
integração. Os graus dos polinômios utilizados estão relacionados ao número de 
incógnitas nodais de cada elemento, devendo também atender a certos requisitos de 
continuidade das variáveis de campo a serem satisfeitos nos nós e nas fronteiras 
entre elementos imediatamente vizinhos. 
 
3ª) Construção das matrizes elementares. Uma vez escolhidos o tipo e número 
de elementos e as funções de interpolação, devemos estabelecer as relações 
matriciais expressando o comportamento (relações de causa-efeito), em termos de 
propriedades físicas e geométricas, para cada elemento, individualmente. Em outras 
palavras, procede-se à formulação em nível elementar. Para tanto, podem ser 
utilizados os seguintes processos: 
 
• processo direto, que é baseado no método da rigidez da análise estrutural, através 
do qual são obtidas as relações matriciais entre forças e deslocamentos nodais, a 
partir das relações de equilíbrio de forças e compatibilidade de deslocamentos. 
Procedimento similar pode ser utilizado na modelagem de problemas 
unidimensionais de transmissão de calor. Embora o processo direto só seja 
conveniente no tratamento de problemas mais simples, ele tem a grande vantagem 
de ser de fácil entendimento, permitindo clara interpretação física do significado das 
relações matriciais obtidas, interpretação esta que pode ser estendida a problemas 
mais complexos. 
 
• processo variacional, que é baseado no Cálculo Variacional e envolve a busca dos 
pontos críticos – geralmente pontos de mínimo – de um funcional associado ao 
problema estudado. De acordo com este processo, as relações matriciais em nível 
Introdução ao Método dos Elementos Finitos D.A. Rade 
 5
elementar resultam da imposição da condição de estacionaridade do funcional 
associado ao problema. Em problemas de Mecânica dos Sólidos, por exemplo, este 
funcional pode representar a energia de deformação ou a energia potencial 
complementar. O processo variacional, embora mais complicado que o processo 
direto sob o ponto de vista teórico, permite estender o MEF a problemas mais 
complexos. Contudo, sua aplicabilidade é limitada aos problemas regidos por 
princípios variacionais, que estabelecem a existência de funcionais. 
 
• processo dos resíduos ponderados. Este é um procedimento ainda mais versátil 
que os dois procedimentos anteriores, sendo baseado integralmente em operações 
matemáticas. O processo opera diretamente sobre as equações diferenciais que 
governam o problema e prescinde da existência de um funcional ou de um princípio 
variacional. 
 
4ª) Montagem das matrizes elementares para obtenção das matrizes 
globais. Para caracterizar o comportamento do sistema completo, resultante da 
associação dos vários elementos, devemos agrupar as matrizes de cada um dos 
elementos de uma forma adequada. Em outras palavras, devemos combinar as 
equações matriciais expressando o comportamento dos elementos individuais para 
formar as equações matriciais que descrevem o comportamento do sistema em todo 
o domínio. Este processo é conhecido como montagem das matrizes globais. No 
processo de montagem, impõe-se a condição que em cada nó onde vários elementos 
estão interconectados, os valores das variáveis de campo são os mesmos para cada 
elemento compartilhando aquele nó. 
 
No final deste processo, as equações matriciais globais devem ser modificadas para 
satisfazer as condições de contorno do problema. A ordem das matrizes globais 
coincide com o número total de incógnitas nodais. Este número é chamado número 
de graus de liberdade do modelo. 
 
5ª) Imposição dos carregamentos externos e das condições de contorno. As 
equações matriciais globais devem ser modificadas para satisfazer as condições de 
contorno do problema, que expressam o fato que alguns valores das incógnitas 
nodais são prescritos. Assim, por exemplo, em problemas de transferência de calor, 
os valores da temperatura em alguns pontos do contorno podem ser previamente 
conhecidos. Da mesma forma, deve-se alterar as equações globais para leva em 
conta que, em alguns nós, cargas externas conhecidas (forças, fluxos de calor, etc.) 
são aplicadas. Ao final deste processo, o número total de incógnitas nodais 
remanescentes define o chamado número de graus de liberdade do modelo. 
 
6ª) Resolução do sistema de equações. Ao final do processo de montagem das 
matrizes globais, o modelo matemático do problema estará representado por um 
conjunto de equações, que podem ser lineares ou não lineares, algébricas ou 
diferenciais, dependendo da natureza do problema enfocado. Estas equações devem 
então ser resolvidas numericamente para a determinação dos valores das variáveis 
de campo nos pontos nodais. Neste processo de resolução, procedimentos numéricos 
apropriados, implementados sob a forma de rotinas computacionais, devem ser 
utilizados. 
 
Introdução ao Método dos Elementos Finitos D.A. Rade 
 6
7ª) Realização de cálculos complementares. Em várias situações, cálculos 
complementares devem ser realizados para a determinação de grandezas 
dependentes das variáveis de campo, determinadas na etapa precedente. Assim, por 
exemplo, nos problemas de Mecânica dos Sólidos, uma vez determinados os 
deslocamentos, cálculos adicionais são necessários para a determinação das 
deformações (utilizando as relações deformação-deslocamento) e das tensões 
(utilizando as relações tensão-deformação). 
 
 
1.2 - Domínios de aplicação do MEF 
 
 Os problemas passíveis de tratamento pelo MEF podem ser divididos em três 
categorias: 
 
• problemas de equilíbrio. Esta é a classe de problemas cuja solução é independente 
do tempo. São exemplos os problemas da Mecânica dos Sólidos envolvendo a 
determinação de tensões e deformações em elementos estruturais submetidos a 
carregamentos estáticos e os problemas da Mecânica dos Fluidos tratando da 
determinação de distribuições de pressão, velocidade em regime permanente e os 
problemas de Transferência de Calor em regime permanente. Para este tipo de 
problema, o processo de discretização através do MEF conduz a um modelo 
matemático representado por um conjunto de equações algébricas, que podem ser 
lineares ou não lineares. 
 
• problemas de autovalor. Nesta classe de problemas, o modelo matemático obtido é 
representado por um conjunto de equações lineares homogêneas, caracterizado pela 
dependência em relação a um parâmetro, cuja resolução conduz a um conjunto de 
autovalores e autovetores. São exemplos os problemas que tratam da determinação 
de freqüências naturaise modos de vibração de meios sólidos e fluidos, além de 
cargas de flambagem de elementos estruturais. No primeiro caso os autovalores 
correspondem às freqüências naturais e os autovetores associam-se aos modos 
naturais de vibração; no segundo, os autovalores correspondem às cargas de 
flambagem e os autovetores dizem respeito aos campos de deslocamentos 
correspondentes. 
 
• problemas de propagação. Os problemas de propagação são aqueles em que se 
busca caracterizar a evolução das variáveis de campo em função do tempo. É o caso 
típico de fenômenos que se desenvolvem em regime transitório. Os seguintes 
exemplos podem ser mencionados: determinação do movimento de sistemas 
estruturais submetidos a cargas de impacto e determinação de distribuições de 
temperatura geradas por fluxos de calor variáveis. 
 
 Alguns exemplos de aplicações práticas do MEF são ilustrados a seguir. 
 
 
 
Introdução ao Método dos Elementos Finitos D.A. Rade 
 7
 
 
 
 
Figura 1.3 - Modelo de EF para análise aeroelástica de um Airbus A-320 (100.000 g.d.l.) 
(extraído de J.F. Imbert, "Analyse des Structures par Eléments Finis) 
 
 
 
 
 
 
 
 
 
Introdução ao Método dos Elementos Finitos D.A. Rade 
 8
 
(b) 
 
Figura 1.4 - Análise térmica por EF de uma tubulação em forma de estrela 
(de Desai C.S. e Abel J.F., "Introduction to the Finite Element Method") 
 
 
 
Introdução ao Método dos Elementos Finitos D.A. Rade 
 9
 
 
 
(c) 
 
Figura 1.5 - Análise por EF do escoamento de um fluido ideal em torno de um cilindro 
posicionado entre duas placas paralelas 
(de Desai C.S. e Abel J.F., "Introduction to the Finite Element Method") 
 
 
 
 
 
 
Introdução ao Método dos Elementos Finitos D.A. Rade 
 10
 
 
 
Figura 1.6 - Análise por EF do resfriamento de um lingote metálico 
(fonte: http://www.comsol.com (Comsol Group)) 
 
 
 
 
Figura 1.7 - Análise dinâmico de freios automotivos pelo MEF 
(fonte: http://www.simulia.com (Dassault Systems)) 
 
 
 
 
 
 
 
 
 
 
 
 
Introdução ao Método dos Elementos Finitos D.A. Rade 
 11
 
Figura 1.8 – Simulação numérica do processo de soldagem pelo MEF 
(fonte: http://www.simulia.com (Dassault Systems)) 
 
 
 
 
 
Figura 1.9 – Simulação do comportamento da bexiga humana pelo MEF 
(fonte: http://www.comsol.com (Comsol Group)) 
 
 
 
Introdução ao Método dos Elementos Finitos D.A. Rade 
 12
 
Figura 1.10 – Caracterização do fluxo magnético em um gerador pelo MEF 
(fonte: http://www.comsol.com (Comsol Group)) 
 
 
 
 
 
Figura 1.11 – Simulação de impacto em palhetas de uma turbina pelo MEF 
(fonte: http://www.ansys.com (Ansys, Inc.)) 
 
 
 
 
 
 
 
Introdução ao Método dos Elementos Finitos D.A. Rade 
 13
 
 
 
Figura 1.12 – Simulação de “crash-test” pelo MEF 
(fonte: http://www.ansys.com (Ansys, Inc.)) 
 
 
1.3 - Sobre as limitações do método dos elementos finitos 
 
 Embora o MEF seja reconhecidamente uma ferramenta útil e eficiente para a 
análise de diversos tipos de problemas de Engenharia, é importante estar atento às 
limitações do método, lembrando que ele fornece modelos matemáticos aproximados 
para representar o comportamento de sistemas físicos. Conforme será evidenciado 
no desenvolvimento do curso, em geral, a elaboração de um modelo de EF envolve a 
admissão de uma série de simplificações que terão efeito direto sobre a precisão das 
previsões do modelo. Deve também ser levado em conta que, na maioria das vezes, 
sob considerações de custo e limitações de recursos computacionais, busca-se 
estabelecer um compromisso entre a complexidade do modelo e a precisão 
considerada satisfatória. 
Algumas fontes de incerteza inerentes à modelagem por EF são: 
 
 • a não consideração de certos tipos de efeitos físicos, tais como não 
linearidades, histerese, amortecimento, etc. 
 • erros de discretização, devidos à impossibilidade de se obter uma perfeita 
representação de domínios de geometria complexa utilizando os tipos de elementos 
disponíveis. 
 • conhecimento impreciso dos valores de alguns parâmetros físicos e/ou 
geométricos que são utilizados na elaboração do modelo (ex.: módulo de elasticidade, 
densidade, condutividade térmica, viscosidade, etc.) 
 • dificuldade de modelar efeitos localizados, tais como junções parafusadas e 
rebitadas. 
 • erros oriundos do processo de resolução numérica. 
 
Introdução ao Método dos Elementos Finitos D.A. Rade 
 14
Visando a validação de modelos de elementos finitos, um procedimento 
recomendável é a confrontação de suas previsões com dados provenientes de outros 
tipos de análises, em particular, de análises experimentais. Várias técnicas foram 
desenvolvidas, principalmente no âmbito da Dinâmica Estrutural, objetivando a 
correção sistemática de modelos de elementos finitos a partir da confrontação com 
dados experimentais (consultar, por exemplo, o livro de Friswell e Mottershead). 
 Outro fator a ser considerado é que a elaboração de modelos de EF de 
problemas complexos é, na maioria das vezes, um processo interativo, fazendo apelo 
ao conhecimento do Engenheiro acerca do próprio método e do problema em estudo. 
 
 
1.4 - Bibliografia 
 
• Bathe, K. J., Finite Element Procedures in Engineering Analysis, Prentice-Hall, 
1982. 
• Cook, R.D., Concepts and Applications of Finite Element Analysis, John Wiley & 
Sons, 1974. 
• Cook, R., Finite Element Modeling for Analysis, John Wiley & Sons, 1995. 
• Desai, C.S., Abel, J. F., Introduction to the Finite Element Method, Van Nostrand 
Reinhold Company, 1972. 
• Imbert, J. F., Analyse des Structures par Eléments Finis, 3ième édition, Cépauès 
Editions, 1991. 
• Friswell, M.I., Mottershead, J.E., Finite Element Model Updating in Structural 
Dynamics, Kluwer Academic Publishers, 1995. 
• Moaveni, S., Finite Element Analysis. Theory and Application with ANSYS, 
Prentice-Hall, 1999. 
•Zienkiewicz, O.C., The Finite Element Method, 3rd edition, McGraw-Hill Book Co., 
1977. 
•Zienkiewicz, O.C., Morgan, K, Finite Elements and Approximation, John Wiley & 
Sons, 1983 
•Huebner, K.H., Thornton, E.A., The Finite Element Method for Engineers, John 
Wiley & Sons, 1982 
 
 
 
D.A.Rade Formulação do MEF pelo Processo Direto 
15 
 
 
 
CAPÍTULO 2 
 
FORMULAÇÃO DO MÉTODO DOS ELEMENTOS FINITOS 
PELO PROCESSO DIRETO 
 
 No Capítulo 1 havíamos mencionado que as relações matriciais traduzindo o 
comportamento de cada subdomínio (elemento) podiam ser obtidas empregando um 
dos três métodos: Método Direto, Método Variacional e Método dos Resíduos 
Ponderados. 
Neste capítulo vamos ilustrar, com o auxílio de exemplos simples, a utilização 
do Processo Direto, que tem a vantagem de proporcionar facilitar a interpretação 
física para as diversas etapas da elaboração de um modelo de elementos finitos. 
Embora sua utilização seja viável apenas no tratamento de problemas 
unidimensionais simples, os conceitos derivados deste método podem ser estendidos, 
com vantagem, a problemas mais complexos. 
 
2.1 – Análise estática de sistemas compostos por associações de molas 
lineares. 
 
 Um dos problemas mais elementares que podem ser examinados sob oenfoque do MEF é a análise estática do sistema constituído por um conjunto de 
molas lineares, com constantes de rigidez ik , i=1, 2, ..., como aquele exemplificado 
na Figura 2.1. Nesta figura, if e iu designam, respectivamente, as forças externas 
aplicadas e os deslocamentos dos nós. O problema consiste em determinar os 
deslocamentos nodais, dados os valores das forças aplicadas. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Figura 2.1 
 
 
 
 
1k 2k 3k 
nó 1 nó 2 nó 3 nó 4 
2f 3f 4f 1f 
1u 2u 3u 4u 
D.A.Rade Formulação do MEF pelo Processo Direto 
16 
 
 
 
 
2.1.1 – Obtenção das equações de equilibro em nível elementar 
 
 Cada mola é identificada com um elemento finito. Para obter a matriz de 
rigidez para um elemento genérico, este elemento é considerado isoladamente, 
conforme mostrado na Figura 2.2, onde os índices E e D designam as grandezas 
associadas aos nós das extremidades esquerda e direita, respectivamente. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Figura 2.2 
 
 
 Admitindo que todo o conjunto e, portando, cada elemento, encontra-se em 
equilíbrio, podemos escrever: 
 
 EiDi ff −= 
 
 Sabemos ainda que para as molas lineares, o alongamento é proporcional à 
força aplicada em suas extremidades, sendo a constante de proporcionalidade o 
coeficiente de rigidez. Assim, escrevemos: 
 
 ( )EiDiiEi uukf −−= 
 ( )EiDiiDi uukf −= 
 
 As duas equações acima podem ser postas na seguinte forma matricial: 
 
 ⎪⎭
⎪⎬⎫⎪⎩
⎪⎨⎧=⎪⎭
⎪⎬⎫⎪⎩
⎪⎨⎧⎥⎦
⎤⎢⎣
⎡
−
−
D
i
E
i
D
i
E
i
ii
ii
f
f
u
u
kk
kk (2.1) 
 
 A equação acima pode ainda ser escrita na forma compacta: 
 
E
iu Diu 
ik 
E
if 
D
if 
D.A.Rade Formulação do MEF pelo Processo Direto 
17 
 
 ( )[ ] ( ){ } ( ){ }eieiei FUK = i=1,2,... (2.2) 
 
 
onde ( ){ } [ ] ( ){ } [ ]TDiEieiTDiEiei ffFeuuU == são, respectivamente, os vetores de 
deslocamentos nodais e forças nodais em nível elementar e ( )[ ]eiK é denominada 
matriz de rigidez elementar. Sobre esta matriz, podemos observar: 
 
 • o elemento genérico ( )( )mneiK representa a força de reação aplicada no nó m 
quando provocamos um deslocamento unitário no nó n, mantendo bloqueado o nó m. 
 
• é uma matriz simétrica ( )[ ] ( )[ ]Teiei KK = . Levando em conta a interpretação 
dada acima, a simetria da matriz de rigidez traduz o Princípio de Reciprocidade de 
Maxwell-Betti, aplicável a sistemas mecânicos lineares. Isto significa que a força de 
reação que surge no nó m quando provocamos um deslocamento unitário no nó n, 
mantendo bloqueado o nó m é idêntica à força de reação que surge no nó n quando 
provocamos um deslocamento unitário no nó m, mantendo bloqueado o nó n. 
 
• é uma matriz singular (não inversível). Isto porque, como não foram 
introduzidas as condições de contorno em nível elementar, a relação (2.2) deve 
contemplar a existência de uma configuração de equilíbrio sem deformação da mola ( )1+= ii uu e sem o aparecimento de forças nos nós ( )01 == +ii ff . 
 
• é uma matriz semi-definida positiva { } ( )[ ] { } { } { }⎟⎠⎞⎜⎝⎛ ≠∀≥ 00 x,xKx eiT 
 
 
2.1.2 – Obtenção das equações de equilíbrio em nível global 
 
Após a obtenção das equações que descrevem o comportamento de cada 
elemento, isoladamente, devemos considerar o fato que os elementos estão, na 
realidade, interconectados nos pontos nodais. Fisicamente, a interconexão significa 
que deve haver, nos nós compartilhados por mais de uma mola, equilíbrio das forças 
e compatibilidade de deslocamentos. Considerando dois elementos vizinhos, i e i+1, 
ilustrados na Figura 2.3, estas condições são expressas segundo: 
 
11 ++ =+ iEiDi fff (equilíbrio do nó i+1) (2.3) 
 
11 ++ == iEiDi uuu (compatibilidade de deslocamentos no nó i+1) (2.4) 
 
 
 Para imposição destas condições, vamos primeiramente desenvolver a 
equação matricial (2.1) para os elementos i e i+1: 
 
D.A.Rade Formulação do MEF pelo Processo Direto 
18 
 
E
i
D
ii
E
ii fukuk =− (2.5) 
 
D
i
D
ii
E
ii fukuk =+− (2.6) 
 
E
i
D
ii
E
ii fukuk 11111 +++++ =− (2.7) 
 
D
i
D
ii
E
ii fukuk 11111 +++++ =+− (2.8) 
 
Somando as equações (2.7) e (2.8) e introduzindo nas três equações 
resultantes as equações (2.3) e (2.4), obtemos as relações. 
 
E
iii
E
ii fukuk =− +1 (2.9) 
 
( ) iDiiiiiEii fukukkuk =−++− ++++ 1111 (2.10) 
 
D
i
D
iiii fukuk 1111 ++++ =+− (2.11) 
 
Retornando à notação matricial, as equações (2.9) a (2.11) são dispostas sob a 
forma: 
 
⎪⎭
⎪⎬
⎫
⎪⎩
⎪⎨
⎧
=
⎪⎭
⎪⎬
⎫
⎪⎩
⎪⎨
⎧
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢
⎣
⎡
−
−+−
−
+
+
+
+
++
++
D
i
i
E
i
D
i
i
E
i
ii
iiii
ii
f
f
f
u
u
u
kk
kkkk
kk
1
1
1
1
11
11
0
0
 (2.12) 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Figura 2.3 
 
 
E
iu Diu 
ik 
E
if 
E
iu 1+ Diu 1+ 
1+ik 
E
if 1+ Dif 1+ 
D
if 
nó i nó i+1 nó i+1 nó i+2 
D.A.Rade Formulação do MEF pelo Processo Direto 
19 
 
Considerando o exemplo ilustrado na Figura 2.1, aplicando sucessivamente o 
mesmo procedimento aos dois pares de molas que compartilham os nós 2 e 3, 
obtemos a seguinte relação matricial: 
 
 
 
 
⎪⎪⎭
⎪⎪⎬
⎫
⎪⎪⎩
⎪⎪⎨
⎧
=
⎪⎪⎭
⎪⎪⎬
⎫
⎪⎪⎩
⎪⎪⎨
⎧
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
−
−+−
−+−
−
4
3
2
1
4
3
2
1
33
3322
2211
11
00
0
0
00
f
f
f
f
u
u
u
u
kk
kkkk
kkkk
kk
 (2.13) 
 
 A equação acima pode ainda ser escrita na forma compacta: 
 
 ( ) ( ){ } ( ){ }g g gK U F⎡ ⎤ =⎣ ⎦ (2.14) 
 
onde ( ){ } [ ] ( ){ } [ ]1 2 3 4 1 2 3 4eT Tg gU u u u u F f f f f= = são, respectivamente, os 
vetores de deslocamentos e forças em nível global e ( )gK⎡ ⎤⎣ ⎦ é denominada matriz de 
rigidez global. Sobre esta matriz, podemos fazer as mesmas observações 
anteriormente apresentadas para as matrizes de rigidez elementares: 
 
 • o elemento genérico ( )g
mn
K⎡ ⎤⎣ ⎦ representa a força de reação aplicada no nó m 
quando provocamos um deslocamento unitário no nó n, mantendo bloqueado o nó m. 
 
• é uma matriz simétrica ( ) ( )
Tg gK K⎡ ⎤ ⎡ ⎤=⎣ ⎦ ⎣ ⎦ . 
 
• é uma matriz singular. 
 
• é uma matriz semi-definida positiva { } ( ) { } { } { }( )0, 0T gx K x x⎡ ⎤ ≥ ∀ ≠⎣ ⎦ 
 
O processo de obtenção da matriz de rigidez global a partir das matrizes 
elementares e denominado montagem da matriz global. Um procedimento mais 
geral de montagem, bem adaptado à implementação computacional, pode ser 
formulado através da introdução de matrizes de transformação, compreendendo as 
seguintes etapas: 
 
1ª) Para cada elemento, estabelecemos relações matriciais expressando as 
transformações dos deslocamentos nodais em nível elementar com os deslocamentos 
nodais em nível global.. Assim, no exemplo considerado: 
 
D.A.Rade Formulação do MEF pelo Processo Direto 
20 
 
[ ]( )
⎪⎪⎭
⎪⎪⎬
⎫
⎪⎪⎩
⎪⎪⎨
⎧
=⎪⎭
⎪⎬⎫⎪⎩
⎪⎨⎧ ×
4
3
2
1
42
u
u
u
u
T
u
u
iD
i
E
i i=1,2,32ª) Pré-multiplicamos as matrizes de rigidez elementares por [ ]TiT e a pós-
multiplicamos por [ ]iT , para obter as matrizes elementares expandidas, com a 
mesma dimensão da matriz de rigidez global: 
 
 ( )[ ]( ) [ ]( ) ( )[ ]( ) [ ]( )42222444 ×××× = ieiTiei TKTK i=1,2,3 (2.15) 
 
3ª) Adicionamos as matrizes elementares expandidas para obter finalmente a 
matriz de rigidez global. No exemplo: 
 
( ) ( )3
1
g e
i
i
K K
=
⎡ ⎤ ⎡ ⎤=⎣ ⎦ ⎣ ⎦∑ (2.16) 
 
Detalhamos, a seguir, o procedimento de construção da matriz de rigidez 
global para o exemplo considerado: 
 
• elemento nº 1: 
 
⎪⎪⎭
⎪⎪⎬
⎫
⎪⎪⎩
⎪⎪⎨
⎧
⎥⎦
⎤⎢⎣
⎡=⎪⎭
⎪⎬⎫⎪⎩
⎪⎨⎧
4
3
2
1
1
1
0010
0001
u
u
u
u
u
u
D
E
 
 
( )[ ] [ ] ( )[ ] [ ]
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
−
−
=⎥⎦
⎤⎢⎣
⎡⎥⎦
⎤⎢⎣
⎡
−
−
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
==
0000
0000
00
00
0010
0001
00
00
10
01
11
11
11
11
1111
kk
kk
kk
kk
TKTK eTe 
 
 
• elemento nº 2: 
 
⎪⎪⎭
⎪⎪⎬
⎫
⎪⎪⎩
⎪⎪⎨
⎧
⎥⎦
⎤⎢⎣
⎡=⎪⎭
⎪⎬⎫⎪⎩
⎪⎨⎧
4
3
2
1
2
2
0100
0010
u
u
u
u
u
u
D
E
 
 
D.A.Rade Formulação do MEF pelo Processo Direto 
21 
 
( )[ ] [ ] ( )[ ] [ ]
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
−
−=⎥⎦
⎤⎢⎣
⎡⎥⎦
⎤⎢⎣
⎡
−
−
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
==
0000
00
00
0000
0100
0010
00
10
01
00
22
22
22
22
2222 kk
kk
kk
kk
TKTK eTe 
 
• elemento nº 3: 
 
⎪⎪⎭
⎪⎪⎬
⎫
⎪⎪⎩
⎪⎪⎨
⎧
⎥⎦
⎤⎢⎣
⎡=⎪⎭
⎪⎬⎫⎪⎩
⎪⎨⎧
4
3
2
1
3
3
1000
0100
u
u
u
u
u
u
D
E
 
 
( )[ ] [ ] ( )[ ] [ ]
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
−
−=⎥⎦
⎤⎢⎣
⎡⎥⎦
⎤⎢⎣
⎡
−
−
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
==
33
3333
33
3333
00
00
0000
0000
1000
0100
10
01
00
00
kk
kkkk
kk
TKTK eTe 
 
 
( ) ( ) ( ) ( )
1 2 3
g e e eK K K K⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤= + +⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ 
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
−
−+−
−+−
−
=
33
3322
2211
11
00
0
0
00
kk
kkkk
kkkk
kk
 (2.17) 
 
 
2.1.3 – Imposição das condições de contorno 
 
As equações (2.13) devem ser modificadas para levar em conta que os valores 
dos deslocamentos nos nós extremos (1 e 4) podem ser prescritos. Admitamos que as 
condições de contorno sejam as seguintes: 
 
11 uu = (2.18.a) 
 
44 uu = (2.18.b) 
 
onde 1u e 4u são os valores conhecidos. 
 
 Desenvolvendo (2.13) considerando as relações (2.18), obtemos o seguinte 
conjunto de equações: 
 
 
D.A.Rade Formulação do MEF pelo Processo Direto 
22 
 
 
( )
( )
( )
( )
⎪⎪
⎪⎪
⎩
⎪⎪
⎪⎪
⎨
⎧
=+−
+=++−
−=−+
=−
⇒
⎪⎪
⎪⎪
⎩
⎪⎪
⎪⎪
⎨
⎧
=+−
=−++−
=−++
=−
44333
43333222
11232221
12111
44333
34333222
23222111
12111
fukuk
ukfukkuk
ukfukukk
fukuk
fukuk
fukukkuk
fukukkuk
fukuk
 (2.19) 
 
 Determinamos os deslocamentos nodais incógnitos 2u e 3u resolvendo 
simultaneamente a segunda e a terceira equações do sistema (2.19): 
 
⎭⎬
⎫
⎩⎨
⎧
−
−
⎥⎦
⎤⎢⎣
⎡
+−
−+=
⎭⎬
⎫
⎩⎨
⎧⇒
⎭⎬
⎫
⎩⎨
⎧
−
−=
⎭⎬
⎫
⎩⎨
⎧⎥⎦
⎤⎢⎣
⎡
+−
−+ −
433
112
1
322
221
3
2
433
112
3
2
322
221
ukf
ukf
kkk
kkk
u
u
ukf
ukf
u
u
kkk
kkk (2.20) 
 
 
 Os valores das forças de reação nos nós 1 e 4 são finalmente determinados 
através da primeira e quarta equações do sistema (2.19): 
 
⎪⎩
⎪⎨
⎧
+−=
−=
43334
21111
ukukf
ukukf
 (2.21) 
 
 É importante observar que, após a imposição das condições de contorno, a 
matriz de rigidez modificada, deve ser inversível para que a operação indicada em 
(2.20) possa ser efetuada. Este será sempre o caso quando as condições de contorno 
impostas forem suficientes para impedir que o sistema mecânico seja 
cinematicamente variável ou, em outras palavras, que possa se movimentar sem 
que haja deformação de pelo menos uma das molas. 
 Um procedimento sistemático para imposição das condições de contorno e 
cálculo das forças de reação, bem adaptado à implementação computacional, 
consiste em particionar os graus de liberdade em dois conjuntos: graus de liberdade 
livres e graus de liberdade impostos. No nosso exemplo, estes dois conjuntos seriam: 
 
 • graus de liberdade livres: ( ){ } [ ]2 3 TgU u u=A ; ( ){ } [ ]2 3 TgF f f=A 
 
 • graus de liberdade impostos: ( ){ } [ ]1 4 TgU u u=A ; ( ){ } [ ]1 4 TgF f f=A 
 
 Após reordenação das equações, o sistema (2.13) pode ser expresso sob a 
forma: 
 
 
( ) ( )
( ) ( )
( ){ }
( ){ }
( ){ }
( ){ }
g gg g
i
g gg g
i ii ii
U FK K
U FK K
⎡ ⎤ ⎧ ⎫ ⎧ ⎫⎡ ⎤ ⎡ ⎤⎣ ⎦ ⎣ ⎦ ⎪ ⎪ ⎪ ⎪⎢ ⎥ =⎨ ⎬ ⎨ ⎬⎢ ⎥⎡ ⎤ ⎡ ⎤ ⎪ ⎪ ⎪ ⎪⎢ ⎥ ⎩ ⎭ ⎩ ⎭⎣ ⎦ ⎣ ⎦⎣ ⎦
A AAA A
A
 
D.A.Rade Formulação do MEF pelo Processo Direto 
23 
 
 
 Do sistema acima tiramos duas equações matriciais: 
 
 ( ) ( ){ } ( ) ( ){ } ( ){ }g g g g gi iK U K U F⎡ ⎤ ⎡ ⎤+ =⎣ ⎦ ⎣ ⎦AA A A A (2.22) 
 
( ) ( ){ } ( ) ( ){ } ( ){ }g g g g gi ii i iK U K U F⎡ ⎤ ⎡ ⎤+ =⎣ ⎦ ⎣ ⎦A A (2.23) 
 
De (2.22) obtemos os valores dos deslocamentos nodais correspondentes aos 
graus de liberdade livres e, em seguida, a partir de (2.23) calculamos as forças de 
reação correspondentes aos graus de liberdade impostos: 
 
 ( ){ } ( ) ( ){ } ( ) ( ){ }( )1g g g g gi iU K F K U−⎡ ⎤ ⎡ ⎤= −⎣ ⎦ ⎣ ⎦A AA A A (2.24) 
 
( ){ } ( ) ( ){ } ( ) ( ){ }g g g g gi i ii iF K U K U⎡ ⎤ ⎡ ⎤= +⎣ ⎦ ⎣ ⎦A A (2.25) 
 
 
2.2 – Exemplo resolvido utilizando o MATLAB® 
 
 Nesta seção apresentamos a resolução do problema mostrado na Figura 2.1, 
particularizado para a seguinte situação: 
 
 41 20 10k = × N/m 42 30 10k = × N/m 43 50 10k = × N/m 
 
 10002 =f N 1 3 0f f= = N 
 
 041 == uu 
 
 A Figura 2.4 mostra o código implementado em ambiente MATLAB® . A 
solução obtida para o problema é a seguinte: 
 
 32 3,64 10 mu −= × 33 2,73 10 mu −= × 
 
 N37271 ,f −= N72724 ,f −= 
 
 
 
 
 
 
 
 
 
 
 
 
D.A.Rade Formulação do MEF pelo Processo Direto 
24 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Programa para análise estática por EF de sistemas constituídos por molas lineares 
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Domingos Alves Rade 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 
% ENTRADA DE DADOS 
 
% valores dos coeficientes de rigidez das molas 
 k_values = 1e4*[ 20 30 10 ]; % valores em N/m 
% número de nós 
 nb_node = 4; 
% matriz de conectividade 
 mat_conect=[1 2 ; 2 3 ; 3 4 ] 
% condições de contorno nos gdl impostos 
 cond_cont=[1 0 ; 4 0]; 
% forças externas aplicadas nos gld livres 
 forcas_aplic= [ 2 1000 ; 3 0]; % valores em Newtons 
 nb_ele=length(k_values); 
% CONSTRUÇÃO DAS MAT. ELEMENTARES E MONTAGEM DA MATRIZ DE 
% RIGIDEZ GLOBAL 
 K_global=zeros(nb_node); 
 for ii=1:nb_ele 
 K_elementar=k_values(ii)*[1 -1 ; -1 1]; 
 mat_ident=eye(nb_node); 
 mat_transf=[mat_ident(mat_conect(ii,1),:); mat_ident(mat_conect(ii,2),:)]; 
 K_global=K_global+mat_transf'*K_elementar*mat_transf;end 
%IMPOSIÇÃO DAS CONDIÇÕES DE CONTORNO PELO MÉTODO DO 
%PARTICIONAMENTO DA MATRIZ DE RIGIDEZ 
% identificação dos gdl livres e gdl impostos 
 gdl_livres=forcas_aplic(:,1); 
 gdl_impostos = cond_cont(:,1); 
% construção das submatrizes de rigidez 
 K_ll=K_global(gdl_livres,gdl_livres); 
 K_li=K_global(gdl_livres,gdl_impostos); 
 K_ii=K_global(gdl_impostos,gdl_impostos); 
% construção dos vetores de forças nos gdl livres e de deslocamentos nos gdl 
% impostos 
 f_liv=forcas_aplic(:,2); 
 d_imp=cond_cont(:,2); 
% CÁLCULO DOS DESLOCAMENTOS NOS GDL LIVRES E FORÇAS DE REAÇÃO % NOS 
GDL IMPOSTOS 
 d_liv=inv(K_ll)*(f_liv-K_li*d_imp) 
 f_imp=K_li*d_liv+K_ii*d_imp 
 
 
Figura 2.4 
 
D.A.Rade Formulação do MEF pelo Processo Direto 
 
 25
2.3 – Análise estática de treliças planas. 
 
 A idéia de modelar estruturas como uma série de elementos finitos surgiu 
como uma extensão dos métodos matriciais utilizados para análise de treliças. Estas 
estruturas consistem de barras interconectadas em certos pontos considerados como 
rótulas perfeitas, de modo que apenas forças axiais podem ser transmitidas às 
barras, não sendo considerados efeitos de flexão e cisalhamento. Um exemplo de 
treliça é ilustrado na Figura 2.5. 
 
 
 
 
 
 
 
 
 
 
 
Figura 2.5 
 
 
 Na modelagem, cada barra da treliça é tratada com um elemento finito, e os pontos 
de conexão das barras são associados aos nós do modelo de EF. 
 Para determinar as relações forças-deslocamentos para a treliça completa, 
devemos, primeiramente, determinar estas relações para um elemento genérico, 
considerado isoladamente, levando em conta sua orientação, e proceder, em seguida, 
à montagem das equações matriciais que levam em conta a compatibilidade de 
deslocamentos e o equilíbrio de forças nos nós. 
 Considerando a Figura 2.5, que mostra um elemento genérico i, orientado 
segundo um ângulo iα , cujas propriedades relevantes são: comprimento iL , área de 
seção transversal iA e módulo de elasticidade iE . Vamos utilizar os conceitos 
elementares da Resistência dos Materiais para obter as relações forças-
deslocamentos em nível elementar. 
Para cada elemento, definimos dois sistemas de referência ilustrados na 
Figura 2.6: 
 
- o sistema de referência global X-Y, comum a todos os elementos. 
 
- um sistema de referência local x-y, cuja orientação é determinada pela 
orientação de cada barra. 
 
 
 
 
 
 
 
 
D.A.Rade Formulação do MEF pelo Processo Direto 
 
 26
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Figura 2.6 
 
 
 Deve ser notado, na Figura 2.6, que cada um dos nós apresenta duas 
componentes independentes de deslocamento nas direções X e Y, havendo, portanto, 
dois graus de liberdade por nó. Além disso, as forças aplicadas nos nós são 
orientadas na direção da barra. 
 
Da condição de equilíbrio do elemento, temos: 
 
 ( ) ( )E Di ix xF F= − (2.26) 
 
 A relação força-deslocamento, derivada da lei de Hooke para o elemento 
considerado permite escrever: 
 
( ) ( ) ( )Ei iD E xi i ix x i i
F L
L u u
E A
Δ = − = − (2.27) 
 
 ( ) ( ) ( )Di iD E xi i ix x i i
F L
L u u
E A
Δ = − = (2.28) 
 
 
 As duas últimas equações acima podem ser dispostas na seguinte forma 
matricial. 
 
Y 
X 
y 
x 
iL
i
α 
E 
D 
ii A,E 
( ) ( ),D Di ix xu F 
( ) ( ),E Ei ix xu F 
( ) ( )XEiXEi F,u
( ) ( )YEiYEi F,u 
( ) ( ),D Di iY Yu F 
( ) ( ),D Di iX Xu F 
D.A.Rade Formulação do MEF pelo Processo Direto 
 
 27
 
( )
( )
( )
( )
i i i i E E
i ii i x x
D Di i i i i ix x
i i
E A E A
u FL L
E A E A u F
L L
⎡ ⎤− ⎧ ⎫ ⎧ ⎫⎢ ⎥ ⎪ ⎪ ⎪ ⎪⎢ ⎥ =⎨ ⎬ ⎨ ⎬⎢ ⎥ ⎪ ⎪ ⎪ ⎪−⎢ ⎥ ⎩ ⎭ ⎩ ⎭⎣ ⎦
 (2.29) 
 
Devemos observar que, na equação acima os deslocamentos e forças nodais 
estão expressas em termos de suas componentes nas direções dos eixos locais. Como 
cada um dos elementos possui, em geral, orientação diferente da dos demais, antes 
da montagem das equações matriciais globais é necessário expressar as relações 
forças-deslocamentos em termos das componentes no sistema de referência global. 
Isto é feito introduzindo a seguinte transformação para os deslocamentos, que pode 
ser facilmente extraída da Figura 2.4: 
 
( ) ( ) ( )
( ) ( ) ( )
( )
( )
( )
( )
( )
( )
cos sen cos sen 0 0
0 0 cos sencos sen
E
i X
E E E E E
i i i i i i ix X Y x Yi i
D D D D Di ii i i i i i ix X Y x X
D
i Y
u
u u u u u
u u u u u
u
α α α α
α αα α
⎧ ⎫⎪ ⎪⎪ ⎪⎧ ⎫= + ⎪ ⎪⎡ ⎤⎪ ⎪ ⎪ ⎪⇒ =⎨ ⎬ ⎨ ⎬⎢ ⎥⎣ ⎦= + ⎪ ⎪ ⎪ ⎪⎩ ⎭ ⎪ ⎪⎪ ⎪⎪ ⎪⎩ ⎭
 
 
ou, em forma compacta: 
 
 
( )
( )
( )
( )
( )
( )
i
E
i X
E E
i ix Y
D D
i ix X
D
i Y
u
u u
T
u u
u
α
⎧ ⎫⎪ ⎪⎪ ⎪⎧ ⎫ ⎪ ⎪⎪ ⎪ ⎪ ⎪⎡ ⎤=⎨ ⎬ ⎨ ⎬⎣ ⎦⎪ ⎪ ⎪ ⎪⎩ ⎭ ⎪ ⎪⎪ ⎪⎪ ⎪⎩ ⎭
 (2.30) 
 
 
onde: 
 
 
[ ] ⎥⎦⎤⎢⎣⎡= iiii sencossencosT i ααααα 00 00 (2.31) 
 
 
Introduzindo (2.30) em (2.29) e pré-multiplicando a equação resultante por [ ]T
i
Tα , obtemos: 
 
D.A.Rade Formulação do MEF pelo Processo Direto 
 
 28
[ ]T
i
Tα
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
−
−
i
ii
i
ii
i
ii
i
ii
L
AE
L
AE
L
AE
L
AE
[ ]
i
Tα
( )( )( )( ) ⎪⎪⎭
⎪⎪⎬
⎫
⎪⎪⎩
⎪⎪⎨
⎧
Y
D
i
X
D
i
Y
E
i
X
E
i
u
u
u
u
= [ ]T
i
Tα
( )
( )
E
i x
D
i x
F
F
⎧ ⎫⎪ ⎪⎨ ⎬⎪ ⎪⎩ ⎭
 (2.32) 
 
 Com base na Figura 2.4, podemos desenvolver o membro do lado direito de 
(2.32) sob a forma: 
 
[ ]T
i
Tα
( )
( )
E
i x
D
i x
F
F
⎧ ⎫⎪ ⎪⎨ ⎬⎪ ⎪⎩ ⎭
 =
⎥⎥
⎥⎥
⎦
⎤
⎢⎢
⎢⎢
⎣
⎡
i
i
i
i
cos
cos
sen
cos
α
α
α
α
0
0
0
0 ( )
( )
E
i x
D
i x
F
F
⎧ ⎫⎪ ⎪⎨ ⎬⎪ ⎪⎩ ⎭
=
( )( )( )( ) ⎪⎪⎭
⎪⎪⎬
⎫
⎪⎪⎩
⎪⎪⎨
⎧
=
⎪⎪⎭
⎪⎪⎬
⎫
⎪⎪⎩
⎪⎪⎨
⎧
Y
D
i
X
D
i
Y
E
i
X
E
i
i
D
i
i
D
i
i
E
i
i
E
i
F
F
F
F
senF
cosF
senF
cosF
α
α
α
α
 (2.33) 
 
 Assim, a equação (2.32) assume a seguinte forma, representando a relação 
forças-deslocamentos expressa no sistema de referência global: 
 
 
⎥⎥
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢⎢
⎢
⎣
⎡
−
−−
−−
−−
iiiiii
iiiiii
iiiii
iiiiii
i
ii
sensencossensencos
sencoscossencoscos
sensencossensencos
sencoscossencoscos
L
AE
αααααα
αααααα
αααααα
αααααα
22
22
22
22 ( )( )( )( ) ⎪⎪⎭
⎪⎪⎬
⎫
⎪⎪⎩
⎪⎪⎨
⎧
Y
D
i
X
D
i
Y
E
i
X
E
i
u
u
u
u
=
( )( )( )( ) ⎪⎪⎭
⎪⎪⎬
⎫
⎪⎪⎩
⎪⎪⎨
⎧
Y
D
i
X
D
i
Y
E
i
X
E
i
F
F
F
F
 
 
(2.34) 
 
ou: 
 
 ( )[ ] ( ){ } ( ){ }eieiei FUK = i=1,2,... (2.35) 
 
onde: 
 
 
( )[ ] =eiK
⎥⎥
⎥⎥
⎥
⎦
⎤
⎢⎢
⎢⎢
⎢
⎣
⎡
−
−−
−−
−−
iiiiii
iiiiii
iiiii
iiiiii
i
ii
sensencossensencos
sencoscossencoscos
sensencossensencos
sencoscossencoscos
L
AE
αααααα
αααααα
αααααα
αααααα
22
22
22
22
 (2.36) 
 
é a matriz de rigidez elementar do elemento de treliça. 
 Devemos observar que cada elemento possui quatro graus de liberdade (2 
graus de liberdade por nó, correspondendo aos deslocamentos nas direções X e Y. 
 Uma vez desenvolvida a expressão geral da matriz de rigidez para um 
elemento de barra genérico, podemos utilizar os procedimentos já detalhados na 
Seção 2.1 para a montagem da matriz de rigidez global – traduzindo a 
compatibilidadede deslocamentos e equilíbrio de forças nos nós – e a imposição das 
D.A.Rade Formulação do MEF pelo Processo Direto 
 
 29
condições de contorno. No que segue, estas etapas serão detalhadas a partir de um 
exemplo numérico. 
 
 
2.4 – Exemplo resolvido usando o MATLAB®. 
 
 
 Consideremos a treliça plana constituída de três barras, ilustrada na 
Figura 2.7(a), cujas características físicas e geométricas são dadas na Tabela 2.1. 
Interessa-nos determinas os deslocamentos nodais, as reações nos apoios e as 
tensões em cada uma das barras. 
 A Figura 2.7(b) mostra a numeração escolhida para os nós e elementos, bem 
como a indicação das componentes das forças e deslocamentos nodais no sistema de 
referência global adotado. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
(a) 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
(b) 
 
Figura 2.7 
 
 
45o 100 KN 
Y 
X 
3
2
1 
1
2
3 ( )Xu3 
( )Yu3 
( )Xu2 
( )Yu2 
( )XF3 
( )YF3 
( )YF2 
( )XF2 
( )Xu1 
( )Yu1 ( )YF1 
( )XF1 
D.A.Rade Formulação do MEF pelo Processo Direto 
 
 30
Tabela 2.1 
Elemento A(cm2) E(N/m2) L(m) Conectividade 
(nós) 
α (graus) 
 
1 
 
 
32,3 
 
6,9×1010 
 
2,54 
 
2,3 
 
90 
 
2 
 
 
38,7 
 
20,7×1010 
 
2,54 
 
1,2 
 
0 
 
3 
 
 
25,8 
 
20,7×1010 
 
3,59 
 
1,3 
 
135 
 
 O programa MATLAB® utilizado e os resultados obtidos são apresentados na 
Figura 2.5. 
 Deve ser observado que, após o cálculo das forças e deslocamentos nodais, as 
tensões nas barras foram calculadas de acordo com o seguinte procedimento: 
 Utilizando a equação (2.28), escrevemos: 
 
 
( ) ( ) ( )Di D Ex ii i ix xi i
F E u u
A L
σ ⎡ ⎤= = −⎢ ⎥⎣ ⎦ 
 
 Os valores de ( )Di xu e ( )Ei xu são calculados a partir dos deslocamentos nodais 
expresso no sistema de referência global utilizando a equação (2.30). 
Os resultados obtidos para este exemplo são os seguintes: 
 
• Deslocamentos 
 
( )1 0,6858mmXu = − ( )2 0,9100 mmXu = − ( )2 0,8059 mmYu = − 
 
• Reações de apoio 
 
( )2 70,711KNYF = ( )3 70,711KNXF = ( )3 0YF = 
 
• Tensões nas barras 
 
1 21,819MPaσ = 2 18,271MPaσ = 3 38,760MPaσ = − 
 
 
 
 
 
 
 
 
 
D.A.Rade Formulação do MEF pelo Processo Direto 
 
 31
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Programa para análise estática por EF de treliças planas 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Domingos Alves Rade Última modificação: 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% ENTRADA DE DADOS 
 % propriedades físicas e geométricas das barras [A(m^2) E(N/m^2) L(m) alfa(grau)] 
 propriedades = [ 32.3e-4 6.9e10 2.54 90 
 38.7e-4 20.7e10 2.54 180 
 25.8e-4 20.7e10 2.59 135]; 
 % número total de graus de liberdade 
 nb_gdl = 6; 
 % matriz de conectividade (graus de liberdade ordenados segundo: 
 % 1: nó 1, direção x 
 % 2: nó 1, direção y 
 % 3: nó 2, direção x 
 % 4: nó 2, direção y 
 % 5: nó 3, direção x 
 % 6: nó 3, direção y 
 mat_conect=[ 3 4 5 6 
 1 2 3 4 
 1 2 5 6]; 
 % condições de contorno nos gdl impostos 
 cond_cont=[2 0 
 5 0 
 6 0]; 
 % forças externas aplicadas nos gld livres(valores em Newtons) 
 forcas_aplic= [ 1 0 
 3 -100000*cos(pi/4) 
 4 -100000*sin(pi/4)]; 
 [nb_ele,dummy]=size(propriedades); 
% CONSTRUÇÃO DAS MATRIZES ELEMENTARES E MONTAGEM DA MATRIZ DE RIGIDEZ GLOBAL 
K_global=zeros(nb_gdl); 
for ii=1:nb_ele 
 coef=propriedades(ii,1)*propriedades(ii,2)/propriedades(ii,3); 
 alfa=propriedades(ii,4)*pi/180; 
 K_elementar=coef*[1 -1 ; -1 1]; 
 mat_transf=[cos(alfa) sin(alfa) 0 0 
 0 0 cos(alfa) sin(alfa)]; 
 K_elementar=mat_transf'*K_elementar*mat_transf; 
 mat_ident=eye(nb_gdl); 
 mat_transf1=[mat_ident(mat_conect(ii,1),:); mat_ident(mat_conect(ii,2),:); ... 
 mat_ident(mat_conect(ii,3),:); mat_ident(mat_conect(ii,4),:)]; 
 K_global=K_global+mat_transf1'*K_elementar*mat_transf1; 
end 
% IMPOSIÇÃO DAS CONDIÇÕES DE CONTORNO PELO MÉTODO DO PARTICIONAMENTO 
 %identificacao dos gdl livres e gdl impostos 
 gdl_livres=forcas_aplic(:,1); 
 gdl_impostos = cond_cont(:,1); 
 % construção das submatrizes de rigidez 
 K_ll=K_global(gdl_livres,gdl_livres); 
 K_li=K_global(gdl_livres,gdl_impostos); 
 K_ii=K_global(gdl_impostos,gdl_impostos); 
 % construção dos vetores de forças nos gdl livres e de deslocamentos nos gdl 
impostos 
 f_liv=forcas_aplic(:,2); 
 d_imp=cond_cont(:,2); 
 % CÁLCULO DOS DESLOCAMENTOS NOS GDL LIVRES E FORÇAS DE REAÇÃO NOS GDL IMPOSTOS 
 disp('Deslocamentos nos gdl livres (em milímetros)'); 
 d_liv=inv(K_ll)*(f_liv-K_li*d_imp)*1000 
 disp('Reacoes de apoio (em kN) '); 
 f_imp=(K_li'*d_liv/1000+K_ii*d_imp)/1000 
 des_global=zeros(nb_gdl,1); 
 des_global(cond_cont(:,1))=cond_cont(:,2); 
 des_global(forcas_aplic(:,1))=d_liv/1000; 
 % CÁLCULO DAS TENSÕES NAS BARRAS 
for ii=1:nb_ele 
 alfa=propriedades(ii,4)*pi/180; 
D.A.Rade Formulação do MEF pelo Processo Direto 
 
 32
 mat_transf=[cos(alfa) sin(alfa) 0 0 
 0 0 cos(alfa) sin(alfa)]; 
 des_local=mat_transf*des_global(mat_conect(ii,:)); 
 sigma(ii)=propriedades(ii,2)/propriedades(ii,3)*(des_local(2)-des_local(1)); 
end 
 disp('Tensões nas barras (em MN/m^2) '); 
 sigma/1e6 
 
Figura 2.8 
 
 
 
D.A.Rade Formulação do MEF pelo Processo Direto 
33 
 
 
2.5 - Aplicação do processo direto a outros tipos de problemas unidimensionais 
 
A formulação do MEF pelo processo direto, ilustrado na seção 2.1 para 
sistemas formados por associações de molas lineares pode ser utilizado para a 
resolução de vários outros tipos de problemas unidimensionais em regime 
permanente nas áreas de Mecânica dos Sólidos, Mecânica dos Fluidos e 
Transferência de Calor. Para cada um destes tipos de problemas o processo 
consiste em determinar, a partir da aplicação das leis Físicas pertinentes, as 
equações equivalentes às equações de equilíbrio das molas. A partir daí, o 
procedimento de montagem e resolução das equações globais segue o mesmo 
procedimento estudado anteriormente. A obtenção das equações equivalentes é 
ilustrada a seguir. 
 
 
2.5.1 - Barras em solicitação axial 
 
 
Consideramos aqui problemas formados por barras solicitadas por forças 
axiais, como exemplificado na Figura 2.9, na qual, para cada segmento da barra, 
, ,i i iE A l indicam, respectivamente, o módulo de elasticidade, área de seção 
transversal e comprimento. ei iN U designam, respectivamente, as forças axiais e 
os deslocamentos aplicados nos nós que delimitam cada segmento. 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Figura 2.9 
 
1 (1) 
(2) 
(3) 
N1 
U1 
N2 
U2 
N3 
U3 
N4 
U4 
E1, A1, l1 
E2, A2, l2 
E3, A3, l3 
2 3 4 
(i) 
NiE 
Ui E 
Ni D 
Ui D 
Ei, Ai, li 
i i+1 
ki= EiAi/li 
fiE 
ui E 
fi D 
ui D 
D.A.Rade Formulação do MEF peloProcesso Direto 
34 
 
Da Resistência dos Materiais sabemos que para cada segmento em 
equilíbrio, considerado isoladamente, admitindo comportamento linear elástico, 
valem as seguintes relações: 
 
 D D Ei ii i i i i i i i
i
E AN A E A U U
l
     
 
 E D E Di ii i i i
i
E AN N U U
l
    
 
As duas equações acima por se postas na seguinte forma matricial: 
 
1 1
1 1
E E
i ii i
D D
i i i
U NE A
l U N
               
 (2.38) 
 
 Comparando as equações (2.1) e (2.38) concluímos que cada segmento da 
barra solicitada axialmente pode ser considerado como uma mola linear cuja 
constante de rigidez equivalente é dada por 
 
i i
i
i
E Ak
l
 
 
Desta forma, o procedimento de montagem das equações globais e 
imposição das condições de contorno e carregamentos é feito da mesma forma 
adotada para os sistemas de molas considerados na Seção 2.1. 
 
2.5.2 Eixos sujeitos a torção 
 
Consideramos aqui sistemas formados por eixos solicitados por momentos 
torsores (torques), conforme mostrado na Figura 2.10, na qual, para cada 
segmento, , ,i i iG J l indicam, respectivamente, o módulo de cisalhamento, momento 
polar de inércia da seção transversal e comprimento. Além disso, ei i  
designam, respectivamente, os torques aplicados e ângulos de torção das seções 
transversais nos nós que delimitam cada segmento. 
 
 
 
 
 
 
 
 
 
 
 
 
D.A.Rade Formulação do MEF pelo Processo Direto 
35 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Figura 2.10 
 
Da Resistência dos Materiais sabemos que para cada segmento em 
equilíbrio, considerado isoladamente, admitindo comportamento linear elástico, 
valem as seguintes relações: 
 
 D D Ei ii i i
i
G J
l
    
 
 E D E Di ii i i i
i
G J
l
       
 
As duas equações acima por se postas na seguinte forma matricial: 
 
1 1
1 1
E E
i ii i
D D
i i i
G J
l
 
 
               
 (2.39) 
 
 Comparando as equações (2.1) e (2.39) concluímos que cada segmento da 
barra solicitada em torção pode ser considerado como uma mola linear cuja 
constante de rigidez equivalente é dada por 
 
i i
i
i
G Jk
l
 
 
Para o caso de eixos de seção circular de raio iR , o momento polar de 
inércia é dado por: 
 
4
2
i
i
RJ  
 
G1, J1, l1 G2, J2, l2 
T1 ,1 T2 ,2 
T3 , 3 
(2) (1) 
ki = GiJi,/li 
TiE 
iE 
TiD 
iD 
Gi, Ji, li 
TiE , i E TiD 
(1) 
, i D 
D.A.Rade Formulação do MEF pelo Processo Direto 
36 
 
 
O procedimento de montagem das equações globais e imposição das 
condições de contorno e carregamentos é feito da mesma forma adotada para os 
sistemas de molas considerados na Seção 2.1. 
 
 
2.5.3 Problemas de escoamento em regime permanente 
 
 Nesta seção tratamos o problema de escoamento unidimensional em 
regime permanente em tubulações que podem ser consideradas constituída pela 
associação de vários segmentos, como ilustrado na Figura 2.11(a) que mostra a 
vazão total Q e as vazões que percorrem os segmentos, indicadas por qi , i=1,2, ... 
O problema consiste em determinar os valores destas vazões, bem como das 
pressões nos nós que delimitam os segmentos, os quais são também indicados na 
Figura 2.11. 
 
 
 
 
 
 
 
 
 
 
 
 
 
Figura 2.11(a) 
 
 
 A Figura 2.11(b) mostra um segmento arbitrário, delimitado por dois nós E 
e D, nos quais estão indicadas as respectivas pressões Eip e Dip e vazões Eiq e Diq . 
 
 
 
 
 
 
 
 
 
Figura 2.11(b) 
 
 
Em virtude da hipótese de escoamento em regime permanente, para cada 
segmento i deve-se ter: 
 
A 
E H
D Q Q 
q1 
q2 q2
q3
q4
F G
B Cq1
q1
q3
1/Ri 
E
iq
D
iq
E
ip
D
ip
i 
E
iq
E
ip D
ip 
D
iq 
D.A.Rade Formulação do MEF pelo Processo Direto 
37 
 
D E
i i iq q q  (2.41) 
 
 Além disso, sabe-se da Mecânica dos Fluidos que: 
 
D E
i i i ip p R q  (2.42) 
 
onde Ri é o coeficiente de resistência, que depende essencialmente do tipo de 
fluido, da geometria da seção transversal e do acabamento da superfície interna 
do duto. 
 Combinando as equações (2.41) e (2,42) escrevemos: 
 
 
 1E E Di i i
i
q p p
R
 
 
 
 1D E Di i i
i
q p p
R
   
 
ou: 
 
1 11
1 1
E E
i i
D D
ii i
q p
Rq p
             
 (2.43) 
 
 
 Como nos casos anteriores, a comparação de (2.1) e (2.43) permite concluir 
que o problema de escoamento unidimensional em regime permanente pode ser 
modelado da mesma forma que o sistema de molas lineares considerado na Seção 
2.1, sendo a rigidez de mola equivalente dada por 1/ Ri. 
 
 
2.5.4 Transferência de calor unidimensional em regime permanente 
 
 
Tratamos aqui problema de transferência de calor unidimensional em 
regime permanente em uma barra constituída pela associação de vários 
segmentos, como ilustrado na Figura 2.12(a), onde Ki, Ai e li designam, 
respectivamente, o coeficiente de condutividade térmica, a área da seção 
transversal e o comprimento dos segmentos que compõem a barra. O problema a 
ser resolvido consiste em determinar a distribuição de temperatura T(x) para um 
conjunto de condições de contorno impostas nas extremidades 1 e 4. 
A Figura 2.12(b) mostra um segmento genérico i da barra e as grandezas a 
ele associadas. 
 
 
 
 
 
 
D.A.Rade Formulação do MEF pelo Processo Direto 
38 
 
 
 
 
 
 
 
 
 
 
 
Figura 2.12(a) 
 
 
 
 
 
 
 
 
 
Figura 2.12(b) 
 
 
 Recordamos abaixo as equações pertinentes aos dois mecanismos de 
transferência de calor (condução e convecção) que consideramos estar envolvidos 
no problema em questão. 
 
 
 
Condução (Lei de Fourier): 
 
 
 
 
 
i i i
dTq K A
dx
  (2.44.a)
 
 
Convecção: 
 
 
 
 
 
  i i i Pq h A T T  (2.44.b)
 
 
qi 
D 
Ki, Ai, li EiT
E 
i 
D
iT
E
iq
D
iq
E D
i iT T
Ki 
TiE 
qiE 
TiD 
qiD 
Q1 
T1 
Q4 
K1, A1, l1 
K2, A2, l2 
K3, A3, l3 
1 3 4 
T2 T3 T4 
2 
x 
T 
TP 
qi 
D.A.Rade Formulação do MEF pelo Processo Direto 
39 
 
Análise em nível elementar 
 
E
i iq q Di iq q  
 
 E Di ii i i
i
K Aq T T
l
  Di iq q  
 
Assim, para o elemento i: 
 
 
1 1
1 1
E E
i ii i
D D
ii i
q TK A
lq T
             
 (2.45) 
 
 Como nos casos anteriores, a comparação de (2.1) e (2.44) mostra que o 
problema de transferência de calor unidimensional em regime permanente pode 
ser modelado da mesma forma que o sistema de molas lineares considerado na 
Seção 2.1, sendo a rigidez de mola equivalente dada por KiAi/li. 
 Uma especificidade a ser considerada neste problema é que são possíveis 
combinações de três tipos de condição de contorno aplicadas nas extremidades da 
barra: 
 
1) valor do calor imposto; 
2) valor da temperaturaimposta; 
3) condição de contorno do tipo convecção com o valor da temperatura T 
imposto; 
 
A condição de contorno do tipo temperatura imposta deve ser tratada da 
mesma forma que os deslocamentos impostos, apresentada na Seção 2.1.3 e a 
condição de valor imposto é análoga à aplicação de uma força de valor conhecido 
nos problemas estruturais. 
O tratamento das condições de contorno do tipo convecção requer algumas 
transformações adicionais das equações, uma vez que este tipo de condição de 
contorno envolve as temperaturas das extremidades da barra, conforme pode ser 
visto na Equação (2.44.b). 
Considerando, para exemplificação, a barra discretizada com três 
elementos finitos, as equações globais resultantes do procedimento de montagem, 
tomam a forma: 
 
D.A.Rade Formulação do MEF pelo Processo Direto 
40 
 
1 1 1 1
1 1
1 11 1 1 1 2 2 2 2
1 1 2 2 2
33 3 3 32 2 2 2
4 42 2 3 3
3 3 3 3
3 3
0 0
0
0
0
0
0 0
K A K A
l l
q TK A K A K A K A
l l l l T
TK A K AK A K A
q Tl l l l
K A K A
l l
                                       
 (2.46) 
 
Admitamos que na extremidade 1 tenhamos uma condição de contorno do 
tipo temperatura imposta e na extremidade 4 tenhamos uma condição de 
contorno de convecção, as quais são expressas segundo : 
 
1 1T T 
  4 4 3 4q h A T T  
 
 Introduzindo estas expressões em (2.46), obtemos: 
 
1 1 1 1
1 1
1 11 1 1 1 2 2 2 2
1 1 2 2 2
33 3 3 32 2 2 2
4 3 2 2 3 3 4
3 3 3 3
4 3
3 3
0 0
0
0
0
0
0 0
K A K A
l l
q TK A K A K A K A
l l l l T
TK A K AK A K A
h A T l l l l T
K A K A h A
l l

                                        
 
 
O sistema de equações acima pode ser decomposto em dois conjuntos de 
equações cuja resolução fornece o valor das temperaturas desconhecidas T2, T3 e 
T4 e o calor transferido através da seção transversal no nó 1, segundo: 
 
1 1 2 2 2 2
1 2 2
2
3 3 3 32 2 2 2
3
2 2 3 3
4 3 4
3 3 3 3
4 3
3 3
0
0
0
0
K A K A K A
l l l
T
K A K AK A K A T
l l l l
h A T T
K A K A h A
l l

                              
 
 
 
D.A.Rade Formulação do MEF pelo Processo Direto 
41 
 
1
21 1 1 1
1
1 1 3
4
0 0
T
TK A K Aq
l l T
T
              
 
D.A. Rade Fundamentos de Cálculo Variacional 
 42
CAPÍTULO 3 
 
FUNDAMENTOS DE CÁLCULO VARIACIONAL 
 
É possível, e muito conveniente, considerar o procedimento de obtenção de 
modelos de elementos finitos como um processo variacional. Através deste 
processo procura-se estabelecer as equações do movimento ou de equilíbrio em 
nível elementar e em nível global mediante a busca dos pontos estacionários (ou 
pontos críticos, ou pontos extremos) de certas funções, denominadas funcionais, 
que são funções escalares de funções que definem a distribuição da variável de 
campo sobre o domínio considerado. Neste contexto o Cálculo Variacional 
constitui o conjunto de técnicas matemáticas destinadas à determinação dos 
pontos estacionários de funcionais. 
 Os chamados Princípios Variacionais são leis que estabelecem que a 
configuração de equilíbrio estático ou de equilíbrio dinâmico de um dado sistema 
mecânico devem corresponder a pontos estacionários de determinados funcionais 
associados ao problema. Em diversos problemas da Mecânica do Contínuo, os 
funcionais representam a energia do sistema. As equações de equilíbrio ou do 
movimento podem então ser obtidas desenvolvendo, matematicamente, as 
condições de estacionaridade destes funcionais. 
 Neste Capítulo, apresentaremos os fundamentos do Cálculo Variacional e 
ilustraremos os desenvolvimentos apresentados através de exemplos baseados no 
uso de Princípios Variacionais da Mecânica, notadamente o Princípio da Energia 
Potencial Estacionária e o Princípio de Hamilton. Tais princípios serão 
apresentados detalhadamente no Capítulo 5. 
 
 
3.1 Funcionais envolvendo uma variável dependente e uma variável 
independente 
 
 Consideremos inicialmente o problema de determinação dos pontos 
estacionários do seguinte funcional: 
 
       dxxu,xu,xFuI x
x  21 (3.1) 
 
 Admitiremos ainda que a função u deva satisfazer as seguintes condições 
de contorno: 
   11 uxu  ;   22 uxu  (3.2) 
 
Nas equações acima, x designa a variável independente,  xu é a variável 
dependente, que representa a variável de campo associada ao problema, u indica 
a derivada de  xu em relação à variável independente,     xu,xu,xF  indica 
genericamente uma relação funcional envolvendo x,  xu e  xu e 1u e 2u 
D.A. Rade Fundamentos de Cálculo Variacional 
 43
designam valores prescritos, assumidos pela variável dependente nos extremos 
do intervalo em que o problema é definido. 
Designando doravante por u=u(x) a função que corresponde à solução 
estacionária exata do funcional definido em (3.1), o procedimento variacional 
consiste em definir as chamadas funções de comparação, representadas sob a 
forma: 
      xuxuxu~  21 xxx  (3.3.a) 
 
com: 
    xxu   (3.3.b) 
 
onde 1 é um parâmetro escalar arbitrário e  x é uma função contínua 
arbitrária definida no intervalo 21 xxx  , que assume valor nulo nos extremos 
deste: 
 
   01 x   02 x (3.4) 
 
Em decorrência de (3.4), as funções de comparação satisfazem as condições 
de contorno do problema: 
   11 uxu~  ;   22 uxu~  (3.5) 
 
 Em (3.3.b),  x é entendida como uma função arbitrária que representa 
uma perturbação infinitesimal acrescentada à solução exata u(x) para cada valor 
da variável independente x, conforme ilustra a Figura 3.1, onde ficam 
evidenciadas as condições (3.4) e (3.5). 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
Figura 3.1 
 
x 2x 1x
u, u 
 xy   
 xu 
 xu~ 
1u 
2u 
D.A. Rade Fundamentos de Cálculo Variacional 
 44
 Na seqüência do processo variacional, devemos introduzir a solução 
variada (3.3) no funcional dado por (3.1) e impor a condição de estacionariedade 
do funcional, que determina que sua variação I , associada a qualquer variação 
arbitrária u deve ser nula. Esquematicamente, este processo é ilustrado como 
segue: 
 
 
0 IIIFFuu  
 
 
Buscaremos, primeiramente, expressar a variação da função 
    xu,xu,xF  , designada por     xu,xu,xF  , resultante da substituição de  xu 
por    xuxu  . Assim, definimos: 
            xu,xu,xFuu,uu,xFxu,xu,xF   (3.6) 
 
 Em decorrência de (3.3.b), temos: 
 
  xu   (3.7) 
 
 Introduzindo (3.3.b) e (3.7) em (3.6), obtemos: 
               xu,xu,xFxu,xu,xFxu,xu,xF   (3.8) 
 
 Negligenciando os termos de ordem igual e superior à segunda em  , 
procedemos à expansão do primeiro termo do lado direito de (3.8) em série de 
Taylor, obtendo: 
 
             

 
 x
u
Fx
u
Fxu,xu,xFxu,xu,xF  (3.9) 
 
 Introduzindo (3.9) em (3.8), temos: 
 
       F FF x,u x ,u x x x
u u
               (3.10) 
 
A variação do funcional I correspondente à variação da função F é definida 
segundo: 
 
   FIFFII   =           dxxu,xu,xFxu,xu,xFx
x 21  
      dxxu,xu,xFx
x  21      dxxu,xu,xFxx  21  (3.11) 
 
 Introduzindo (3.10) em (3.11), temos: 
 
D.A. Rade Fundamentos de Cálculo Variacional 
 45
        2
1
x
x
F FI x,u x ,u x x x dx
u u
               (3.12) 
 
A condição necessária para a estacionaridade do funcional I é que sua 
variação deve anular-se, ou seja,      0 xu,xu,xI . Assim, com base em (3.12), 
devemos ter: 
 
   2
1
0
x
x
F Fx x dx
u u
        (3.13) 
 
Procedemos, em seguida, à integração por partes da segunda parcela do 
integrando: 
 
       





 2
1
2
1
2
1
x
x
x
x
x
x
dxx
u
F
dx
dx
u
Fdxx
u
F  (3.14) 
 
Introduzindo (3.14) em (3.13), temos: 
 
   22
1 1
0
xx
x x
F F d Fx x dx
u u dx u
                    (3.15) 
 
Em decorrência de (3.4), a primeira parcela de (3.15) se anula: 
 
  02
1




 x
x
x
u
F  , (3.16) 
 
de modo que: 
 
 2
1
0
x
x
F d F x dx
u dx u
          (3.17)
 
Neste ponto, evocamos o seguinte lema: 
 
Lema: Se 1x e 2x ( 1x < 2x ) são constantes e G(x) é uma função contínua de x 
no intervalo  21; xx e se 
 
    02
1
 dxxGx
x
x
 
 
para toda e qualquer função continuamente diferenciável  x que satisfaz 
D.A. Rade Fundamentos de Cálculo Variacional 
 46
 
   01 x   02 x , 
 
conclui-se que   0xG identicamente no intervalo  21; xx . 
 
Com base neste lema podemos afirmar que a equação (3.17) se verifica se: 
 
0F d F
u dx u
       (3.18) 
 
 A equação (3.18) é denominada Equação de Euler-Lagrange associada ao 
funcional dado por (3.1). Por outro lado, a equação (3.16) indica todas as possíveis 
combinações de condições de contorno impostas nos extremos do intervalo  21; xx . 
No caso do funcional analisado, as condições de contorno que satisfazem (3.16) 
são: 
 
 • Em x= 1x :   01 x (3.18.a) 
 
0
1


xxu
F (3.18.b) 
 
 • Em x= 2x :   02 x (3.18.c) 
 
0
2


xxu
F (3.18.d) 
 
 Neste ponto, introduzimos a seguinte classificação para as condições de 
contorno: 
 
a) condições de contorno geométricas ou de Dirichlet: são aquelas que 
envolvem diretamente a função de comparação e suas derivadas. 
 
b) condições de contorno naturais ou de Neumann: são aquelas 
representadas por expressões que envolvem as derivadas da função 
    xu,xu,xF  em relação a um de seus argumentos ou derivadas de 
 x . 
 
De acordo com estas definições, as condições de contorno indicadas pelas 
equações (3.18.a) e (3.18.c) são condições de contorno geométricas, ao passo que 
aquelas indicadas pelas equações equações (3.18.b) e (3.18.d) são condições de 
contorno naturais. 
 
 
 
 
D.A. Rade Fundamentos de Cálculo Variacional 
 47
Observações: 
 
• a equação de Euler-Lagrange (3.18) apresenta-se sob a forma de uma 
equação diferencial, cuja resolução, complementada pelas condições de contorno 
(3.2), fornece a função u(x) que corresponde ao ponto estacionário do funcional 
(3.1). 
 
• para um determinado ponto estacionário do funcional, pode-se 
determinar sua natureza (máximo ou mínimo), a partir da análise do sinal da 
segunda variação do funcional I. A segunda variação é obtida pela expansão de  uuI  em torno de  uI utilizando uma série de Taylor limitada aos termos de 
segunda ordem em u . Para o funcional em questão, escrevemos: 
 
        uu
uu
Fu
u
Fu
u
Fu
u
Fu
u
Fu,u,xFuu,uu,xF 




 
2
2
2
2
2
2
2
2
1
2
1
 
ou: 
 
       u,u,xFu,u,xFu,u,xFuu,uu,xF  2 (3.19) 
 
com: 
 
  u
u
Fu
u
Fu,u,xF 

  (3.20) 
 
      uu
uu
Fu
u
Fu
u
Fu,u,xF 


 
2
2
2
2
2
2
2
2
2
1
2
1 (3.21) 
 
Introduzimos agora a segunda variação do funcional I: 
 
    dxuu,uu,xFuuI x
x  21  (3.22) 
 
Introduzindo (3.19) em (3.22), temos: 
 
         2
1
2
1
2
1
22
x
x
x
x
x
x
dxu,u,xFdxu,u,xFdxu,u,xFuuI  (3.23) 
 
Levando em conta (3.1) e (3.11), reescrevemos a equação acima sob a 
forma: 
 
     IIuIuuI 2  (3.24) 
 
onde: 
 
D.A. Rade Fundamentos de Cálculo Variacional 
 48
   2
1
222
x
x
dxu,u,xFI  (3.25) 
 
é a segunda variação do funcional I. 
 
Introduzindo (3.21) em (3.25) e levando em conta (3.3.b), escrevemos: 
 
  


 



2
1
2
2
2
2
2
2
2
22 2
2
1 x
x
dx
uu
F
u
F
u
FI  (3.26) 
 
 Pode-se mostrar que um dado ponto estacionário do funcional I é um ponto 
de mínimo se 02 I e um ponto de máximo se 02 I para qualquer função 
arbitrária  . 
 
 
Exemplo 1: O problema de equilíbrio de uma corda tensionada 
 
Consideremos a corda elástica de comprimento L ilustrada na Figura 3.2, 
sujeita a uma tração axial constante T e a um carregamento transversal 
distribuído (força por unidade de comprimento) q(x). Nesta mesma figura u=u(x) 
indica o campo de deslocamento transversais correspondente à posição de 
equilíbrio da corda. 
Pode-se mostrar que a energia potencial total do sistema, incluindo a 
energia de deformação da corda e o trabalho da carga distribuída é dada pelo 
seguinte funcional: 
 
     dxquuTuI L   0
2
2
 (a) 
 
 
 
 
 
 
 
 
 
 
 
 
 
Figura 3.2 
 
 
 Confrontando a equação (a) com a equação (3.1), reconhecemos: 
L
q(x) 
x 
y 
u(x) 
D.A. Rade Fundamentos de Cálculo Variacional 
 49
 
     quuTu,uF  2
2
 (b) 
 
 Após o cálculo das derivadas parciais necessárias, as equações (3.18) e 
(3.16), nos fornecem, respectivamente, a equação de Euler-Lagrange e os termos 
de contorno do problema: 
  
T
xqu  (c) 
 
  00  LuT  (d) 
 
 Admitindo    constanteqxq  , integrando duas vezes (c) e introduzindo as 
condições de contorno determinadas pelos vínculos mostrados na Figura 3.1: 
 
   00 u   0Lu , 
 
obtemos a seguinte expressão para o campo de deslocamentos transversais da 
corda em equilíbrio: