Método Dos ElementosFinitos em Análise deEstruturas

Luiz Eloy Vaz

Obrigado por adquirir estee-book

Esta obra é acompanhada de conteúdo complementar.Para acessá-lo, encaminhe a confirmação de compradeste e-book para pin@elsevier.com.br, solicitando seucódigo de acesso.

Sumário

Capa

Folha de rosto

Obrigado por adquirir este e-book

Cadastro

Copyright

Agradecimentos

Prefácio

Capítulo 1. Introdução

Capítulo 2. Fundamentos matemáticos

2.1 Aproximação De Funções

2.2 Integração Numérica

2.3 Representação Paramétrica De Um Quadrilátero

Capítulo 3. A evolução do método dos deslocamentos

3.1 Método Básico

3.2 Método Clássico

3.3 Método Da Análise Matricial

3.4 Método De Castigliano

3.5 Princípio Dos Deslocamentos Virtuais

3.6 Método Da Mínima Energia Potencial Total

3.7 Método De Rayleigh-Ritz

3.8 O MEF Para Vigas

3.9 O Método Dos Resíduos Ponderados De Galerkin

3.10 Generalização Do MEF

Capítulo 4. Problemas de estado plano

4.1 Introdução

4.3 Elementos Da Família Serendipity

4.4 Elementos Da Família De Lagrange

4.5 Exemplos De Problemas De Estado Plano

Capítulo 5. Sólidos de revolução ou axissimétricos

5.1 Introdução

5.2 Elemento Da Família Serendipity De 4 Nós

5.3 Exemplo De Sólido De Revolução, Placa Circular Vazada

Capítulo 6. Sólidos tridimensionais

6.1 Introdução

6.2 Elemento Tetraedro

6.3 Elemento Hexaedro

6.4 Exemplo De Barra Tracionada Modelada Com SólidoTridimensional, Elemento Hexaedro

Capítulo 7. Placas à flexão

7.1 Introdução

7.2 Teorias De Placa À Flexão

7.3 Elemento Retangular De Placas À Flexão Pela Teoria DeKirchhoff

7.4 Elemento Da Família Serendipity Pela Teoria De Mindlin

7.5 Exemplos De Placa À Flexão

Capítulo 8. Análise de estabilidade

8.1 Introdução

8.2 Obtenção Da Carga Crítica Em Pilares Via Solução DasEquações Diferenciais

8.3 Método Aproximado De Rayleigh-Ritz Para Cálculo DaCarga Crítica Em Pilares

8.4 MEF Para O Cálculo Da Carga Crítica Em Pilares

8.5 MEF Para Cálculo Da Carga Crítica Em Placa À Flexão

8.6 Exemplos De Análise De Estabil idade Por ElementosFinitos

Capítulo 9. Análise dinâmica de estruturas

9.1 Introdução

9.2 Equação De Equilíbrio Em Análise Dinâmica

9.3 Matriz De Massa Do Elemento De Viga

9.4 Matriz De Massa Do Elemento Triangular CST

9.5 Matriz De Massa Do Elemento Serendipity QuadrilateralDe 4 Nós

9.6 Frequências E Modos De Vibração Naturais

9.7 Matrizes De Amortecimento

9.8 Análise Modal De Estruturas Para Vibrações Forçadas

9.9 Análise Dinâmica Por Algoritmo De Integração Direta

9.10 Exemplos De Análise De Vibrações Livres

Capítulo 10. Análise com comportamento não linear domaterial

10.1 Sistema De Equações De Equilíbrio Não Linear

10.2 Solução De Sistemas De Equações Não Lineares

10.3 Exemplo De Aplicação Em Treliça

10.4 Análise Não Linear Detalhada Da Treliça

10.5 Análise Não Linear Alternativa

10.6 Exemplo De Análise Não Linear Da Treliça Com AFormulação Do Item 10.6

Índice de figuras

Referências bibliográficas

Cadastro

Preencha a ficha de cadastro no finaldeste livro e receba gratuitamenteinformações sobre os lançamentos e aspromoções da Editora Campus/Elsevier.Consulte também nosso catálogocompleto, últimos lançamentos emwww.elsevier.com.br

Copyright

© 2011, Elsevier Editora Ltda.

Todos os direitos reservados e protegidos pela Lei nº9.610, de 19/02/1998.Nenhuma parte deste livro, sem autorização prévia porescrito da editora, poderá ser reproduzida outransmitida se jam quais forem os meios empregados:eletrônicos, mecânicos, fotográficos, gravação ouquaisquer outros.

Copidesque e revisão: Globaltec Editora Ltda.Editoração Eletrônica: Globaltec Editora Ltda.

Elsevier Editora Ltda.Conhecimento sem FronteirasRua Sete de Setembro, 111 – 16º andar20050-006 – Centro – Rio de Janeiro – RJ – Brasil

Rua Quintana, 753 – 8º andar04569-011 – Brooklin – São Paulo – SP

Serviço de Atendimento ao Cliente0800-0265340sac@elsevier.com.br

ISBN 978-85-352-3929-4ISBN (versão eletrônica): 978-85-352-6655-9

Nota: Muito zelo e técnica foram empregados naedição desta obra. No entanto, podem ocorrer erros dedigitação, impressão ou dúvida conceitual. Emqualquer das hipóteses, solicitamos a comunicação aonosso Serviço de Atendimento ao Cliente, para quepossamos esclarecer ou encaminhar a questão.Nem a editora nem o autor assumem qualquerresponsabilidade por eventuais danos ou perdas apessoas ou bens, originados do uso desta publicação.

CIP-Brasil. Catalogação-na-fonteSindicato Nacional dos Editores de Livros, RJ

V495m Vaz, Luiz EloyMétodo dos elementos finitos em análisede estruturas /Luiz Eloy Vaz. – Rio de Janeiro: Elsevier,2011.

Inclui bibliografia

ISBN 978-85-352-3929-4

1. Método dos elementos finitos.2. Teoria das estruturas.

10-5751 CDD: 620.0015CDU: 62

Agradecimentos

Aos meus pais, Milton e Alice , pelo amor e carinho. Aeles, o meu reconhecimento pelo exemplo, pela firmeorientação e por não terem poupado esforços para meproporcionar uma boa formação.

À minha esposa Regina, engenheira como eu, queem importantes momentos da minha vida profissionalnão hesitou em sacrificar temporariamente seusestudos e sua carreira para me acompanhar nodoutorado na Alemanha e no pós-doutorado no País deGales. Sem sua generosidade e apoio este livro nãoexistiria. A ela, minha gratidão e amor.

Aos meus mestres da graduação e pós-graduação.Na graduação da UFRJ, mestres como os professoresJosé Luiz Cardoso, Ignacio de Loyola Benedicto Ottonie Benjamin Ernani Dias despertaram meu interessepela análise e pelo projeto de estruturas. No mestradoda Coppe/UFRJ, os professores Fernando Luis LoboCarneiro e Fernando Venâncio Filho aguçaram meuinteresse pela pesquisa. Meu agradecimento especialao professor Venâncio que me iniciou no Método dos

Elementos Finitos e abriu meus olhos para a suaenorme potencialidade. Aos professores José OliveiraPedro, John Argyris e Ernest Hinton, que mereceberam, respectivamente, para um estágio noLaboratório Nacional de Engenharia Civil de Lisboa,para o doutorado na Universidade de Stuttgart e para opós-doutorado na Universidade de Wales em Swansea,minha profunda gratidão. Eles foram fundamentaispara o meu amadurecimento acadêmico. Um especialcarinho eu guardo pelo professor Kaspar Willam, daUniversidade de Stuttgart, pela dedicada orientação eapoio durante a minha tese de doutorado. Hoje , oprofessor Kaspar Willam é professor na Universidadede Boulder, no Colorado.

Aos meus colegas e parceiros em co-orientações eprojetos de pesquisa. Devido à variedade dos temas demeu interesse e por ter trabalhado em trêsimportantes universidades, como a PUC-RJ, a UFRJ e aUFF, e les são numerosos e de perfil diversificado. Nãoposso deixar de citar os professores Eurípedes doAmaral Vargas Jr., Luiz Fernando Martha, Marta deSousa Lima Velasco e Giuseppe Guimarães Barbosa,da PUC-Rio, os professores Sergio Hampshire , ClaudiaEboli e José Herskovits, da UFRJ, a professora SilvanaMaria Bastos Afonso, da UFPE, e , mais recentemente, oProfessor Emil Sanches, da UFF. Eles ajudaram aampliar meus horizontes ao despertar meu interessepor novos temas de pesquisa.

Ao Ivan Menezes, coordenador de projetos doTecgraf PUC-Rio e meu ex-orientando de mestrado.

Sua cuidadosa leitura dos manuscritos e valiosassugestões o tornam praticamente um coautor do livro.

Ao Paul Antezana, pela competente colaboração naedição do texto.

Aos meus alunos de graduação e pós-graduação emeus orientandos de mestrado e doutorado. Elesforam o grande incentivo para meu contínuoaprendizado e crescimento acadêmico. Suas dúvidas equestionamentos me forçaram a compreender osconceitos com mais profundidade e clareza e aprocurar um aperfeiçoamento didático.

A todos os referidos e a muitos outros que não foramcitados, meu sincero “muito obrigado”. Espero queeste livro este ja à altura da valiosa contribuição detodos.

À editora Elsevier, especialmente a André GerhardWolff e Vanessa Vilas Bôas Huguenin, pela confiançadepositada no meu trabalho e pela oportunidade depublicar esta obra.

Prefácio

Este livro surgiu das notas de aulas que preparei paraa disciplina Método dos Elementos Finitos que vem sendoministrada por mim há cerca de 10 anos para alunos degraduação da especialidade de estruturas do curso deEngenharia Civil da Escola Politécnica da UFRJ.

Ao ser indicado para lecionar a disciplina me depareicom a dificuldade de escolher um livro-texto. Osmateriais disponíveis propunham-se a ser umaexcelente fonte de consulta para quem já conhecia ométodo, mas não uma ferramenta para iniciar umaluno de Engenharia que se interessasse pelo tema.Algumas vezes, e les usavam conhecimentosmatemáticos que não eram do domínio dos alunos degraduação — como cálculo variacional — paraapresentar o tema; outras vezes, por serem muitoextensos e detalhados, dificultavam a compreensão daessência do método.

Esta obra tem a intenção de fornecer ao le itor, se jaele um aluno de graduação, de pós-graduação ou umengenheiro em um primeiro contato com o assunto, um

texto compreensível para aqueles que tiveram umaformação básica na área de análise de estruturas. Porformação básica nessa área considero conhecimentosem análise de estruturas hiperestáticas, resistênciados materiais e fundamentos da teoria da elasticidade.Alguns conhecimentos matemáticos que são tratadosnos cursos básicos de Engenharia, mas, em geral, nãocom a profundidade necessária ao estudo do método,como integração numérica, são revistos no início dolivro.

Estou convencido de que a vasta difusão do uso decomputadores nos projetos de Engenharia e a grandedisponibilidade de programas comerciais para análisede estruturas pelo Método dos Elementos Finitos tornam oensino do método nos cursos de graduaçãoindispensável. Este livro pretende ser uma estradamenos sinuosa e íngreme para todos aqueles quepretendam entrar no universo dos elementos finitos.

C A P Í T U L O 1

IntroduçãoO Método dos Elementos Finitos (MEF) para a análisede estruturas ganhou projeção internacional a partir demeados dos anos cinquenta do século XX com ostrabalhos independentes e quase simultâneos doprofessor John Argyris, que trabalhava no ImperialCollege em Londres, e de um grupo de engenheiros daBoeing liderados pelo professor Ray W. Clough.

No entanto, um trabalho sobre o problema de torçãode Saint-Venant do matemático alemão RichardCourant, publicado em 1943, é considerado até hoje opioneiro do método. Na época em que foi publicado,esse trabalho não teve, todavia, grande repercussão.Talvez esse fato possa ser atribuído ao pouco apelodos métodos numéricos em um momento em que aindústria de computadores estava em faseembrionária.

Não se pode, contudo, falar do desenvolvimento e dadivulgação do método sem citar o prof. O. C.Zienkiewicz que trabalhou desde 1961 no campus deSwansea da Universidade do País de Gales, no ReinoUnido. Seu livro publicado em 1967, intitulado “The

Finite Element Methods for Engineering” ficouconhecido no meio acadêmico como “The Book”. Olivro criou uma legião de seguidores do método emtodo o mundo.

No Brasil, a primeira tese sobre o MEF foi defendidana Coppe-UFRJ, em 1970. Ela foi apresentada peloengenheiro Alcebíades Vasconcelos e foi desenvolvidaem parte no Laboratório de Engenharia Civil de Lisboa.Alcebíades desenvolveu um programa para a análisede estruturas de estado plano com o uso do elementotriangular CST, resolveu alguns problemas a cujasolução se chega por meio da Teoria da Elasticidade ecomparou os resultados obtidos pelo programa com osfornecidos pela Teoria da Elasticidade. O primeirocurso sobre o método foi ministrado também naCoppe-UFRJ pelo professor Fernando Venâncio Filhoem 1971.

O MEF foi um desenvolvimento natural daformulação em deslocamentos da análise matricial deestruturas reticuladas impulsionado pelo crescimentodo uso de computadores nas universidades, centros depesquisa e na grande indústria. A semelhança entre osdois métodos consiste no uso comum dos conceitos dematriz de rigidez de elemento, montagem (assembly, eminglês) da matriz de rigidez da estrutura a partir dacontribuição das matrizes de rigidez dos elementos edo conceito de cargas equivalentes nodais. O MEFdistingue-se do seu precursor pela sua maiorgeneralidade e por suas raízes nos métodos de energiae nos métodos aproximados. A análise matricial de

estruturas reticuladas sistematizou o método clássicodos deslocamentos e unificou a metodologia para aanálise de diferentes tipos de estruturas reticuladas,tais como treliças planas e espaciais, vigas e grelhas epórticos planos e espaciais. O MEF, porém, foi bemmais além, ele pode ser usado para se formular tantoproblemas de análise de estruturas reticuladas, comotambém de estruturas contínuas bi e tridimensionais.Sua generalidade não parou por aí, sua aplicação, quese iniciou em análise estática de estruturas decomportamento linear elástico, foi estendida à análiseestática de estruturas com não linearidade física egeométrica e à análise dinâmica de estruturas. Eletambém saiu da esfera da análise de estruturas epenetrou em outras áreas, como a engenhariageotécnica, a interação fluido-mecânica e as análisesde fluxo térmico e hidráulico.

Na área de análise de estruturas, a formulação doMEF pode ser feita a partir do Princípio da MínimaEnergia Potencial Total, do Método de ResíduosPonderados ou do Princípio dos DeslocamentosVirtuais. Ele usa os conceitos de “discretização” docontínuo e de “matriz de interpolação” que fornece osdeslocamentos em um ponto no interior do elementoem função de seus deslocamentos nodais. O termodiscretização se refere a um modelo com um númerofinito (discrete, em inglês) de incógnitas (deslocamentosnos nós do modelo) para a análise de meios contínuosem contraposição a uma análise com um númeroinfinito de variáveis como as feitas pela Teoria da

Elasticidade que usam funções contínuas, ou se ja, cominfinitas incógnitas como solução.

Hoje em dia, existem inúmeros programascomerciais altamente sofisticados que fazem os maisdiversos tipos de análise pelo Método dos ElementosFinitos, tais como o SAP, o Ansys, o Abaqus, o Nastranetc. No Departamento de Mecânica Aplicada eEstruturas da Escola Politécnica da UFRJ, está emdesenvolvimento o sistema Salt sob a coordenação doprofessor Silvio de Souza Lima. O programa tem sidolargamente utilizado na elaboração de diversostrabalhos de fim de curso de alunos do departamento.No Tecgraf, na PUC-Rio, há o sistema Mtool comgerador automático de malhas.

Em minha opinião, a difusão do uso do MEF nasempresas e universidades tornou obrigatória aintrodução de um curso sobre o método nas disciplinasde graduação em engenharias civil, mecânica, naval eaeronáutica.

Este livro tem como objetivo servir de base para adisciplina “Introdução ao Método dos ElementosFinitos” que seria ministrada em um curso degraduação em Engenharia Civil na ênfase deEstruturas. O material é adequado para um curso de 16semanas com 3 horas semanais.

O Capítulo 2 faz uma revisão aprofundada de algunsfundamentos matemáticos já vistos no ciclo básico deEngenharia necessários ao longo do curso, comointegração numérica.

O Capítulo 3 mostra a evolução do Método dos

Deslocamentos, desde as formulações clássicas paraestruturas reticuladas até o MEF, visto como umaevolução do Método de Rayleigh-Ritz.

O Capítulo 4 trata das formulações do método para aanálise de estruturas planas, apresentando asformulações do elemento CST, de elementos dasfamílias Serendipity e de Lagrange.

O Capítulo 5 apresenta formulações do método paraanálise de sólidos axissimétricos ou sólidos derevolução, mostrando as formulações de algunselementos, como o Triangular de três nós e elementosda família Serendipity.

O Capítulo 6 aborda formulações do método paraanálise de sólidos tridimensionais, desenvolvendo asformulações de alguns elementos, como o elementotetraedro e o hexaedro.

No Capítulo 7, são estudados elementos para aanálise de placas à flexão, como o elemento retangular,baseado na Teoria de Kirchhoff, próprio para a análisede placas delgadas e os elementos da famíliaSerendipity, baseados na Teoria de Mindlin eapropriados à análise de placas espessas.

O Capítulo 8 trata do problema do cálculo do fator decarga crítica em estruturas. Formulações da matriz derigidez geométrica são apresentadas para estruturasde pórticos planos e de placas, assim como exemplosnuméricos.

O Capítulo 9 contempla o estudo de análise dinâmicaem estruturas. É apresentada a formulação para seobter as frequências e os modos próprios de estruturas

em vibrações livres a partir da matriz de rigidez e damatriz de massa consistente para alguns elementosfinitos. A obtenção da matriz de amortecimentotambém é tratada. Finalmente, são estudadas a análisemodal e a análise por algoritmo de integração direta deNewmark de estruturas submetidas a vibraçõesforçadas. Exemplos referentes a todos os itens sãoapresentados.

O Capítulo 10 aborda a análise de estruturas comcomportamento não linear do material. O conceito dematriz de rigidez tangente é apresentado e umexemplo é resolvido com o uso do Método de Newton-Raphson.

Espero com esse texto facilitar o aprendizado desseapaixonante e revolucionário tema que é o Método dosElementos Finitos.

Prof. Luiz Eloy Vaz

Professor titular em Análise de Estruturas pela UFRJaté 2008

Professor adjunto da UFF a partir de 2009

Demo version limitation, this page not show up.

C A P Í T U L O 3

A evolução do método dosdeslocamentosO Método dos Elementos Finitos (MEF) tratado nestelivro pertence à família do Método dos Deslocamentosou Método da Rigidez onde deslocamentos sãoescolhidos como incógnitas. Todos os membros dessafamília se caracterizam por ter como equaçãofundamental a equação de equilíbrio cujas incógnitassão deslocamentos generalizados. Entendem-se aquipor deslocamentos generalizados, grandezascinemáticas, tais como, deslocamentos lineares,rotações etc.

Os membros dessa família formam uma árvoregenealógica, com novos métodos gerados a partir dosmétodos mais antigos. De certa maneira, a evoluçãodo método ao longo do tempo segue as le is daevolução de Darwin, com mutação e seleção. Os novosmembros da família desses métodos herdam ascaracterísticas de seus antecessores, mas sofrempequenas mudanças que só são bem sucedidas seforem bem adaptadas às condições existentes. Um

exemplo disso é que a Análise Matricial de Estruturas(AME) e o MEF só tiveram larga aceitação quando oscomputadores atingiram uma fase de elevado grau dedesenvolvimento, apesar de este último ter surgidoantes dessa fase.

Este capítulo procura mostrar como se deu aevolução do Método dos Deslocamentos, desde asprimeiras formulações até o MEF. É surpreendenteverificar como as mudanças conceituais são pequenasem comparação ao enorme crescimento do potencialdo método.

3.1 Método básicoA análise de estruturas usa três equações básicas,nomeadamente equações de compatibilidade, deequilíbrio e constitutivas, também chamadas derelação tensão-deformação. O método dosdeslocamentos caracteriza-se por usar a equação deequilíbrio como equação fundamental, ou se ja, aquelade onde são obtidas as incógnitas primárias doproblema, a partir das quais, todas as outras respostasserão obtidas. As incógnitas primárias são osdeslocamentos por meio dos quais é possível obterdeformações, tensões, resultantes de tensões etc.

O método básico da família do método dosdeslocamentos consiste em manipular as trêsequações básicas da análise de estruturas de modo acolocar todas as informações disponíveis nas equações

de equilíbrio com deslocamentos livres comoincógnitas. O número de deslocamentos livres échamado de grau de liberdade da estrutura.

Neste item e em outros que seguem, a estruturaapresentada na Figura 3.1 é utilizada para ilustrar aresolução do método. Trata-se de uma treliça planasimples com quatro barras e dois graus de liberdade,os deslocamentos horizontal e vertical do nó C.

FIGURA 3.1 Treliça com 2 graus de liberdade.

As equações de compatibilidade relacionam

(3.1)

grandezas cinemáticas, nesse caso os deslocamentosnodais livres d1 e d2 na direção horizontal e vertical comalongamentos/encurtamentos δi das barras i. Osdeslocamentos são supostos positivos com os sentidosindicados na Figura 3.1. Os alongamentos serãoconsiderados positivos e os encurtamentos negativos.As expressões para osδi das quatro barras são obtidasprojetando-se os deslocamentos nodais nas direçõesdas barras, assim:

A segunda equação de compatibilidade relaciona osalongamentos/encurtamentos das barrasδi com asdeformações longitudinais εi. Da resistência dosmateriais:

(3.2)

(3.3)

Como os comprimentos das barras são:

Chega-se a:

(3.6)

(3.4)

(3.5)

Para efeito de simplificação, a le i constitutivausada nesse trabalho será a le i de Hooke, Assim, paracada barra, i vale :

Ou, em termos de esforços normais Ni,

Onde E é o módulo de elasticidade do material, A, aárea de seção transversal (as duas grandezas

(3.7)

supostas constantes para todas as barras), Ni o esforçonormal e Li o comprimento da barra i.

Substituindo-se para cada barra i, δi dado em (3.1) em(3.6), obtém-se:

As equações de equilíbrio são obtidas para asdireções horizontal e vertical no nó C. Os sentidos dasforças axiais Ni que atuam nas barras i, são admitidos aprincípio como de tração. Para se escrever asequações de equilíbrio, valem, no entanto, os sentidosindicados na Figura 3.2.

(3.8)

FIGURA 3.2 Equilíbrio do nó C.

As equações de equilíbrio são: Na direção horizontal:

Na direção vertical:

(3.9)

(3.10)

(3.11)

Substituindo-se as expressões (3.7) em (3.9) emanipulando-as, obtém-se:

A expressão (3.10) é a equaçãofundamental do método dos deslocamentos para aanálise da treliça plana da Figura 3.1. Matricialmente,ela pode ser reescrita como:

Cuja solução é :

(3.12)

(3.13)

(3.14)

Com os deslocamentos d1 e d2 é possível obter agoratodas as respostas da estrutura em termos dealongamento/encurtamento, na expressão (3.1),deformações, em (3.4), tensões, em (3.5), e esforçosnormais Ni, em (3.7). Tais valores estão indicados aseguir:

(3.15)

(3.16)

3.2 Método clássicoO método clássico é essencialmente o mesmo que ométodo básico. Sua contribuição foi no sentido desistematizar, ou se ja, organizar, ou ainda criar umametodologia que possa ser aplicada da mesma forma atodas as estruturas.

O método usa os conceitos de estados auxiliares ede superposição de efeitos. Inicialmente, devem-seidentificar os graus de liberdade da estrutura. Emseguida, um estado auxiliar j é criado para cada graude liberdade impondo-se um valor unitário para o graude liberdade dj, enquanto os outros são mantidos nulos.

Resultantes das forças internas resistentes que atuamnas barras aparecem nas direções dos graus deliberdade. A força interna na direção i devido aodeslocamento unitário na direção do grau de liberdadedj é chamada de coeficiente de rigidez k ij. Além disso,um estado auxiliar 0 é criado para as cargas atuantescom todos os graus de liberdade mantidos fixos. Asforças resultantes que atuam nos nós na direção dograu de liberdade dj nesse estado são denominadascargas nodais fj.

Como os estados auxiliares não são autoequilibradoso equilíbrio é conseguido com a superposição deefeitos. Assim, somando-se os produtos das forçasinternas resultantes (nas direções dos graus deliberdade) correspondentes a cada estado auxiliar j pordj, a soma deve ser igual às forças aplicadas (nasdireções dos graus de liberdade) no estado auxiliar 0.Em termos físicos, isso significa que os deslocamentosque surgem na direção dos graus de liberdade djdevem ser tais que as forças internas equilibrem asforças aplicadas.

A aplicação das ideias descritas no exemplo do item3.1 ajuda a esclarecer o método.

Estado auxiliar 1, d1 = 1.

FIGURA 3.3 Termos k11 e k21 da matriz de rigidez datreliça.

Estado auxiliar 2, d2 = 1.

FIGURA 3.4 Termos k21 e k22 da matriz de rigidez datreliça.

Para se obter os coeficientes k ij (força internaresultante na direção i devida a um deslocamentounitário na direção j) procede-se da seguinte maneira:inicialmente, calculam-se osalongamentos/encurtamentos das barras dij(alongamento/encurtamento na barra i devido a umdeslocamento unitário na direção do grau de liberdadedj) de forma análoga ao que foi fe ito para se obter os

(3.17)

(3.18)

alongamentos/encurtamentos em (3.1). Para o estado auxiliar 1.

Para o estado auxiliar 2.

(3.19)

(3.20)

Utilizando-se a relação constitutiva é possívelcalcular os esforços normais nas barras Nij (esforçonormal na barra i devido a um deslocamento unitáriona direção do grau de liberdade dj) com uma expressãoanáloga a (3.6).

Assim: Para o estado auxiliar 1.

Para o estado auxiliar 2.

(3.21)

(3.22)

Os coeficientes de rigidez k ij (esforço na direção ipara um deslocamento unitário na direção j) sãocalculados utilizando-se as equações de equilíbrio nonó C. Assim, das equações de equilíbrio na direçãohorizontal e vertical da Figura 3.5, da correspondente ad1 = 1 obtém-se, respectivamente, os coeficientes k11 ek21.

Para o estado auxiliar 1, Figura 3.5a.

Para o estado auxiliar 2, Figura 3.5b.

(3.23)

(3.24)

O estado auxiliar 0, fornece:

FIGURA 3.5 Forças no nó C para d1 = 1 e d2 = 1.

A superposição de efeitos, que deve garantir oequilíbrio das forças resistentes e aplicadas, podeagora ser escrita como:

(3.25)

(3.26)

ou com os valores da estrutura sendo analisada:

A expressão (3.26) é idêntica à expressão(3.11), como não poderia deixar de ser. Desse modo, asrespostas das estruturas obtidas pelo método básicodadas pelas expressões de (3.12) a (3.16) serão asmesmas.

3.3 Método da análise matricial3.3.1 Formulação Da AnáliseMatricialA análise matricial de estruturas reticuladassistematizou as operações matemáticas da análise deestruturas fazendo uso da álgebra matricial que operacom vetores e matrizes. Ela introduziu diversosconceitos novos na análise de estruturas. Toda asistematização se baseia na ideia de sistema local esistema global de coordenadas. Com esse conceitodefinido, é possível estabelecer matrizes de rigidez de

elemento nos sistemas local e global, assim comovetores de forças nodais de elemento nos sistemaslocal e global. A partir das contribuições das matrizesde rigidez e dos vetores de forças nodais de elementono sistema global, pode-se montar a matriz de rigidezbem como o vetor de forças nodais da estrutura.Deslocamentos nodais também são definidos nossistemas local e global. Uma equação de equilíbrio daestrutura no sistema global fornece os deslocamentosnodais. Uma vez obtidos os deslocamentos nodais daestrutura, as forças atuantes nas extremidades doselementos podem ser determinadas.

O sistema local de coordenada é definido quando seescolhe os nós inicial e final do elemento. Na Figura 3.6,os nós 1 e 2 são, respectivamente, o nó inicial e o nófinal do elemento ou barra. O eixo x local fica entãodefinido na direção da barra e com sentido positivo de1 para 2. O eixo y é perpendicular a x, com o vetor dosentido positivo fazendo 90 graus a partir de x nosentido anti-horário. O sistema global dado pelos eixosX e Y é definido usualmente da seguinte maneira: Xtem direção horizontal e sentido positivo da esquerdapara a direita, e Y tem direção vertical e sentidopositivo de baixo para cima. O sistema global não éobrigatoriamente o definido anteriormente, podendoser escolhido outro que se ja mais conveniente.

FIGURA 3.6 Graus de liberdade no sistema global elocal.

A estrutura de treliça plana tratada até aqui tem doisgraus de liberdade por nó. Ao nó 1 são associados osdeslocamentos 1 e 2 e ao nó 2, os deslocamentos 3 e 4.A Figura 3.6 indica os sentidos positivos dos 4componentes do vetor de deslocamentos dl, no sistemalocal, e dg, no sistema global. O ângulo a define arotação do eixo da barra em relação ao sistema global.Associados aos vetores de deslocamentos, são criadostambém os vetores de forças nodais fl, no sistemalocal, e fg, no sistema global.

Os vetores dos deslocamentos de elemento no

(3.27)

(3.28)

(3.29)

(3.30)

sistema local dl e global dg podem ser relacionados pelamatriz de rotação R, como indicado a seguir:

Ou, sucintamente:

Como o trabalho é um escalar que independe dosistema de coordenadas, e le deve ser o mesmo nossistemas local e global.

Substituindo (3.28) em (3.30), obtém-se:

(3.33)

(3.31)

(3.32)As expressões (3.28) e ((3.32)) formam o princípio da

contragradiência que pode ser enunciado como: “Seuma matriz transforma deslocamentos globais emlocais, sua transposta transforma forças locais emglobais.”

A matriz de rigidez do elemento de treliça plana nosistema local para o elemento m, Kl,m é dada em (3.33).Ela é obtida da definição dos coeficientes de rigidezk l,m(ij). O coeficiente k l,m(ij) significa a força na direção dodeslocamento local i para um deslocamento unitárioaplicado na direção do deslocamento local j, mantendoos outros deslocamentos locais nulos.

Onde Em é o módulo de elasticidade do material, Am aárea da seção transversal e Lm o comprimento da barra

(3.34)

(3.35)

(3.36)

m. A equação de equilíbrio da barra que relacionadeslocamentos, forças e a matriz de rigidez no sistemalocal de coordenadas é dada por:

Ou, sucintamente:

A matriz de rigidez do elemento m no sistema globalde coordenadas pode ser obtida como explicado aseguir. Substituindo-se (3.28) em ((3.35)), obtém-se:

Pré-multiplicando-se ambos os lados de (3.36) por , chega-se a:

(3.37)

(3.38)

(3.39)

(3.40)

Usando (3.32), obtém-se:

Onde,

A partir da matriz de rigidez e das forças nodais decada elemento k no sistema global é fe ita então amontagem da matriz de rigidez K e das forças nodais fglobais da estrutura em função da conexão entre oselementos (incidência), obtendo-se a equação deequilíbrio global da estrutura.

Sendo d os deslocamentos da estrutura no sistemaglobal de coordenadas.

Uma vez obtido d, é possível calcular osdeslocamentos nodais de cada elemento no sistemaglobal e girar esses deslocamentos para o sistemalocal via (3.28) e calcular as forças de extremidade

finais em cada elemento no sistema local via (3.35).

3.3.2 Aplicação Da Análise MatricialA aplicação das ideias descritas no exemplo do item 3.1ajuda a esclarecer o método.

A treliça plana estudada nesse item é reproduzidamais uma vez na Figura 3.7.

FIGURA 3.7 Treliça plana com 2 graus de liberdade.

O sentido positivo do eixo local x das barras édefinido como: Barra 1: do nó A para o nó C; Barra 2: donó B para o nó C; Barra 3: do nó D para o nó C; Barra 4:do nó E para o nó C (Figura 3.8).

FIGURA 3.8 Sistemas de coordenadas locais dasbarras.

Os comprimentos Lm da barra m são: ,

e .As matrizes Rm das quatro barras são:

(3.41)Usando a expressão (3.39) para se obter asmatrizes de rigidez de elemento no sistema global e

(3.42)

somando apenas os termos referentes às duas últimaslinhas e colunas de cada matriz (isso se explica porqueos nós iniciais de todas as barras estão vinculados e ,portanto seus deslocamentos são nulos), obtém-se amatriz de rigidez da estrutura no sistema global Krelativa aos dois graus de liberdade do nó C, dada em(3.42). A equação de equilíbrio da estrutura é a mesmajá obtida em (3.26), o que conduz aos mesmosresultados.

3.4 Método de CastiglianoO Método de Castigliano é assim chamdo emhomenagem ao segundo teorema de Carlo AlbertoCastigliano, que, em 1873, demonstrou que a derivadada energia de deformação de uma estrutura emrelação ao deslocamento di é igual a força externa daestrutura na mesma direção. A demonstração foi fe ita

para estruturas com comportamento linear elástico,mas ela é válida também para materiais e lásticos nãolineares. Nesse item, a demonstração será estendida aestruturas de material e lástico não linear.

Esse teorema representou um importante passo nodesenvolvimento da análise de estruturas porque elemostrou um novo caminho, baseado em teoremas deenergia, para se formular um método para análise deestruturas. Esse caminho levou ao MEF.

3.4.1 Energia De DeformaçãoPara efeito de simplificação, a apresentação doSegundo Teorema de Castigliano será feita aqui para ocaso particular de uma estrutura de treliça. Nesse tipode estrutura, somente uma componente dedeformação e de tensão atua no elemento de barra,nomeadamente, a deformação e a tensão normallongitudinal, ou se ja, trata-se de um problemaunidimensional para efeito da relação tensão xdeformação. Se ja a relação tensão x deformaçãoapresentada na Figura 3.9. A solicitação externa levoua tensão atuante até o valor final que correspondeà deformação final na barra m da treliça.

(3.43)

FIGURA 3.9 Energia de deformação específica U0 dabarra m.

A energia de deformação específica na barra mé definida como:

O adjetivo “específica” deve-se ao fato de ser, emtermos de unidades, trabalho por unidade devolume.

(3.44)

(3.45)

(3.46)

A energia de deformação da barra m, um, é obtidaintegrando-se no volume da barra.

Para se obter a energia de deformação U re lativa atoda a treliça, somam-se os Um de todas as barras, de 1a nb, onde nb é o número de barras da estrutura.

Onde é a deformação final da barra m. Como adeformação final da barra, depende doalongamento/encurtamento longitudinal final da barra

, como expresso em (3.2), que, por sua vez, depende dos deslocamentos nodais finais dasextremidades da barra no sistema global decoordenadas como exemplificado em (3.4), aexpressão (3.45) pode ser reescrita como:

(3.47)

Onde n é o número de graus de liberdade daestrutura de treliça.

A energia de deformação da estrutura correspondefisicamente à energia armazenada na estruturaquando ela se deforma, caso não haja perda deenergia, ou se ja, para um sistema conservativo. Essaenergia é responsável pela volta da estrutura a suaconfiguração inicial, antes da aplicação das cargas,quando estas são retiradas da estrutura.

3.4.2 Trabalho ExternoO trabalho externo total W em uma estrutura de treliçaplana pode ser obtido somando-se os trabalhosexternos W i referentes aos graus de liberdade i daestrutura.

Onde n, como anteriormente, é o númerode graus de liberdade da estrutura. A Figura 3.10esclarece.

FIGURA 3.10 Trabalho externo associado ao grau deliberdade i.

3.4.3 Segundo Teorema De CastiglianoSubstituindo doravante a notação do deslocamentofinal por d para efeito de simplificação, a energia dedeformação (3.46) e o trabalho externo (3.47) em umaestrutura de treliça plana, como visto nos itens 3.4.1 e3.4.2, podem ser escritos como uma função do vetordos deslocamentos nodais finais da estrutura nosistema global de coordenadas d com n componentes.

Expandindo-se W(d) em série de Taylor até o termode primeira ordem, é possível expressar o incremento

(3.48)

(3.50)

(3.49)

(3.51)

de W(d) como:

Procedendo-se da mesma maneira paraU(d), obtém-se:

Pelo princípio da conservação de energiaem sistemas conservativos, todo trabalho externorealizado é armazenado na estrutura em termos deenergia de deformação. Assim, o incremento detrabalho externo é igual ao incremento de energia de

(3.53)

(3.54)

(3.55)

(3.52)

deformação, logo:

Ou se ja,

Ou, ainda, para uma variação arbitrária δd,

O teorema da integral de Newton diz que:

Logo, utilizando-se esse teorema, pode-se escrever:

(3.57)

(3.58)

(3.56)Onde, como foi redefinido no início desseitem, di em (3.56) é o valor final da variáveldeslocamento nodal ui e fi é a força final associada aodeslocamento di.

Com o uso de (3.54) e (3.56), obtém-se finalmente aexpressão do Segundo Teorema de Castigliano:

Ou, grupando-se todas as equações (3.59)correspondentes aos n graus de liberdade em uma sóequação:

Observa-se que o termo à esquerda da expressão(3.58) corresponde ao vetor das forças internasresistentes, doravante denominado fr(d), e o termo àdireita, corresponde ao vetor das forças solicitantes,doravante denominado fs.

(3.59)

(3.60)

A expressão (3.59) fornece um método de análise deestruturas denominado Método de Castigliano. Aexpressão fornece n equações que permitem obter as nincógnitas do problema, ou se ja, os n deslocamentosnodais di, i = 1, ..., n. Se a estrutura tiver umcomportamento linear, as equações (3.59) fornecem umsistema de n equações algébricas lineares, caso ocomportamento se ja não linear, n equações nãolineares são obtidas. O sistema de n equações nãolineares pode ser resolvido, por exemplo, pelo métodode Newton-Raphson para se obter as n incógnitas doproblema, ou se ja, os n deslocamentos nodais di, i = 1,..., n.

A aplicação do método na análise da treliça plana daFigura 3.1 ajuda a esclarecer as expressões descritasanteriormente.

3.4.4 Aplicação Do Método DeCastiglianoA lei de Hooke para materiais linear-elásticos permiteescrever:

A energia de deformação específica U0 pode ser

(3.61)

(3.62)

escrita em função da deformação final da barra m.Empregando-se novamente a notação pararepresentar o valor final da grandeza εm, chega-se a:

A energia de deformação Um para a barra mvale :

Usando as equações de compatibilidadepara a treliça da Figura 3.1 descritas em (3.1) eabandonando mais uma vez, para efeito desimplificação, o sobrescrito – para representar valoresfinais das variáveis, obtém-se:

(3.63)

(3.64)

(3.65)

E as expressões dos comprimentos das barras dadasem (3.3), podem-se escrever:

(3.66)

(3.67)

(3.68)

(3.69)

(3.70)

Usando-se (3.46) para se obter a energia dedeformação total da estrutura, obtém-se:

Aplicando-se agora a expressão (3.57) doSegundo Teorema de Castigliano, obtém-se:

Ou, ainda,

(3.71)Que é idêntica a (3.26).

3.5 Princípio dos deslocamentosvirtuais3.5.1 Incrementos Da Energia DeDeformaçãoO princípio dos trabalhos virtuais será demonstradoneste item para estruturas de treliça. Uma barra detreliça m é carregada até que a deformação final se ja atingida como indicado na Figura 3.11. A tensãoatuante correspondente é ( ). A energia dedeformação específica produzida na barra é .Imagine agora que um incremento de tensão δσm se jaaplicado à barra a partir de . Um incremento dedeformação δεm correspondente ocorre na barra.

(3.72)

FIGURA 3.11 Incremento de energia de deformaçãoespecífica DU0,m da barra m.

O incremento total da energia de deformaçãoespecífica Δ correspondente à aplicação de δσmpode ser escrito como:

ou,

(3.75)

(3.76)

(3.73)

(3.74)

onde

Os termos e são denominadosincremento de primeira e de segunda ordem de ,respectivamente. O termo de primeira ordemcorresponde à área do retângulo vertical hachuradorepresentado na Figura 3.11. O termo de segundaordem corresponde à área do triângulo maior namesma figura. A área em cinza corresponde ao errocometido no cálculo do incremento total erro .

Como a energia de deformação da barra m da treliçaUm é obtida pela integração no volume da barra daenergia de deformação específica, obtém-se:

(3.79)

(3.80)

(3.77)

(3.78)

Logo,

ou

onde

A energia de deformação de toda estrutura com mbarras pode ser obtida somando-se a energia dedeformação de todas as barras, assim:

(3.84)

(3.85)

(3.81)

(3.82)

(3.83)

Logo,

ou

onde

As expressões (3.84) e (3.85) podem sergeneralizadas para o caso em que há váriascomponentes de tensão, por exemplo, σx, σy e τxy, e dedeformação, por exemplo, εx, εy e γxy atuando em um

(3.86)

(3.87)

elemento infinitesimal do elemento m da estrutura comn e lementos. Nesse caso pode-se escrever:

Onde σm, δσm e δεm representam, respectivamente, osvetores das componentes de tensão atuantes, dosincrementos das componentes de tensão atuantes edos incrementos das componentes de deformação noelemento m.

3.5.2 Incrementos Do Trabalho ExternoOs incrementos do trabalho externo podem ser obtidospelo raciocínio análogo ao desenvolvido no itemanterior para a energia de deformação.

Uma força é aplicada em um dado grau de liberdadei até produzir um deslocamento final comorepresentado na Figura 3.12. A força atuantecorrespondente à é . O trabalho externo

produzido correspondente ao grau de liberdade i é W i.Imagine agora que um incremento de força δfi é

aplicado à força . Um incremento de deslocamentoδdi ocorre no grau de liberdade correspondente.

FIGURA 3.12 Incremento de trabalho externo DWi.

O incremento total do trabalho externo ΔW icorrespondente à aplicação de δfi no grau de liberdadei pode ser escrito como:

(3.89)

(3.91)

(3.88)

(3.90)

ou

onde

Os termos e são denominadosrespectivamente incremento de primeira e de segundaordem de W i. O termo de primeira ordem correspondeà área do retângulo vertical hachurado representadona Figura 3.12. O termo de segunda ordemcorresponde à área do triângulo maior na mesmafigura. A área em cinza corresponde ao erro cometidono cálculo do incremento total erroW i.

O trabalho externo correspondente a toda a treliça

(3.96)

(3.92)

(3.93)

(3.94)

(3.95)

com n graus de liberdade pode ser obtido somando-seo trabalho externo de todos os graus de liberdade,assim:

Logo,

ou, ainda,

onde,

As expressões (3.95) e (3.96) podem ser escritas

(3.98)

(3.97)

usando-se vetores:

Onde , δd e δf representam, respectivamente, osvetores das forças solicitantes nodais finais, dosincrementos dos deslocamentos nodais e dosincrementos das forças nodais.

3.5.3 Formulação Do Princípio DosDeslocamentos VirtuaisO princípio dos deslocamentos virtuais baseia-se noprincípio de conservação de energia. Seu enunciado éo seguinte: “Para toda estrutura, o incremento deprimeira ordem da energia de deformação é igual aoincremento de primeira ordem do trabalho externo.” Aaplicação do princípio não se limita a sistemasconservativos. Matematicamente, e le pode serexpresso por:

(3.100)

(3.99)Para o caso geral em que há várias componentes de

tensão e deformação atuando em um elementoinfinitesimal de um elemento m de uma estrutura com nelementos, a expressão (3.99) pode ser escrita como:

As grandezas δεm e δd em (3.100) são cinemáticas,virtuais e compatíveis enquanto que as grandezas e são ditas estáticas, reais e em equilíbrio. O termovirtual é sinônimo de potencial, ou se ja, pode vir aacontecer, não real. As grandezas δεm e δd estãorelacionadas por equações de compatibilidade já queas componentes de δd produzem as componentes deδεm. As grandezas reais e estão relacionadas porequações de equilíbrio já que as tensões reais sãoproduzidas pelas forças reais .

3.5.4 Exemplo Da Aplicação Do PrincípioDos Deslocamentos VirtuaisInicialmente serão deduzidas as equações decompatibilidade entre as deformações virtuais δεm das

(3.101)

barras m e os deslocamentos virtuais nodais δdi dosgraus de liberdade i.

As expressões são análogas às expressões (3.4),substituindo-se as grandezas reais por grandezasvirtuais.

As tensões reais são expressas emfunção dos deslocamentos reais. Elas podem serobtidas por meio de novas expressões (3.4)multiplicadas pelo modo de elasticidade E paratransformar deformação em tensão pela le i de Hooke.

(3.102)

(3.103)

Substituindo (3.101) e (3.102) na expressão(3.100) e integrando-se no volume de cada barra, ouseja, multiplicando-se por A Lm, pois as tensões sãoconstantes no volume de cada barra m, e considerandoque o termo à direita em (3.100) vale P δd1, chega-se a:

Como δd1 e δd2 são arbitrários, deve-seter:

(3.105)

(3.104)

(3.106)

Ou, matricialmente,

que é , de novo, a mesma expressão (3.26)que conduz aos mesmos resultados anteriores emtermos de deslocamentos nodais di nos graus deliberdade i e de mesmos alongamentos/encurtamentosδm, deformações εm, tensões σm e esforços normais Nmnas barras m conforme obtido no item 3.1.

3.6 Método da mínima energiapotencial total3.6.1 Energia Potencial Total

(3.107)

(3.108)

A energia potencial total Π(d) é definida para sistemasconservativos como:

Onde U(d) é a energia de deformação da estrutura,como definido em (3.44) e (3.46), e Wp(d) é o trabalhopotencial das forças externas, dado por:

Novamente, os sobrescritos –, utilizados pararepresentar valores finais das variáveis são retiradospara efeito de simplificação. Em sistemasconservativos, U(d) é a energia que traz a estrutura devolta à configuração inicial caso as forças externasse jam retiradas da estrutura. Wp(d) é o trabalhopotencial, ou se ja, aquele que seria realizado caso aestrutura voltasse a sua configuração inicial e ascargas permanecessem atuando sobre ela. Assim, Π(d)é a energia total necessária para trazer de volta aestrutura a sua configuração inicial com as cargasatuando sobre ela.

3.6.2 O Princípio Da Mínima EnergiaPotencial Total

(3.109)

(3.110)

O princípio da mínima energia potencial total enunciaque os deslocamentos d de uma estrutura em equilíbrioestável tornam mínima a energia potencial total daestrutura. Em outras palavras, uma estrutura que estáem equilíbrio estável se deformou de modo a gastar omínimo de energia potencial total. Matematicamente, acondição de primeira ordem de mínimo de uma funçãoé dada por:

Ao combinar as expressões (3.58), (3.59), (3.107) e(3.109) pode-se escrever:

Observe que a expressão (3.110) éidêntica à expressão (3.59). Isso significa que osdeslocamentos da estrutura em equilíbrio estável dsatisfazem a equação de equilíbrio (3.59) e minimizama energia potencial total. Quando se usa a expressão(3.110) para obter os deslocamentos d da estrutura, diz-se que a estrutura foi calculada pelo método da mínimaenergia potencial total.

(3.111)

(3.112)

3.6.3 Aplicação Do Princípio Da MínimaEnergia Potencial TotalA energia de deformação total da estrutura foi obtidano item 3.4.4, (vide expressão (3.68)), ou se ja:

A energia potencial total é dada por:

Aplicando o princípio da mínima energia potencialtotal, obtém-se:

ou

(3.113)As expressões (3.112) e (3.113) são,respectivamente, idênticas às expressões (3.67) e (3.68)e conduzem à mesma solução em termos dedeslocamentos d1 e d2, bem como dealongamentos/encurtamentos, deformações, tensões eesforços normais que dependem de d1 e d2.

3.7 Método de Rayleigh-RitzO método de Rayleigh-Ritz representou um grandepasso na evolução do método dos deslocamentos, poiscontribuiu decisivamente para o aparecimento do MEF.O método de Rayleigh-Ritz é , na essência, o método doprincípio da mínima energia potencial total, mas, apequena modificação introduzida nesse último permitiuum grande avanço. Para uma melhor compreensão dométodo, o exemplo da treliça usado até aqui vai sersubstituído por um novo exemplo de análise de umaviga em balanço representada na Figura 3.13.

(3.114)

FIGURA 3.13 Viga em balanço de inércia variável.

Para fazer a análise da viga da Figura 3.13 pelométodo do princípio da mínima energia potencial totalé preciso, inicialmente, obter a expressão para aenergia de deformação de uma viga. A viga,supostamente, deve satisfazer a hipótese de Bernoulli(1705), a qual considera que “seções transversais retaspermanecem planas e normais à tangente ao eixofletido da viga”. O deslocamento vertical do eixo daviga ao longo do comprimento é descrito pela funçãoν(x). Da resistência dos materiais, sabe-se que adeformação longitudinal ε(x,y) no ponto da seção x ecota y é dada por:

Sendo,

(3.115)

(3.117)

(3.116)

(3.118)

A energia de deformação específica de um materiallinear elástico com módulo de elasticidade E, é dadapor:

Para um ponto da seção x e cota y da vigaà flexão:

A energia de deformação da viga pode ser obtidapor:

ou

(3.119)

(3.120)

(3.121)

onde L é o comprimento da viga e I o momento deinércia da seção da viga, dado por:

Como no exemplo em estudo, a inércia da seçãovaria ao longo do comprimento, a energia potencialtotal da viga pode ser obtida por:

Observando a expressão (3.121), verifica-se que a energia potencial total da viga Π é função dafunção que descreve a deformação do eixo da vigaν(x), ainda desconhecida. Uma função de função édenominada um funcional. Esse problema difereradicalmente do problema resolvido no item 3.6.3,quando a estrutura a ser resolvida era uma treliça e Π,dado em (3.111), era uma função das variáveis d1 e d2.

Do ponto de vista matemático o problema anteriorda treliça era um problema de minimização de uma

(3.122)

função de duas variáveis. O problema da viga é umproblema de minimização de um funcional da funçãoν(x). Trata-se agora de encontrar a função ν(x) e nãomais apenas as variáveis d1 e d2 que minimizam Π.Esse é um problema clássico de cálculo variacional, esua solução está fora do escopo deste livro.

Como então resolver o problema da viga à flexão? Éaqui que surge a ideia básica do método de Rayleigh-Ritz: a função ν(x) que representa a elástica da viga édescrita por uma função aproximadora.

As funções aproximadoras devem satisfazer asseguintes condições:

a) Devem ser funções polinomiais outrigonométricas que satisfaçam às condições decontorno em deslocamento da viga.

b) Devem ter derivadas contínuas até a ordem n-1,sendo n a maior ordem de derivação da funçãono funcional Π (no caso n = 2).

c) Devem ser definidas em todo o domínio doproblema.

A solução “exata” para o deslocamento d naextremidade livre da viga da Figura 3.13 é 1875.

Primeira tentativa:A primeira função aproximadora adotada é um

polinômio de segundo grau.

(3.123)

(3.124)

(3.125)

(3.126)

Vale observar que a função satisfaz às condições decontorno em deslocamento do problema:

Substituindo

na expressão (3.121), e integrando-se, chega-se a:

Vale observar que agora Π é uma função doparâmetro α1 e não mais da função ν(x). Isso significaque o problema a ser resolvido é um problema demínimo de função e não mais de mínimo de umfuncional. Essa é a contribuição do método aproximadode Rayleigh-Ritz.

Aplicando-se agora o princípio da mínima energiapotencial total, o qual afirma que a configuraçãodeformada minimiza a energia potencial total de uma

(3.127)

(3.128)

(3.129)

(3.130)

(3.131)

estrutura em equilíbrio estável, obtém-se:

logo

e, portanto,

Observa-se que o erro no cálculo de d em relação àsolução exata é muito grande:

Da resistência dos materiais sabe-se que:

(3.132)

(3.133)

Assim, no trecho (a),

A Figura 3.14 compara os momentos dasolução aproximada e da solução correta (vigaisostática). Os momentos são constantes ao longo de xnos dois trechos porque ν(x) é uma função do segundograu.

FIGURA 3.14 Diagrama de momentos na vigaassociado a ν(x) definido em (3.128).

Observação: a solução é ruim tanto em termos dedeslocamentos quanto em termos de momentos. Aaproximação dos momentos é ainda pior porque ela é

(3.134)

(3.135)

obtida de derivadas de funções aproximadoras.Segunda tentativa:No problema estudado a solução é muito simples

porque a viga é isostática. No caso de uma vigaaltamente hiperestática de vários vãos com inérciasdiferente em cada vão e cargas distribuídas, a soluçãonão é trivial e não estará disponível para se saber se asolução aproximada é boa ou não. Nesse caso, oprocedimento a seguir é usar uma funçãoaproximadora mais “rica” e verificar a mudança naresposta. Quando, ao se refinar a solução, a respostanão melhora significativamente, a solução anterior jápode ser considerada boa.

Na segunda tentativa, a função aproximadora é umpolinômio do terceiro grau dado por:

Vale observar que a função satisfaz às condições decontorno em deslocamento (3.123) e rotação (3.124).

Substituindo

em (3.121) e integrando-se, chega-se a:

(3.137)

(3.138)

(3.139)

(3.140)

(3.136)Vale observar que P agora é uma funçãodos parâmetros α1 e α2. Aplicando-se o princípio damínima energia potencial total, obtém-se:

Que fornece,

Logo,

(3.142)

(3.141)

(3.143)

(3.144)

(3.145)

Usando-se (3.131), chega-se a:

A comparação entre os momentos dasolução aproximada e da solução exata (vigaisostática) está apresentada na Figura 3.15.

FIGURA 3.15 Diagrama de momentos na vigaassociado a ν(x) definido em (3.140).

Observações:1) A solução melhorou significativamente em

termos de deslocamentos, mas continua ruimem termos de momentos. Não é coincidênciaque o deslocamento na extremidade livre se jainferior ao da solução exata, pois a aproximaçãotorna a estrutura mais rígida.

2) O problema na descontinuidade no diagrama demomentos na solução aproximada continua. Adescontinuidade acontece porque ν(x) e ,consequentemente, sua segunda derivada, écontínua no domínio enquanto que a rigidez EI édescontínua em x = 5.

3) O problema identificado revela uma limitação dométodo de Rayleigh-Ritz que é o de trabalharcom apenas uma função contínua no domínio.Para se superar o problema é preciso usar duasfunções, uma no trecho (a) e outra no trecho (b),impondo condições de continuidade em x = 5

(3.146)

(3.147)

(3.148)

(3.149)

para ν(x) e para sua primeira derivada emrelação a x, mas, liberando a curvatura para serdescontínua.

Terceira tentativa:Serão usadas duas funções cúbicas aproximadoras,

uma para o trecho (a) e outra para o trecho (b):

Vale observar que a função νa(x) satisfazàs condições de contorno em deslocamento definidasem (3.123) e (3.124). Além disso, serão impostas asseguintes condições de continuidade em x = 5.

Essas duas condições permitem reduzir o número deparâmetros incógnitos de 6 para 4. Os parâmetros α5 eα6, por exemplo, podem ser escritos em função dosoutros parâmetros.

(3.150)

(3.151)

(3.152)

(3.153)

Aplicando-se o princípio da mínima energia potencialtotal, obtém-se:

É possível obter os parâmetros a1, a2, a3 e a4 que,substituídos em (3.146) e (3.147), fornecem:

(3.154)

(3.155)

(3.156)

Nota-se que

é a solução exata para o deslocamento naextremidade livre. O diagrama de momentoscorrespondentes às expressões (3.154) e (3.155)também é exato.

Observações:a) O uso de duas funções aproximadoras νa(x) e

νb(x) permitiu obter a solução exata do problemaporque foi possível representar adescontinuidade que existe na derivada segundada função elástica em x = 5. O procedimentousado na terceira tentativa foi o de melhorar aprecisão da solução usando duas funçõesaproximadoras, uma para cada trecho da viga,em vez de continuar a aumentar o grau dopolinômio da função ν(x) no domínio de 0 a L.Mesmo usando um polinômio do quarto graupara ν(x) não se pode obter a solução exata

porque haverá ainda uma descontinuidade nasegunda derivada de ν(x) o causará umadescontinuidade no diagrama de momento, umavez que há uma descontinuidade na rigidez EI daviga.

b) Posto como está, o método de Rayleigh-Ritzainda não é um método dos deslocamentos, nosentido clássico, porque as incógnitas não são osdeslocamentos.

c) Com o uso de duas funções no domínio o métododeu um grande passo para se aproximar dométodo dos elementos finitos. Na verdade odomínio foi “discretizado” em dois subdomínios,ou elementos.

d) Para transformar definitivamente o método deRayleigh-Ritz no MEF, o método de Rayleigh-Ritzprecisa substituir as incógnitas a i pelos graus deliberdade da estrutura di.

3.8 O MEF para vigasA ideia básica do modelo de elementos finitos para aestrutura consiste em usar funções aproximadoras,descritas em subdomínios ou elementos finitos, paradescrever os campos de deslocamento da estrutura. Amelhora da solução deve ser obtida com o uso de maissubdomínios ou elementos e não apenas com o uso depolinômios de mais alto grau. Para sistematizar asoperações matemáticas do problema, as funções

(3.157)

aproximadoras devem ser descritas em cadasubdomínio por funções de interpolação previamentedefinidas.

Para o trecho de viga de comprimento L que,posteriormente, será denominado elemento finito deviga, representado na Figura 3.16, escreve-se,inicialmente, a função aproximadora de terceiro grauem função dos parâmetros do polinômio α i, i = 1, ..., 4.

FIGURA 3.16 Elemento finito de viga.

Para se escrever a função aproximadora em função

dos deslocamentos nodais, as seguintes condições decontorno são impostas de acordo com a Figura 3.17:

FIGURA 3.17 Funções de forma ou de interpolação doelemento finito de viga.

(3.158)

(3.159)

(3.160)

Ou, matricialmente,

A solução de (3.159) fornece os α i emfunção dos deslocamentos nodais di. Substituindo os α iobtidos da solução de (3.159) em (3.157), chega-se a:

sendo,

(3.161)

As funções são denominadas funções deinterpolação de viga. Qualquer função ν(x) quedescreva a elástica de um trecho de viga pode serescrita em função das funções ϕi(x). As funções deinterpolação de viga têm um significado cinemáticoque é comum a todas as funções de interpolação detodos os elementos finitos: “ϕi(x) representa o campode deslocamentos no interior do elemento para di = 1mantendo-se todos os outros deslocamentos dj = 0, j ≠i.”

Na Figura 3.17, onde as quatro funções de viga fi(x)estão representadas, pode-se constatar essapropriedade das funções de interpolação.

Para ilustrar o uso das funções de interpolação deviga ϕi(x) na análise de uma viga, o exemplo da Figura3.17 será reanalisado com o uso dessas funções deinterpolação.

(3.162)

(3.163)

Observando que se deve usar L = 5 (comprimento decada trecho) nas expressões de ϕi(x) para se obter asfunções νa(x) e νb (x), pode-se escrever para a viga daFigura 3.18:

FIGURA 3.18 Graus de liberdade da viga em balançomodelada por 2 elementos finitos.

Considerando as duas funções distintasνa(x) e νb(x), respectivamente nos trechos (a) e (b), eintegrando-as em x de 0 a L = 5 em cada trecho(subdomínio do trecho ou elemento finito) e

(3.164)

(3.165)

observando-se que a força P atua no sentido negativoda direção de d3 , pode-se escrever a expressão daenergia potencial total Π(d1, d2, d3, d4) como:

Substituindo na expressão (3.166) EIa, EIb eP pelos seus valores numéricos, efetuando as integraise usando o princípio da mínima energia potencial totalcomo descrito em (3.109), obtém-se:

ou,

(3.166)

(3.167)

(3.168)

(3.169)

Sendo K a matriz de rígidez da viga, d ovetor dos deslocamentos nodais e f o vetordas cargas nodais. Essa solução é “exata” e coincidecom a última solução obtida para o método deRayleigh-Ritz.

Para se obter o sistema de equações equivalente aosistema (3.165), mas para apenas um elemento decomprimento L e rigidez EI, fazendo uso das funções deinterpolação de viga ϕi(x) e com os deslocamentosnodais de di conforme descrito na Figura 3.17, repete-seo procedimento descrito a partir de:

Seguindo os mesmos passos anteriorescom Π(d1, d2, d3, d4) dado agora por:

Chega-se a:

(3.171)

(3.170)Observe que a matriz K obtida em (3.170)é a mesma matriz de rigidez do elemento de viga daanálise matricial de estruturas e que o termo Kij podeser obtido de:

3.9 O método dos resíduosponderados de galerkinComo visto no item anterior, quando existe umfuncional e um correspondente princípio de mínimoassociado a um dado problema de engenharia, o MEFpode ser formulado com as funções que representamos campos incógnitos descritas por funções deinterpolação de variáveis nodais.

(3.172)

(3.173)

(3.174)

Alternativamente, as equações do MEF podem serobtidas diretamente das equações diferenciais doproblema. A vantagem desse enfoque é que o métodopode ser aplicado a uma gama de problemas para osquais não há um funcional disponível.

Se ja um problema unidimensional representado pelaequação diferencial dada a seguir, onde u(x) é umafunção incógnita no domínio do problema.

Com as condições de contorno dadas por:

Uma função aproximadora ua(x) que aproxima u(x) nodomínio do problema e que contenha n parâmetrosincógnitos pode ser escrita como:

Onde Ni(x) são as funções de interpolação dasvariáveis nodais di que passam a ser as incógnitas doproblema.

Matricialmente, a expressão (3.176) pode ser

(3.175)

(3.176)

reescrita como:

A expressão (3.175) representa a funçãoaproximadora aplicada em um subdomínio ou em umelemento no MEF. Ao usar (3.175), o campo u(x) estásendo representado por um elemento apenas. A opçãode usar somente um elemento no domínio é feita aquiapenas para simplificar a apresentação do método dosresíduos ponderados, mas não é uma limitação dométodo. Em geral, vários elementos podem ser usadospara representar o campo das funções incógnitas.

Quando se substitui a função aproximadora ua(x)dada em (3.175) na expressão (3.172), a expressão nãodeve satisfazer a igualdade em todo o domínio doproblema fornecendo o que se costuma chamar defunção resíduo R(x) da solução.

Os melhores valores de di são aqueles que reduzema função resíduo R(x) de uma forma integral nodomínio do problema. Como, no entanto, há nincógnitas di para o problema, são necessárias nequações para obtê-las. Uma maneira de se obter as nequações é usar n funções de ponderação W i(x) e ,

(3.177)

(3.178)

consequentemente, n equações da forma:

A expressão (3.177) é a equação fundamental dométodo dos resíduos ponderados. No método deGalerkin usa-se:

Ou se ja, as funções de ponderação W i(x) são iguaisàs funções de interpolação Ni(x).

3.9.1 Exemplos De Aplicação Do MétodoDe Galerkin3.9.1.1 Equação de equilíbrio de umabarra de treliçaSeja a barra de treliça tracionada com área da seçãotransversal A e módulo de elasticidade do material Erepresentada na Figura 3.20.

FIGURA 3.19 Barra tracionada.

FIGURA 3.20 Funções de interpolação dodeslocamento longitudinal u(x).

A equação de equilíbrio das forças horizontais paraum elemento dx da barra é dada por:

(3.180)

(3.182)

(3.183)

(3.179)

(3.181)

ou

A lei de Hooke fornece:

E a equação de compatibilidade:

logo, substituindo (3.182) em (3.181) e , em seguida,(3.181) em (3.180), chega-se a:

(3.184)

(3.185)

(3.186)

Que representa a equação diferencial do problemano domínio 0 ≤ x ≤ L. A solução de equações diferenciaisconduz a dois tipos de problema, nomeadamente:problema de valor de contorno e problema de valoresiniciais. O problema de valor de contorno só necessitada especificação das condições de contorno para suasolução, enquanto o de valores iniciais precisa tambémdas condições iniciais das variáveis definidas noespaço do tempo. A expressão (3.183) representa umproblema de valor de contorno. As condições decontorno naturais desse problema são:

Seja o campo u(x) aproximado pela função deinterpolação ua(x) dada por:

Sendo as funções de interpolação dadaspor,

(3.187)

(3.188)

(3.190)

(3.189)

Os deslocamentos u1 e u2 são os dois parâmetrosincógnitos. A Figura 3.20 ilustra as funções deinterpolação N1(x) e N2(x) da função aproximadora ua(x).

Nesse caso, a expressão (3.176) vale :

ou

Onde cada sobrescrito vírgula representa umaderivada da função em relação a x. O método deGalerkin fornece duas equações:

(3.191)

(3.192)

(3.193)

(3.194)

(3.195)

A derivada do produto de duas funçõesf(x) e g(x) em relação a x, vale :

logo,

Integrando os dois lados de (3.193) de 0 aL, obtém-se:

ou

A expressão (3.195) é conhecida namatemática como técnica de integração por partes.Se jam:

(3.196)

(3.198)

(3.199)

(3.197)

Considerando (3.196), (3.191) e (3.195), chega-se a:

Examinando-se a primeira parcela àdireita de (3.197) e considerando que:

Obtém-se, para i = 1:

(3.203)

(3.200)

(3.201)

(3.202)

E, para i = 2:

A segunda parcela à direita de (3.197)pode ser reescrita como:

onde,

Considerando agora (3.200), (3.201) e (3.202), épossível reescrever (3.197) na forma matricial como:

(3.204)

(3.205)

(3.206)

onde K, d e f são, respectivamente, a matriz derigidez do elemento de treliça no sistema local decoordenadas, o vetor dos deslocamentos nodais e ovetor das cargas nodais, dados por:

O sistema de equações linearesalgébricas dado em (3.204) representa as equações deequilíbrio de uma barra de treliça no seu sistema local.

3.9.1.2 Equações de equilíbrio de umabarra de vigaA equação diferencial de equilíbrio de uma viga semcargas atuantes é dada por:

E é o módulo de elasticidade do material da viga, I éo momento de inércia da seção transversal e ν(x) afunção que descreve os deslocamentos transversais

(3.207)

(3.208)

(3.209)

(3.210)

da viga. As condições de contorno naturais doproblema são:

Usando-se a função aproximadora dada em (3.168),

Onde as funções de interpolação ϕi(x) são as que

(3.211)

estão descritas em (3.161) e repetindo-se oprocedimento análogo ao que foi adotado no itemanterior, ou se ja, aplicando-se o método de Galerkin, épossível chegar a um sistema de equações linearesalgébricas que representa as equações de equilíbrio deuma barra de viga, sendo agora o elemento Kij damatriz de rigidez dado por:

Observa-se que a expressão dada em (3.211)coincide com o resultado obtido em (3.171) para ocoeficiente de rigidez Kij obtido pelo MEF.

3.10 Generalização do MEF3.10.1 Formulação Geral Do MEFO MEF descrito no item anterior será generalizadonesse item de modo que ele possa ser aplicadotambém a estruturas contínuas bi e tridimensionais.Inicialmente, as seguintes hipóteses são introduzidaspara o contínuo:

a) O contínuo é idealizado como formado porelementos com diferentes formas geométricas,como triângulos, quadriláteros, tetraedros etc.

(os elementos finitos), ligados por alguns nóssituados no contorno.

b) Matrizes de interpolação para o elemento m(matriz Nm),cujos termos são funções conhecidascomo funções de interpolação ou de forma, quefornecem os campos de deslocamento (vetor um)no interior dos elementos em função dosdeslocamentos nodais do elemento dm, ou se ja,um = Nm dm.

c) O vetor das deformações no interior doselementos (vetor σm) pode ser obtido porderivação dos campos de deslocamentos um emrelação às coordenadas do sistema gerando aexpressão εm = Bm dm.

d) As tensões no interior dos elementos (vetor σm)são obtidas a partir das deformações por meiode relações constitutivas. Para um corpohomogêneo e um material de comportamentolinear elástico, é possível definir apenas umamatriz constitutiva C que relaciona asdeformações a as tensões no elemento por σm =C εm. Para materiais isotrópicos, os termos damatriz C dependem apenas dos seguintesparâmetros mecânicos do material: E, módulo deelasticidade longitudinal e ν, coeficiente dePoisson.

e) Uma matriz de rigidez (matriz Km) e um vetor decargas equivalentes nodais fm para o elementopodem ser obtidos a partir das matrizes geradas

Nm, Bm e C.f) As matrizes de rigidez e as cargas nodais

equivalentes de cada elemento são combinadasadequadamente de forma a montar a matriz derigidez global Kg e o vetor global de cargasnodais fg da estrutura.

g) Os deslocamentos globais são calculados daequação de equilíbrio global da estrutura Kg dg =fg.

A partir do vetor dos deslocamentos global daestrutura dg é possível recuperar deslocamentosnodais de cada elemento dm, e , em seguida, os calculartodas as respostas da estrutura em termos dedeformações e tensões em qualquer ponto daestrutura fazendo uso das equações anteriores.

A expressão da matriz de rigidez para um elementofinito pode ser obtida por meio da expressão (3.100) doprincípio dos deslocamentos virtuais generalizado. Elaestá repetida a seguir, mas em uma forma mais geralpara que os vetores das forças volumétricas q,superficiais p e das forças nodais f, possam serconsiderados na formulação do problema. Além disso,considera-se que a estrutura possui apenas umelemento, de modo que se adota nb = m = 1 em (3.100).

(3.212)

(3.213)

(3.214)

(3.215)

(3.216)

(3.217)

As integrais em V1 e Γ1 significam,respectivamente, integrais no domínio e nocontorno do elemento 1. Nas expressões a seguir, seráabandonado o subescrito 1 referente ao únicoelemento por questão de simplicidade. Lendo comatenção os itens de (a) a (h) descritos anteriormente,as seguintes expressões podem ser escritas relativasàs grandezas virtuais:

Para as grandezas reais, as expressões são

Substituindo as expressões de (3.213) a (3.217) naexpressão (3.212), chega-se a:

(3.220)

(3.221)

(3.222)

(3.218)

(3.219)

Por ser arbitrário, o vetor dosdeslocamentos virtuais nodais transposto δdt queaparece nos dois lados da expressão (3.215) pode sereliminado da equação, o que resulta em:

onde

Nas expressões mencionadas, K é a matriz de rigidezdo elemento, fq o vetor das forças nodais equivalentesàs cargas de volume, fp o vetor das forças nodaisequivalentes às cargas de superfície e f o vetor dasforças nodais propriamente ditas.

Uma vez calculado d em (3.219), u, ε e σ, podem serobtidos pelas expressões (3.215), (3.216) e (3.217),respectivamente.

3.10.2 Critérios De Convergência DoMEFAs funções de interpolação ou funções de forma damatriz de interpolação N são como funçõesaproximadoras do método de Rayleigh-Ritz. Adiferença entre as duas funções está nos parâmetrosincógnitos. Enquanto no método de Rayleigh-Ritz osparâmetros incógnitos são coeficientes generalizados,no MEF, eles são os deslocamentos nodais. Nos doismétodos, todavia, as funções tentam aproximar assoluções “exatas”. No método de Rayleigh-Ritz,melhores soluções são obtidas quando se usampolinômios de mais alto grau ou uma sérietrigonométrica com mais termos como funçõesaproximadoras. No MEF, resultados mais precisos sãoesperados quando se usa uma malha mais refinada deelementos. Não podemos esquecer, todavia, que emambos os casos as soluções são aproximadas. Umaestrutura modelada por elementos finitos é uma

estrutura, em geral, mais rígida que a estrutura realporque as funções aproximadoras usadas pararepresentar os campos de deslocamento, na maioriadas vezes, não conseguem reproduzir o campo real dedeslocamentos. Assim, elas impõem restrições à livredeformação da estrutura de modo que ela possaminimizar sua energia potencial total. Para garantir queas soluções convirjam para a solução “exata” no MEFalguns critérios devem ser satisfeitos na sua escolha.

Antes de apresentar os critérios de convergência aserem satisfeitos para que a solução aproximada viaMEF convirja para a solução “exata”, é convenientefalar de completidade de um polinômio. Um polinômiof(x,y) no espaço bidimensional será utilizado paraesclarecer o conceito de completidade de um polinômio.O triângulo de Pascal mostrado na Figura 3.21 ilustra ostermos de um polinômio completo do terceiro grau.

(3.223)

FIGURA 3.21 Triângulo da Pascal com os termos dopolinômio completo p(x,y).

O polinômio dado a seguir é um polinômioincompleto do terceiro grau, pois não possui todos ostermos do polinômio completo, ou se ja:

Para que se ja completo ele deveria tertambém os termos x2y e yx2.

Se ja ϕ(x,y,z) uma função aproximadora que descreveo campo de deslocamentos de um elemento eP(ϕ(x,y,z)) a energia potencial total de uma estruturamodelada com elementos formulados com essasfunções. Suponha que Π(ϕ(x,y,z)) contenha derivadasde ordem m de ϕ(x,y,z). Para que a solução se aproximeda solução exata da estrutura quando se refina amalha com esses elementos, alguns critérios precisamser satisfeitos. Esses critérios, denominados critérios

de convergência, são os seguintes: Critério 1: dentro de cada elemento, a função

aproximadora ϕ(x,y,z) precisa conter umpolinômio completo de grau m.

Critério 2: na fronteira entre elementos, devehaver continuidade ϕ(x,y,z) e de suas derivadasaté a ordem m − 1.

Critério 3: se ja uma malha de elementossubmetida a condições de contorno compatíveiscom valores constantes de qualquer dasderivadas de ordem m de ϕ(x,y,z). Então, quandoa malha é refinada, cada elemento devereproduzir essas deformações constantes.

O critério 1 assegura que ϕ será contínua dentro doelemento e é necessário (mas nem sempre suficiente)para garantir a satisfação do critério 3.

O critério 2 é satisfeito para qualquer malha paraelementos compatíveis (conforming elements). Elementosincompatíveis (nonconforming elements) devem se tornarcompatíveis com o refinamento da malha. Muitoselementos bem sucedidos de placa à flexão sãoincompatíveis com relação às rotações na fronteirapara malhas pouco refinadas, mas tornam-secompatíveis com o refinamento da malha.

O critério 3 é satisfeito para a grande maioria doselementos. Ele pode ser usado para testar os critérios 1e 2 no que se denomina “Patch test” na literatura.

3.10.3 Montagem Da Matriz De Rigidez

Global E Do Vetor Global De CargasEquivalentes NodaisO MEF herdou da análise matricial de estruturasreticuladas a técnica de montagem da matriz de rigidezglobal da estrutura a partir da contribuição dasmatrizes de rigidez local de cada elemento da malha.

Isso foi possível porque no modelo de elementosfinitos, assim como na análise matricial de estruturas, aestrutura é representada por elementos conectadosentre si por meio de nós. Na análise matricial deestruturas, os elementos são barras com dois nós nasextremidades, no método dos elementos finitos, e lessão polígonos ou poliedros e os nós estão em geral nocontorno (normalmente nos vértices) dos elementos.As palavras polígono, poliedro e vértice são usadas aquide forma livre, uma vez que alguns elementos podemter lados curvos ou superfícies laterais curvas. Hátambém a possibilidade de haver nós no interior doelemento em alguns elementos finitos, mas esse casotambém se adapta bem à técnica de montagem.

Duas ideias básicas são utilizadas para a montagemda matriz de rigidez global:

a) Associar a cada nó da estrutura graus deliberdade que dependem do número do nó.

Em análise de estruturas bidimensionais, porexemplo, cada nó tem dois graus de liberdade,nomeadamente deslocamento horizontal (e ixo x) edeslocamento vertical (e ixo y). Ao nó de número i naestrutura são, normalmente, associados os graus de

liberdade 2i – 1, na direção horizontal, e 2i, na direçãovertical. Para estruturas tridimensionais em que cadanó tem três graus de liberdade associados,respectivamente, aos eixos x, y e z, os deslocamentosdo nó i são associados aos graus de liberdade 3i – 2, 3i –1 e 3i. A Figura 3.22 esclarece.

FIGURA 3.22 Graus de liberdade associados ao nó i notriângulo e no tetraedro.

b) Associar a numeração local dos nós em cadaelemento à numeração global dos nós a nível global, ouseja, da estrutura.

A associação da numeração local dos nós doelemento com a numeração global dos nós daestrutura é feita pela matriz de incidência Inc. A matrizInc indica como o elemento se conecta com a estrutura.

(3.224)

Supondo-se que o modelo da estrutura tem neelementos e que cada elemento tem n nós, a matriz deincidência teria ne colunas e n linhas. Se ja o elementotriangular de número m com três nós nos vértices doelemento com numeração local 1, 2 e 3 no sentido anti-horário (no próximo capítulo será mostrado porque senumeram os nós do elemento triangular dessamaneira). A escolha de qual será o nó 1 é arbitrária,mas, uma vez definido, as posições dos outros doisficam determinadas pela regra do sentido anti-horário.Na malha de elementos finitos da estrutura os nóscorrespondentes têm numeração 10, 12 e 15, porexemplo. Para esse exemplo, a coluna m da matriz Incteria 10, 12 e 15 na primeira, segunda e terceira linha,respectivamente, como indicado a seguir:

A Figura 3.23 ilustra a relação entre asnumerações local e global dos nós do elemento.

FIGURA 3.23 Relação entre as numerações local eglobal dos nós do elemento m.

Com a regra (a) os graus de liberdade do elemento(local) e global seriam os indicados na Figura 3.24.

FIGURA 3.24 Relação entre os graus de liberdadelocal e global do elemento m.

Com as regras (a) e (b) é possível agora criar umaexpressão que relaciona os 6 graus de liberdade nosistema local com os 6 graus de liberdade do sistemaglobal para os 3 nós do elemento m que serãoarmazenados na matriz de ponteiros dg cujo elementodg i,m representa o grau de liberdade na direção globalcorrespondente ao grau de liberdade do elemento i =1,..,6 do elemento m = 1,...,ne, ou se ja:

(3.225)

(3.226)

Usando os valores de incidência Incdefinidos nacoluna m da expressão (3.224), a expressão (3.225)passa a armazenar os seguintes valores:

Os valores correspondentes aos graus de liberdadeglobais da expressão (3.226) são os valoresrepresentados na Figura 3.25.

(3.227)

Com os graus de liberdade do sistema globalarmazenados na matriz dg pode-se proceder amontagem sistemática da matriz de rigidez global Kg apartir das contribuições das matrizes dos elementos nosistema local Ke como indicado a seguir:

Observa-se que, pela expressão (3.227), otermo da matriz de rigidez local do elemento mda Figura 3.25, por exemplo, será somado ao termoK20,29 da matriz global. Vale notar que esse mesmotermo da matriz de rigidez global pode recebercontribuições de outros elementos da malha e por issoa rigidez global deve ser acumulada com ascontribuições das matrizes de rigidez do sistema locaisdos elementos.

Analogamente, o vetor global de cargas nodaisequivalentes fg deve ser acumulado com ascontribuições dos vetores das cargas equivalentesnodais que atuam no elemento fe pela expressão.

(3.228)

C A P Í T U L O 4

Problemas de estado plano

4.1 IntroduçãoEstruturas de estado plano ou chapas são estruturasbidimensionais, ou se ja, aquelas em que umadimensão, denominada espessura t , normalmentemedida na direção do eixo z do sistema decoordenadas cartesiano, é muito menor do que asoutras duas medidas nas direções dos eixos x e y doplano xy como indicado na Figura 4.1. As forçasatuantes nessas estruturas agem também no plano xy.

Exemplos de estrutura de estado plano em Engenharia Civil são as vigas-parede e uma fatia deespessura constante de uma barragem de gravidadeou de um muro de arrimo.

4.1.1 Equações De CompatibilidadeOs campos de deslocamento dessas estruturas sãou(x,y) e υ(x,y), respectivamente nas direções dos eixosx e y como mostra a Figura 4.1.

FIGURA 4.1 Chapa plana de espessura t com planomédio no plano xy.

As componentes do vetor de deformações deinteresse nesse problema são εx, εy e γxy,nomeadamente deformações longitudinais nasdireções dos eixos x e y e a deformação de distorçãono plano xy.

As equações de compatibilidade são dadas por:

(4.1)

(4.2)

(4.3)

Ou, matricialmente,

ou, sucintamente,

Onde ε é o vetor das deformações, L amatriz operadora de derivação e u o vetor dosdeslocamentos.

4.1.2 Equações constitutivasExistem dois tipos de estruturas planas: as que estãoem estado plano de tensão e as que estão em estadoplano de deformação.

As estruturas em estado plano de tensãocaracterizam-se por apresentarem tensão normal nula,σz = 0, na direção do eixo z normal ao plano xy onde selocaliza sua superfície média. As deformaçõesassociadas à direção do eixo z são livres, nesse caso, εz≠ 0.

As estruturas em estado plano de deformação, aocontrário, têm deformação impedida na direção do eixoz, εz = 0, e tensões normais, σz ≠ 0.

4.1.2.1 Estado plano de tensão

As componentes de tensão que caracterizam umestado plano são σx, σy e τxy, respectivamente, a tensãonormal na direção x, a tensão normal na direção y e atensão cisalhante.

Dois parâmetros são suficientes para descrever ocomportamento mecânico das estruturas em estadoplano, caso elas tenham um comportamento isotrópicoe linear elástico, que são o módulo de elasticidadelongitudinal E e o coeficiente de Poisson ν.

Se um elemento infinitesimal dxdy de uma estruturaem estado plano de tensões no plano xy estiversubmetido ao seguinte estado de tensão longitudinal σx

(4.4)

(4.5)

(4.6)

≠ 0 e σy = τxy = 0, ou se ja, a um estado uniaxial de tensãona direção x, estando a estrutura livre para sedeformar em todas as direções, e la irá se alongar nadireção x segundo a le i de Hooke:

E encurtar-se na direção y, segundo o efeito dePoisson:

Devido à propriedade de isotropia e decomportamento linear elástico do material, quepermite a superposição de efeitos, caso o estado detensão se ja σx ≠ 0 e σy ≠ 0, ou se ja, um estado biaxial detensões, obtém-se:

(4.7)

Se o elemento estiver submetido a dois binários detensões cisalhantes que, por condições de equilíbrio,devem ser autoequilibrados como representado naFigura 4.2, seus ângulos retos irão fechar ou abrir de:

FIGURA 4.2 Representação da deformação porcisalhamento.

O parâmetro G é denominado módulo de

(4.8)

(4.9)

(4.10)

elasticidade transversal e depende de E e ν comoindicado a seguir:

As relações constitutivas podem ser escritasmatricialmente:

ou, sucintamente,

Onde D é a matriz de elasticidade referenteao estado plano de tensão.

A inversa da matriz D re laciona o vetor das tensõescom o vetor das deformações:

(4.12)

(4.11)A equação (4.11) representa o que se chama de lei de

Hooke generalizada, pois ela representa umageneralização para o caso de estado plano de tensãoda lei de Hooke unidimensional. Nessa equação, C = D−1

é denominada matriz constitutiva para o estado planode tensão, dada por:

4.1.2.2 Estado plano de deformação

Como visto anteriormente, o estado plano dedeformação é caracterizada por, σz ≠ 0 e εz = 0, logo, apartir de:

(4.13)

(4.14)

(4.15)

E, considerando εz = 0, obtém-se:

Substituindo (4.14) nas duas primeiras equações(4.13), obtém-se:

ou, ainda, matricialmente:

(4.17)

(4.18)

ou, sucintamente,

Multiplicando ambos os lados da expressão(4.17) por D−1 = C, chega-se a uma equação análoga àexpressão (4.11) onde a matriz C é agora denominadamatriz constitutiva para uma estrutura em de estadoplano deformações, dada por:

4.2 Elemento Triangular DeDeformação ConstanteO elemento triangular de 3 nós representado na Figura4.3 é chamado na literatura de triângulo de deformação

(4.19)

constante (CST, constant strain triangle) por razões queserão explicadas mais adiante.

FIGURA 4.3 Elemento CST (constant strain triangle).

Os campos de deslocamento que descrevem osdeslocamentos no interior do elemento são polinômioslineares em x e y da forma:

(4.20)

(4.21)

ou, matricialmente,

ou, sucintamente,

A escolha de polinômios lineares de 3 termos com 6coeficientes incógnitos a i, 3 para cada campo dedeslocamento, pode agora ser justificada pelas 6condições de contorno a seguir:

(4.22)

(4.23)

que podem ser reescritas usando-se a expressão(4.19), ou se ja:

ou, sucintamente,

(4.26)

(4.24)

(4.25)

(4.27)

ou

onde o vetor d contém os deslocamentos nodais, ouseja:

Substituindo a expressão (4.25) em (4.21), obtém-se:

ou, ainda,

(4.31)

(4.28)

(4.29)

(4.30)

Sendo que,

A matriz N(x, y) tem a forma:

Observando-se as equações (4.28) e (4.30) épossível escrever:

onde, ui e υi são os deslocamentos nodais relativosaos eixos x e y, respectivamente.

As funções de interpolação Ni(x,y) são dadas por:

(4.33)

(4.34)

(4.32)onde a variável área representa a área doelemento triangular cuja expressão pode ser obtidapelo seguinte determinante dado:

Deve-se ter o cuidado de numerar os nós doelemento no sentido anti-horário para que o resultadoda expressão (4.33) se ja positivo.

Como visto anteriormente, a matriz de rigidez de umelemento finito qualquer pode ser obtida por:

(4.35)

(4.36)

(4.37)

(4.38)

onde B é a matriz de compatibilidade cinemática quetransforma deslocamentos nodais em deformações nointerior do elemento

e C é a matriz constitutiva que transforma o vetor dedeformações ε em vetor de tensões σ para o materialde comportamento linear elástico (le i de Hooke) e aintegral é efetuada no volume do elemento.

Substituindo u(x, y) e υ(x, y) dado em (4.31)em (4.3), obtém-se:

onde L é a matriz operadora de derivação definidaem (4.3). A expressão (4.37) pode ser reescrita como:

o que permite concluir que para o elemento emquestão, a matriz B é obtida aplicando-se a matriz

(4.39)

(4.40)

operadora de derivação L à matriz N(x, y),

resultando em:

onde Ni(x, y),x e Ni(x, y),y representam aderivada de Ni(x,y) em relação à x e y,respectivamente, e estão explicitadas a seguir:

(4.41)

(4.42)

A matriz de rigidez é dada por:

no caso de a espessura t ser constante na área doelemento. Nesse contexto, a integração da matriz derigidez é trivial devido ao fato de a matriz B serconstante, ou se ja, independente de x e y para oelemento em questão, o que permite reescrever (4.42)como:

(4.44)

(4.44)Caso a espessura t não se ja constante no

interior do elemento ela pode ser representada poruma interpolação dos valores nodais t i analogamenteàquela adotada para os campos de deslocamento.

Nesse caso, seria necessário recorrer a umaintegração numérica para efetuar a integração. Umaalternativa à integração numérica, que fornece bonsresultados, seria adotar uma espessura média igual aum terço da soma dos valores nodais. Para estruturasem estado plano de deformação adota-se t igual a umaunidade de comprimento.

4.3 Elementos da famíliaSerendipityUm elemento finito é dito isoparamétrico quando asmesmas funções de interpolação são usadas parainterpolar não apenas grandezas cinemáticas(deslocamentos), como é usual nos elementosconvencionais, mas também as grandezas

geométricas (no caso coordenadas). A família desseselementos é chamada de família Serendipity, umareferência ao conto persa infantil Os três príncipes deSerendip. Esta história conta as aventuras de trêspríncipes do Ceilão, atual Sri Lanka, que viviamfazendo descobertas inesperadas, cujos resultadoseles não estavam procurando realmente. Graças à suacapacidade de observação e sagacidade, descobriam“acidentalmente” a solução para dilemas impensados.Essa característica tornava-os especiais e importantes,não apenas por terem um dom especial, mas porterem a mente aberta para as múltiplas possibilidades.

Como será visto mais adiante, para um elementoquadrilateral como o representado na Figura 4.4, oscampos que descrevem as coordenadas cartesianasdevem ser polinômios de 4 termos em coordenadasparamétricas para cada coordenada.

(4.45)

FIGURA 4.4 Elemento isoparamétrico de 4 nós dafamília Serendipity.

Assim, os polinômios paramétricos são:

ou, matricialmente,

(4.46)

(4.47)

ou, sucintamente,

A escolha de polinômios de 4 termos com 8coeficientes incógnitos a i pode agora ser justificadapelas 8 condições de contorno seguintes:

(4.48)

(4.49)

que podem ser reescritas pela expressão (4.46).

ou, sucintamente,

(4.50)

(4.51)

(4.52)

(4.53)

(4.54)

ou

onde o vetor d contém os deslocamentos nodais.Substituindo a expressão (4.51) em (4.47), obtém-se:

ou, ainda,

sendo,

A matriz N(ξ,η) tem a forma:

(4.56)

(4.57)

(4.55)Observando-se as equações (4.55) e (4.53) épossível escrever:

Analogamente, as coordenadas x(ξ,η) e y(ξ,η) podemser escritas em função das coordenadas nodais, ouseja:

onde, xi e yi são as coordenadas nodais e ui e υi sãoos deslocamentos nodais relativos aos eixos x e y,respectivamente.

(4.58)

As funções de interpolação Ni(ξ,η) são dadas por:

As expressões (4.57) permitem mapear um pontoP(ξ, η) do quadrado representado no plano paramétricopara um ponto P(x, y) no quadrilátero representado noplano cartesiano como indicado na Figura 4.4.

A Figura 4.5 apresenta a representação geométricada função de interpolação N1(ξ, η). As outras funçõesde interpolação Ni(ξ, η) são análogas a N1(ξ, η), ou se ja,elas valem 1 no nó i e 0 nos outros nós.

(4.59)

FIGURA 4.5 Função de interpolação N1(x, y).

Seja uma função ϕ(x,y). Se x e y forem definidosconforme as expressões (4.56), a relação entre asderivadas de ϕ quanto às coordenadas cartesianas eas derivadas de ϕ no tocante às coordenadasparamétricas é dada pela regra da cadeia:

(4.60)

(4.61)

ou, matricialmente,

Pode-se definir agora a matriz jacobiana J(ξ, η) como:

e, utilizando (4.56), obtém-se:

(4.62)

(4.63)

(4.64)

(4.65)

ou, matricialmente,

onde os subscritos ξ e η significam,respectivamente, a derivada em relação a ξ e η.

Sucintamente, (4.63) pode ser reescrita como:

Observando (4.60), pode-se deduzir que a inversa damatriz jacobiana Γ(ξ,η), dada por,

(4.66)

(4.67)

transforma derivadas paramétricas de ϕ emderivadas cartesianas de ϕ. Sendo assim, pode-seescrever:

observa-se que as submatrizes em (4.66) temdimensão 2 × 2. Sucintamente, a expressão (4.66) podeser escrita como:

onde, u,c é o vetor que contém as derivadascartesianas das componentes de deslocamentos u e υ,u,p o vetor que contém as derivadas paramétricas dascomponentes de deslocamentos u e υ e Γu a matriz quetransforma derivadas paramétricas dosdeslocamentos em derivadas cartesianas dosdeslocamentos.

As expressões (4.56) permitem escrever:

(4.71)

(4.68)

(4.69)

(4.70)

ou, sucintamente:

sendo d o vetor dos deslocamentos nodais.Foi mostrado no Capítulo 2 que o determinante da

matriz jacobiana é o fator de escala que transforma aárea elementar dξdη no quadrado paramétrico emárea elementar correspondente no quadrilátero doplano cartesiano dA, como indicado a seguir.

Como visto anteriormente, a matriz de rigidez de umelemento finito qualquer pode ser obtida por:

(4.72)

(4.73)

(4.74)

onde B é a matriz de compatibilidade cinemática quetransforma deslocamentos nodais em deformações nointerior do elemento. As componentes do vetor dedeformações podem ser escritas em função do vetorque contém as derivadas dos deslocamentos u(x, y) ev(x, y) em relação as coordenadas x e y, ou se ja:

ou, sucintamente,

Usando-se agora (4.67) e (4.69), a expressão (4.73)pode ser reescrita como:

O que permite concluir que a expressão para amatriz de compatibilidade cinemática B(ξ,η), para o

(4.75)

(4.76)

elemento em questão, vale ,

e a matriz de rigidez pode ser dada por:

sendo t a espessura do elemento. Caso tvarie no interior do elemento, e le pode ser interpoladode forma análoga aos deslocamentos e incorporado àintegral A integração da matriz de rigidez é feita noplano paramétrico por integração numérica porque,para o elemento isoparamétrico, as funções emquestão estão definidas nesse espaço.

A integração da matriz de rigidez é feita porintegração numérica pelo método de Gauss. Se foremusados ng pontos de Gauss com coordenadasparamétricas e e pesos de integração e

, a expressão (4.76) pode ser reescrita como:

(4.76)Caso a espessura t varie no interior doelemento, e la pode ser interpolada de formaanáloga aos deslocamentos e incorporada aosomatório que efetua a integração numérica. Oelemento finito isoparamétrico quadrilateral pode terlados curvos quando possui 8 nós, sendo 4 nos“vértices” do quadrilátero e 4 sobre os “lados” comoindicado na Figura 4.6. Isso acontece porque omapeamento do quadrilátero de 8 nós no planoparamétrico se torna um “quadrilátero” de ladoscurvos no plano cartesiano como será mostradoadiante.

FIGURA 4.6 Elemento isoparamétrico de 8 nós dafamília Serendipity.

Para esse elemento de 8 nós, os polinômiosparamétricos devem ter 8 termos para descrever ocampo de deslocamento u(ξ, η) e 8 termos paradescrever o campo de deslocamento υ(ξ, η).Observando o triângulo de Pascal que representa ostermos de um polinômio no plano ξ,η como indicado naFigura 4.4, podem-se escolher os seguintes 8 termospara os dois campos de deslocamento conformeindicado na Equação (4.78).

FIGURA 4.7 Triângulo da Pascal e termos dopolinômio completo p(ξ, η).

A escolha dos 6 primeiros termos do polinômio éobvia, pois eles formam um polinômio de segundo graucompleto. Para completar os 8 termos, é precisoescolher mais 2 termos do terceiro grau. A justificativa

(4.78)

(4.79)

para a escolha dos 2 termos ξ2η e ξη2 pode ser dadapelo fato de serem simétricos e mais neutros dos queos termos ξ3 e η3 e por conterem as duas variáveis ξ eη. A representação matricial fornece:

ou, sucintamente,

A escolha de polinômios de 8 termos com 16coeficientes incógnitos a i pode agora ser justificadapelas 8 condições de contorno seguintes:

(4.81)

que podem ser reescritas usando-se a expressão(4.79), como:

(4.82)

(4.83)

(4.84)

(4.85)

ou, sucintamente,

ou

Onde o vetor d é o vetor dos deslocamentos nodais.Substituindo a expressão (4.84) em (4.80), obtém-se:

(4.89)

(4.86)

(4.87)

(4.88)

Ou, ainda,

sendo

A matriz N(ξ,η) tem a forma:

Observando as equações (4.87) e (4.88) épossível escrever:

(4.90)

(4.91)

Analogamente, as coordenadas x(ξ, η) e y(ξ, η)podem ser escritas em função das coordenadasnodais, ou se ja:

onde, xi e yi são as coordenadas nodais e ui e υi sãoos deslocamentos nodais relativos aos eixos x e y,respectivamente.

As funções de interpolação Ni(ξ, η) são dadas por:

As expressões (4.90) permitem mapear um pontoP(ξ,η) do quadrado representado no plano paramétricopara um ponto P(x,y) no “quadrilátero” representadono plano cartesiano como indicado na Figura 4.6.

Na Figura 4.8, a representação geométrica da funçãode interpolação N1(ξ,η) para o elemento de 8 nós éapresentada como uma combinação da função deinterpolação N1(ξ,η)do elemento de 4 nós (primeirotermo da expressão de N1(ξ,η) para 8 nós) e dasfunções N5 (ξ,η) e N8(ξ,η). As outras funções deinterpolação Ni(ξ,η) são análogas a N1(ξ,η), ou se ja, e lasvalem 1 no nó i e 0 nos outros 7 nós.

FIGURA 4.8 Funções de interpolação N1(ξ, η), N5(ξ, η)e N8(ξ, η) do elemento isoparamétrico de 8 nós.

A matriz jacobiana definida na expressão (4.61) valeagora:

(4.92)

(4.93)

(4.94)

ou, matricialmente,

onde os subscritos ξ e η significam aderivada em relação a ξ e η, respectivamente.

Sucintamente, (4.93) pode ser reescrita como:

Uma vez obtida J(ξ,η), pode-se chegar a Γ(ξ,η):

(4.95)

(4.96)

(4.97)

(4.98)

e as derivadas cartesianas dos deslocamentos por:

Observe que as submatrizes em (4.96) temdimensão 2 × 2. Sucintamente, (4.96) pode ser escritacomo:

As expressões (4.90) permitem escrever:

onde

(4.100)

(4.99)

(4.101)

e, o vetor d,

Como já visto anteriormente, o determinante damatriz jacobiana é o fator de escala que transforma aárea elementar dξdη no quadrado paramétrico emárea elementar correspondente no quadrilátero doplano cartesiano dA, como indicado a seguir.

(4.102)

(4.103)

(4.104)

A matriz de compatibilidade cinemática é agora dadapor:

E a matriz de rigidez:

sendo t a espessura do elemento. Aintegração da matriz de rigidez é feita no planoparamétrico por integração numérica pelo Método deGauss. Se forem usados ng pontos de Gauss comcoordenadas paramétricas e e pesos deintegração e , a expressão (4.103) pode serreescrita como:

(4.105)

(4.106)

4.4 Elementos da família deLagrangeA família dos elementos lagrangeanos é assimdenominada porque as funções de interpolação desseselementos podem ser facilmente geradas por produtosde polinômios de Lagrange. Polinômios de Lagrange deprimeiro grau podem ser obtidos com o uso dascoordenadas dos pontos notáveis ξ0 = −1 e ξ1 = 1 noeixo ξ, como indicado a seguir:

O polinômio de Lagrange Lξi(ξ) tem a propriedadede valer 1 na coordenada ξi e 0 nas coordenadas ξj,para j ≠ i. A Figura 4.9 esclarece.

(4.107)

(4.108)

FIGURA 4.9 Polinômio de Lagrange Lξ0(ξ) para 2pontos.

Analogamente, polinômios de primeiro grau deLagrange Lη0(η) e Lη1(η) podem ser gerados com o usodas coordenadas dos pontos notáveis η0 = −1 e η1 = 1 noeixo η como indicado a seguir:

Uma forma simples e sistemática de obter funções

(4.109)

(4.110)

(4.111)

de interpolação, também denominadas de funções deforma, para elementos finitos pode ser por meio deprodutos de polinômios lagrangeanos. A função deinterpolação N1(ξ,η) no plano paramétrico para oelemento quadrilátero representado na Figura 4.5 podeser obtida por:

Analogamente, as demais funções Ni(ξ,η),i = 2,..4, são dadas por:

(4.113)

(4.114)

(4.115)

(4.112)Observe que as funções de interpolaçãodo elemento de 4 nós da família deLagrange são idênticas às funções de interpolação doelemento de 4 nós da família Serendipity dadas em(4.58). Os elementos são, portanto, iguais eisoparamétricos.

Os polinômios de Lagrange do segundo grau aquidenominados de Lξ0(ξ), Lξ1(ξ) e Lξ2(ξ) são gerados como uso das coordenadas notáveis ξ0 = −1, ξ1 = 0 e ξ2 = 1no eixo ξ como indicado a seguir:

O polinômio de Lagrange Lξi(ξ) tem a propriedadede valer 1 na coordenada ξi e 0 nas coordenadas ξj,para j ≠ i. A Figura 4.10 esclarece.

(4.116)

(4.117)

FIGURA 4.10 Polinômio lagrangeano Lξ0(ξ) para 3pontos.

Analogamente, polinômios de Lagrange Lη0(η), Lη1(η)e Lη2(η) podem ser gerados com o uso dascoordenadas dos pontos notáveis η0 = −1z, η1 = 0 e η2 = 1no eixo ξ como indicado a seguir:

(4.118)

Procedendo de maneira análoga àquela utilizadacom polinômios de Lagrange de primeiro grau, funçõesde interpolação Ni(ξ,η) no plano paramétrico podemser geradas para o elemento “quadrilátero” de 9 nós.O elemento pode modelar lados curvos. O elementopossui 8 nós no contorno, sendo 4 nos “vértices” do“quadrilátero” e 4 sobre os “lados” e um nó centralcomo indicado na Figura 4.11. O nó central é necessárioporque as funções Ni(ξ,η) valem zero em ξ = 0 e η = 0.

(4.119)

FIGURA 4.11 Mapeamento do ponto P(ξ,η) doquadrilátero no plano paramétrico no ponto P(x, y) do“quadrilátero” de lados curvos no plano cartesiano.

A função de forma N1(ξ,η) no plano paramétrico podeser obtida por:

A função N1(ξ,η) que vale 1 no nó 1 e 0 nos demaisnós, está representada na Figura 4.12.

(4.120)

(4.121)

FIGURA 4.12 Função de interpolação N1(ξ,η) doelemento lagrangeano de 9 nós.

Analogamente, as demais funções Ni(ξ,η), i = 2, ..., 9,são dadas por:

(4.122)

(4.123)

(4.124)

(4.125)

(4.126)

(4.127)Uma vez definidas as funções de interpolação

Ni(ξ,η),a formulação dos elementos de 4 ou 9 nós segueo mesmo padrão do elemento isoparamétrico dafamília Serendipity. Sendo nnos , o número de nós doelemento, pode-se escrever:

(4.128)

(4.129)

(4.130)

(4.131)

e,

onde, xi e yi são as coordenadas nodais do nó irelativas aos eixos horizontal x e vertical y,respectivamente; e ui e υi são os deslocamentos nodaisdo nó i re lativos aos eixos x e y, respectivamente.

A expressão (4.129) pode ser reescrita como:

Onde a matriz N(ξ,η) tem a forma com nnos sendo onúmero de nós do elemento.

E o vetor d é definido como:

(4.132)

(4.133)

Usando a matriz jacobiana J(ξ,η) já utilizada noselementos da família Serendipity,

e , generalizando para um elemento de nnos nós,obtém-se:

(4.134)

(4.135)

(4.136)

ou, matricialmente,

onde os subscritos ξ e η significam aderivada em relação a ξ e η, respectivamente.

Sucintamente, a expressão (4.135) pode ser reescritacomo:

Como já visto anteriormente, a inversa da matrizjacobiana Γ(ξ,η) que como J(ξ,η) tem dimensão 2x2, éobtida da operação:

(4.137)

(4.138)

(4.139)

transforma derivadas paramétricas de uma funçãoem derivadas cartesianas da mesma função. Sendoassim, pode-se escrever:

ou, sucintamente,

Onde, u,c é o vetor que contém as derivadascartesianas das componentes de deslocamentos u e υ,u,p o vetor que contém as derivadas paramétricas dascomponentes de deslocamentos u e υ e Γu a matriz 4x4que transforma derivadas paramétricas dosdeslocamentos em derivadas cartesianas dosdeslocamentos.

As expressões (4.129) permitem escrever:

(4.140)

(4.141)

(4.142)

(4.143)

onde, DNd(ξ,η) é uma matriz que contém asderivadas das funções de interpolação em relação àscoordenadas paramétricas com o seguinte aspecto:

Em analogia a (4.102) a matriz B decompatibilidade cinemática pode ser escrita como:

Como já visto anteriormente, o determinante damatriz jacobiana é o fator de escala que transforma aárea elementar dξdη no quadrado paramétrico emárea elementar correspondente no quadrilátero doplano cartesiano dA, como indicado a seguir.

A partir da matriz B e da expressão (4.143) pode-seobter a matriz de rigidez como nos itens anteriores. Aexpressão para a integral numérica também é análoga

(4.144)

(4.145)

às que foram obtidas para o elemento “Serendipidy”:

4.4.1 Condensação Da Matriz De RigidezOs elementos da família Serendipity e lagrangeana de 4nós com funções de interpolação formadas porpolinômios bilineares (produto de 2 funções lineares)são idênticos. Todavia, o elemento lagrangeano cujasfunções de interpolação são formadas por produtos depolinômios quadráticos gera um elemento de 9 nósdiferente do elemento Serendipity de 8 nós. O elementolagrangeano de 9 nós pode, porém, ser transformadonum elemento de 8 nós com a eliminação do nó 9interno. Usando o índice c para indicar os 8 nós docontorno e o índice i para indicar o nó interno 9, asequações de equilíbrio do elemento podem ser escritascomo:

Explicitando dina segunda equação em (4.147),obtém-se:

(4.146)

(4.147)

(4.148)

(4.144)

A substituição de (4.148) na primeira equação em(4.147), fornece:

ou, sucintamente,

Onde, Kr é a matriz de rigidez reduzida e fr o vetordas cargas nodais reduzidas, onde os graus deliberdade referentes ao nó 9 foram excluídos.

A técnica de condensação reduz o número total denós do elemento lagrangeano “quadrático” para 8, masnão o torna igual ao elemento Serendipity de 8 nós.

4.5 Exemplos de problemas deestado plano4.5.1 Estudo De Uma Viga EmBalanço, Influência Do RefinamentoDa Malha Nos Resultados

Solução pela teoria da elasticidadeConsidere a viga espessa em balanço de largura

estreita da Figura 4.13 com uma carga concentrada Paplicada na extremidade livre. A distribuição detensões na viga é dada por:

FIGURA 4.13 Viga em balanço.

(4.150)

(4.152)

(4.153)

(4.151)

(4.154)

Sendo I o momento de inércia da seção transversaldado por:

o deslocamento vertical total da linha elástica daviga v(x) obtido para y = 0 é dado pela soma da parceladevida à flexão, vb(x), com a parcela decorrente docisalhamento, vs(x), logo:

A parcela devida à flexão é dada por:

(4.155)

(4.157)

(4.158)

(4.156)

e a devida ao cisalhamento:

sendo E o módulo de elasticidade longitudinal, αs arelação entre a tensão (ou deformação) decisalhamento no eixo da viga (y = 0) e a tensão (oudeformação) média de cisalhamento da seçãotransversal e G o módulo de elasticidade transversal.

Os valores de G e αs podem ser obtidos em funçãodo coeficiente de Poisson ν. A expressão para o cálculode G é :

Determinações rigorosas do valor de αs pela teoriada elasticidade fornecem para a seção retangular oseguinte valor:

O deslocamento δ no ponto do eixo neutro e na

(4.159)

extremidade livre (x = L e y = 0) pode ser estimado paravigas em balanço por:

Dados do problema:b = 1; h = 3; L = 12; E = 20000; υ = 0,2; P = 1;Resultados da teoria da elasticidadeTensão no ponto A, x = 0,5 e y = 1,5.σx(0,5; 1,5) = –7,667;Deslocamento na extremidade livre δ.δ = 0,0134;A elástica υ(x) com as parcelas correspondentes à

deformação por flexão, υb(x), e à deformação porcisalhamento, υs(x), estão representadas na Figura 4.14.

FIGURA 4.14 Elásticas da viga em balanço.

Resultados das análises por elementos finitosPara as análises da viga em balanço por elementos

finitos, são geradas 12 malhas, 6 malhas comelementos isoparamétricos de 4 nós e 6 com elementosCST. As malhas são numeradas de 1 a 6 para cada tipode elemento. As malhas dos elementosisoparamétricos de 4 nós estão representadas naFigura 4.15.

a) Malhas de elementos isoparamétricos: as

distribuições dos elementos adotadas para asmalhas são dadas na Tabela 4.1, onde neh é onúmero de elementos na direção horizontal, nevo número de elementos na direção vertical e ne onúmero total de elementos.

Tabela 4.1Distribuições dos elementos na malha

b) Malhas de elementos CST: as malhas deelementos CST são geradas dividindo-se cadaquadrado das malhas de elementosisoparamétricos em 2 triângulos como indicadona Figura 4.16. Dessa forma são geradas asmalhas de 1 a 6 com, respectivamente, 6, 24,36,72, 144 e 480 elementos.

FIGURA 4.15 Malhas de elementos isoparamétricos de4 nós.

FIGURA 4.16 Formação da malha de triângulos.

Observações:a) Todos os resultados são obtidos com o programa

MTOOL.b) As tensões referentes às malhas 1 a 5 são

calculadas no interior dos elementos que contémo ponto A (0.5; 1.5).

c) As tensões referentes à malha 6 são calculadascom a média dos valores dados nos pontos nointerior dos elementos na vizinhança do ponto A(0.5; 1.5) porque o nó dos modelos recai sobreesse ponto.

d) As tensões estão representadas no gráfico emvalor absoluto. Nos resultados dos programaselas aparecem com valor negativo por serem decompressão.

Os deslocamentos verticais na extremidade livre

estão representados na Figura 4.17 para as 6 malhasestudadas. Os pontos i de 1 a 6 no eixo horizontalrepresentam as malhas de 1 a 6. Os valores dTi, e dQi,marcados no eixo vertical, representam,respectivamente, os deslocamentos na extremidadelivre para as malhas i de elementos triangulares equadrilaterais. A linha horizontal trace jada representao valor dado pela teoria da elasticidade dodeslocamento vertical na extremidade livre δ.

FIGURA 4.17 Deslocamentos verticais na extremidadelivre para as malhas i.

Analogamente, a tensão σx no ponto A decoordenadas (0.5; 1.5) está representada na Figura 4.18para as 6 malhas estudadas. Os pontos i de 1 a 6 noeixo horizontal representam as malhas de 1 a 6. Osvalores σTi, e σQi, marcados no eixo vertical,representam, respectivamente, as tensões no ponto Apara as malhas i de elementos triangulares equadrilaterais. A linha horizontal trace jada representao valor dado pela teoria da elasticidade para a tensãono ponto A, σe.

FIGURA 4.18 Tensão σx no ponto A(0,5; 1,5) para asmalhas i.

Como era de se esperar, com o refinamento damalha, os resultados convergem para os resultadosprevistos pela teoria da elasticidade. Os elementosisoparamétricos de 4 nós fornecem melhoresresultados do que os elementos CST, apesar dasmalhas i de elementos CST apresentarem o dobro doselementos das correspondentes malhas de elementosisoparamétricos de 4 nós. Isso se deve ao campo dedeslocamentos mais rico do elemento quadrilateral em

relação ao do elemento CST. Mesmo com 240elementos quadrilaterais e com 480 elementostriangulares os resultados não atingem os valoresteóricos.

4.5.2 Exemplo Da Barra Tracionada,Elemento CSTDados do problema:

módulo de elasticidade, E = 20000; coeficiente de Poisson, υ= 0,2; carga P = 10

largura da seção transversal, b = 1; altura da seçãotransversal, h = 1

comprimento da barra, 2L = 4.Dados do modelo de elementos finitos:nnodes = 6 (número de nós); nelem = 4 (número de elementos)

FIGURA 4.19 Malha adotada.

Solução do problema pela resistência de materiaisA = Área da seção Transversal A = h · b A = 1F = Força resultante axial F = 2P F = 20σ = Tensão axial

σ = 20ε = Deformação axial

ε = 0,001δ = Alongamento de barra δ = ε · 2L δ = 0,004Solução por elementos finitosA solução é apresentada em pseudolinguagem de

programação.Leitura ou definição da matriz das coordenadas

nodaisx. Na linha i da matriz x estão as coordenadas xi(primeira coluna) e yi (na segunda coluna) do nó i.

O vetor das forças nodais f e a matriz de rigidezglobal K são inicializados (zerados). A variável gdl é onúmero total de deslocamentos nodais, livres e fixos,que no caso do problema plano é de 2 vezes o númerode nós.

A notação i = j...n, significa que será criado um ciclo(loop), onde a variável i varia de j até n.

gdl = 2 · nnodesi = 1..gdl

fi = 0j = 1..gdlKi,j = 0A matriz constitutiva C é calculada para o estado

plano de tensão.

A matriz de incidência Incé definida. O elemento Incijda matriz contém o número do nó i do elemento j,sendo os nós numerados no sentido anti-horário.

As vinculações ou restrições da estrutura sãodefinidas. A variável ndirres define o número dedeslocamentos restringidos enquanto o vetor nres comndirres e lementos define as direções que devem ser

vinculadas. A cada nó i são associadas duas direções,nomeadamente, direção 2.i-1 e direção 2.i.

A matriz de rigidez local do elemento m é calculada.m = 1..nelenConstantes da matriz de compatibilidade cinemática

Bm do elemento m são calculadas.

sendo am a área do elemento.A matriz Bm do elemento m é definida.

A matriz de rigidez do elemento Ke(m) do elemento mé definida.

Ke(m) = BTm · C · Bm · am · b

Os ponteiros que relacionam os graus de liberdade ido elemento m (i de 1 a 6) com os graus de liberdade doelemento m na malha da estrutura dg i,m são calculados.

Montagem da matriz de rigidez da estrutura KPara i = 1..6 j = 1..6

Os valores não nulos do vetor das forças nodais sãodefinidos.

Introdução dos vínculos na matriz de rigidez com atécnica do número grande:

k = 1..ndirres

Cálculo do vetor dos deslocamentos nodais daestrutura d:

Cálculo do vetor dos deslocamentos nodais dem doelemento m:

m = 1..nelem

Cálculo do vetor de deformações εm e do vetor detensões σm no elemento m.

εm = Bm · demσm = C εmValores de εm e σm obtidos para o elemento 4.

Deve-se observar que os resultados reproduzem osresultados da resistência dos materiais.

ε = 0,001, s = 20 e d = 0,004

4.5.3 Exemplo De Cálculo De Cargas

Equivalentes Nodais, ElementoIsoparamétrico Bilinear.Programa em pseudolinguagem de programação.

Vetor das cargas equivalentes nodais fc para cargasconcentradas no interior do elemento.

Pontos notáveis e pesos de Gauss:

Coordenadas nodais:

FIGURA 4.20 Dimensões do elemento.

Matriz de coordenadas nodais:

Espessura: t = 1Funções de interpolação e suas derivadas:

Matriz de interpolação N(ξ,η) com as funções deinterpolação:

Matriz com as derivadas paramétricas das funçõesde interpolação DNx(ξ,η) para obtenção da matrizjacobiana J(ξ,η):

Matriz jacobiana J(ξ,η) e inversa da jacobiana Γ(ξ,η):

Cálculo do determinante de J(ξ,η):det J(ξ,η) = |J(ξ,η)|Vetor das cargas equivalentes nodais fc re lativas às

cargas concentradas fcx e fcy aplicadas no ponto P decoordenadas ξP = 1/2 e ηP = 2/3 no interior do elemento.

FIGURA 4.21 Cargas concentradas no ponto P(ξP,ηP).

Como se pode verificar a seguir, as cargasequivalentes nodais formam um sistemaestaticamente equivalente às cargas aplicadas:

Somatório das forças horizontais: fc1 + fc3 + fc5 + fc7 =20

Somatório das forças verticais: fc2 + fc4 + fc6 + fc8 = 30Vetor das cargas equivalentes nodais fq relativas às

cargas distribuídas qx e qy no interior do elemento

FIGURA 4.22 Cargas distribuídas na área do elemento.

qx = 2; qy = 3Resultantes:área = 2

Integração de Gauss na área do elemento:

Como se pode verificar a seguir, as cargasequivalentes nodais formam um sistemaestaticamente equivalente às cargas aplicadas:

Somatório das forças horizontais: fq1 + fq3 + fq5 + fq7 =4

Somatório das forças verticais: fq2 + fq4 + fq6 + fq8 = 6Vetor das cargas equivalentes nodais fp relativas às

cargas distribuídas px e py aplicadas ao longo do bordo ξ= 1 do elemento.

FIGURA 4.23 Cargas distribuídas ao longo do bordo ξ =1 do elemento.

px = 2; py = 3Resultantes:L = comprimento do bordo; L = 1

Integração de Gauss ao longo do bordo ξ =1 doelemento:

Como se pode verificar a seguir as cargasequivalentes nodais formam um sistema

estaticamente equivalente às cargas aplicadas:Somatório das forças horizontais: fp1 + fp3 + fp5 + fp7 =

2Somatório das forças verticais: fp2 + fp4 + fp6 + fp8 = 3

4.5.4 Exemplo De Placa Circular Vazada,Elemento Isoparamétrico Bilinear.Problema de sólido de revolução resolvido comoproblema de estado plano

Dados:re = 20; ri = 10; pe = pi = p = 10; E = 20000; υ = 0,2; t = 1

FIGURA 4.24 Placa circular vazada.

Solução pela teoria da elasticidade

Observação: como as componentes de tensões σr esσθ são ortogonais e do mesmo valor (−10) e t rz é nulo, ocirculo de Mohr em cada ponto se reduz a um ponto noeixo σ, logo σx, σy = −10, txy = 0.

Função de deslocamento ur(r) ao longo de r.

FIGURA 4.25 Gráfico da função ur(r).

Solução por elementos finitosMalha de elementos finitos utilizada para a análise

do problema.

FIGURA 4.26 Malha adotada com uso da duplasimetria do problema.

Dados: nnodes = 25; nelem = 16; gdl = 2 · nnodes; ri = 10; re= 20; Δr = (re– ri)/4;

Inicialização:

Geração das coordenadas nodaisA matriz n(i,j) armazena a numeração dos nós. O

primeiro nó gerado é n(1,1) = 1. Quando i varia de 2 a 5com j mantido igual a 1, os outros nós sobre acircunferência com mesmo r são gerados. Os nós commesmo j formam uma linha reta que passa pela origemdo sistema de coordenadas e tem inclinação constanteem relação ao eixo X. Em seguida, com j = 2, os nóscom r = ri + Δr são gerados e assim sucessivamentecom raios iguais para pontos de mesmo j. A matriz r(i,j)armazena os raios dos nós n(i,j). O vetor α(i) armazenaos ângulos em radianos com o eixo X das retas quecontém os pontos com mesmo i.

Finalmente, pode-se escrever:

Geração da matriz das coordenadas nodais X:i = 1..nnodes

Pontos notáveis e pesos de Gauss:

Matriz de incidência:

Definição das funções de interpolação e suasderivadas:

Definição da matriz DNd(ξ,η) e DNx(ξ,η):

Definição da matriz He da matriz constitutiva Cparao problema de estado plano de tensões:

Ciclo para gerar as matrizes dos elementos:m = 1..nelemCoordenadas nodais do elemento m:

Matriz Xe(m) com as coordenadas nodais doelemento m:

Matriz jacobiana J(ξ,η,m) e sua inversa Γ(ξ,η,m) paraelemento m:

Geração da matriz Γu(ξ,η,m).

Matriz de compatibilidade cinemática B(ξ,η,m) doelemento m:

Calculo do determinante de J(ξ,η,m) para o elementom:

Cálculo da matriz de rigidez no ponto decoordenadas ξ,η do elemento m:

Matriz de rigidez do elemento m:

Montagem da matriz de rigidez da estruturaPonteiros:

i = 1..8j = 1..8

Vinculação pela técnica do número grande:ndirres = 10Direções restringidas:

k = 1..ndirres

Vetor das forças nodais equivalentes f:

Nós da face interna:

i = 2..4

Nós da face externa:

i = 2..4

Cálculo do vetor de deslocamentos nodais d daestrutura:

d = K−1·f

Cálculo dos vetores dos deslocamentos nodais de(m)do elemento m:

m = 1..nelem

Cálculo dos vetores de deformação e de tensão noelemento m:

Resultado para o elemento 1 com coordenadasparamétricas ξ = 0 e η = 0:

Deve-se lembrar que a solução por elementos finitosmodelou o problema como um problema de estadoplano de tensão. Desse modo, os resultados fornecidosem termos de deformação e tensão referem-se acomponentes de deformação εx, εy e γxyecorrespondentes componentes de tensão σx, σy e τxy.Todavia, as componentes de tensão σx, σy e τxycorrespondem exatamente às componentes de tensãoσr, σθ e τrθ do problema axissimétrico resolvido pelateoria da elasticidade uma vez que o círculo de Mohrse degenera em um ponto que torna os resultadosindependentes do sistema de coordenadas adotado,tanto para tensões como para deformações. Pode-seobservar que os resultados da solução por elementosfinitos se aproximam bastante dos resultados teóricosdados pela teoria da elasticidade. O problemaresolvido como problema de estado plano não forneceas componentes de deformação e de tensão na direçãoz.

4.5.5 Exemplo De Barra TracionadaModelada Por Elemento De Lagrange De9 NósCoordenadas paramétricas notáveis:

Polinômios de Lagrange:

Representação gráfica do polinômio de LagrangeLξ2(ξ). Observa-se que ele vale 0 nas coordenadasparamétricas−1 e 0 e 1 na coordenada paramétrica +1:

FIGURA 4.27 Polinômio lagrangeano Lξ2(ξ).

FIGURA 4.28 Elemento langrageano de 9 nós noplano paramétrico.

Funções de interpolação do elemento:

Representações gráficas das funções deinterpolação N2(ξ,η) e N9(ξ,η):

FIGURA 4.29 Função N2(ξ,η) e N9(ξ,η).

Coordenadas nodais:x1 = 1; x2 = 5; x3 = 5; x4 = 1; x5 = 3; x6 = 5; x7 = 3; x8 = 1; x9 =

3y1 = 1; y2 = 1; y3 = 3; y4 = 3; y5 = 1; y6 = 2; y7 = 3; y8 = 2; y9 =

2Parâmetros geométricos e mecânicos: espessura t =

1; módulo de elasticidade E = 1.000; coeficiente dePoisson ν = 0,2; carga distribuída p = +10 no bordo ξ = +1.

Matriz das coordenadas nodais X:

FIGURA 4.30 Malha de 1 elemento e dados doproblema.

Pontos notáveis e pesos de Gauss:

Matriz constitutiva C e matriz H do problema deestado plano de tensões:

Derivadas paramétricas das funções deinterpolação:

Matriz DNd(ξ,η):

Matriz DNx(ξ,η):

Matriz jacobiana

Determinante da matriz jacobiana J(ξ,η).

Inversa da matriz jacobiana Γ(ξ,η).

Matriz Γu(ξ,η).

Matriz de compatibilidade cinemática B(ξ,η).

Cálculo da matriz de rigidez no ponto decoordenadas paramétricas ξ,η.

Partição da matriz de rigidez em submatrizescorrespondentes aos graus de liberdade c do contorno(nós de 1 a 8) e graus de liberdade i correspondentesao nó interior 9:

Matriz de rigidez condensada Krr:

Vetor das forças equivalentes nodais f:

Partição e cálculo do vetor das cargas equivalentesnodais condensadas fr:

Introdução dos vínculos na matriz de rigidezcondensada com a técnica do número grande:

Cálculo do vetor dos deslocamentos condensados dr:

Cálculo do vetores de deslocamentos do contorno dce do nó interno di:

Cálculo do vetor de deslocamentos d:

Cálculo das deformações e tensões no ponto decoordenada paramétrica ξ,η:

O resultado exato é σx = 10, σy = τxy = 0:

Demo version limitation, this page not show up.

C A P Í T U L O 6

Sólidos tridimensionais

6.1 IntroduçãoExemplos de sólidos tridimensionais (3D) emEngenharia Civil são: blocos de estaca, sapatas, blocosde fundações de máquinas, etc. As análises de sólidos3D por elementos finitos são, ainda hoje , poucoutilizadas devido a dificuldade na geração da malha.Ultimamente, grandes avanços têm sido feitos com oaparecimento de programas para a geraçãoautomática de malhas tridimensionais.

6.1.1 Equações De CompatibilidadeOs campos de deslocamento de um sólido são u(x,y,z),υ(x,y,z) e w(x,y,z), respectivamente na direção dos eixosx, y e z.

As componentes de deformação são dadas por:

(6.1)

Em notação vetorial:

(6.2)

ou, matricialmente,

(6.3)

(6.4)

ou, sucintamente:

Em (6.4), é o vetor das deformações, L amatriz operadora de derivação e u o vetor dascomponentes de deslocamentos.

6.1.2 Equações ConstitutivasA lei de Hooke e o efeito de Poisson permitem escrever

(6.5)

(6.6)

(6.7)

na forma matricial:

ou, sucintamente,

A relação inversa pode ser expressa por:

sendo C a matriz constitutiva para um materialisotrópico e linear elástico de uma estrutura 3D, dadapor:

(6.8)6.2 Elemento tetraedroO elemento tetraedro de 4 nós representado na Figura6.1 para problemas de sólidos tridimensionais tambémapresenta deformação constante assim como oelemento triangular de 3 nós para o problema deestado plano como será visto adiante.

FIGURA 6.1 Elemento tetraedro e seus graus deliberdade.

Os campos que descrevem os deslocamentos nointerior do elemento são polinômios lineares em x, y ez, ou se ja:

(6.9)

(6.10)

(6.11)

ou, matricialmente,

ou, sucintamente,

A escolha de polinômios lineares de 4 termos com 12coeficientes incógnitos a i pode agora ser justificadapelas 12 condições de contorno seguintes:

(6.12)

que podem ser reescritas usando-se a expressão(6.9) como:

(6.13)

(6.14)

(6.15)

ou, sucintamente,

ou,

onde o vetor d contém os deslocamentos nodais.

(6.16)

(6.17)

Substituindo a expressão (6.15) em (6.11), obtém-se:

ou, ainda,

(6.21)

(6.18)

(6.19)

(6.20)

sendo,

A matriz N(x,y,z) tem a forma:

Observando as equações (6.18) e (6.20) épossível escrever:

(6.22)

(6.23)

(6.24)

onde ui, υi e wi são os deslocamentos nodais relativosaos eixos x, y e z, respectivamente.

O volume do elemento é representado pela variávelVol cuja expressão pode ser obtida pelo determinante:

Usando mais uma vez a expressão geral para amatriz de rigidez de um elemento finito qualquer, ouseja:

onde B é a matriz de compatibilidade cinemática quetransforma deslocamentos nodais em deformações nointerior do elemento, ou se ja:

(6.25)

(6.26)

(6.27)

(6.28)

C é a matriz constitutiva que transforma o vetor dedeformações e em vetor de tensões s para o materialde comportamento linear elástico (le i de Hooke).

No caso de um problema 3D, ascomponentes do vetor de deformação e são dadaspela expressão (6.3) que está representada a seguir naforma matricial.

onde L é a matriz operadora de derivação.Substituindo-se (6.18) em (6.26) obtém-se:

A expressão (6.27) pode ser reescrita como:

o que permite concluir que para o elemento emquestão vale ,

(6.31)

(6.29)

(6.30)

sendo a matriz B dada por:

onde as submatrizes Bi representam a parcela deBrelativa ao nó i, dada por:

Os coeficientes bi, ci e di são dados por:

(6.32)

(6.33)

(6.34)

Os nós i, j, k e l seguem a sequência “i j k l i j k l .....”.Assim, se a numeração local dos nós 1, 2, 3 e 4

corresponder aos nós globais 7, 9, 12, 15, B1, quecorresponde ao nó 7 global, será formada com ascoordenadas dos nós j = 9, k=12 e l=15. Já B 2,quecorresponde ao nó 9 global,será formada com ascoordenadas dos nós j=12, k=15 e l=7 e assimsucessivamente. A numeração dos nós do elementodeve seguir a seguinte regra: olhando do nó i para osnós do triângulo oposto, os nós devem ser vistos no

(6.35)

sentido horário para que o volume calculado pelaexpressão (6.22) se ja positivo.

A integração da matriz de rigidez é trivial devido aofato de a matriz B ser constante, ou se ja, independentede x, y e z para o elemento em questão, o que permitereescrever (6.23) como:

6.3 Elemento hexaedroO hexaedro é um elemento da família Serendipity deelementos isoparamétricos. Ele está representado naFigura 6.2.

FIGURA 6.2 Elemento hexaedro.

As coordenadas paramétricas dos nós desseelemento são dadas por:

(6.36)

(6.37)

As funções de interpolação são:

As coordenadas de um ponto no interior doelemento podem ser obtidas por interpolação dascoordenadas nodais:

(6.38)

(6.39)

Assim como os deslocamentos em pontos interiorespodem ser obtidos por interpolação dos deslocamentosnodais:

(6.40)

(6.41)

A matriz jacobiana J( ,h,z) é expressa como:

com

Substituindo as expressões (6.38) em (6.40), chega-sea:

(6.42)

(6.43)

(6.44)

(6.45)

ou, sucintamente,

Derivando-se as expressões (6.39) em relação àscoordenadas paramétricas, chega-se a:

ou, sucintamente,

Sabendo que a matriz ξ(ε,ζ,z) transforma derivadasparamétricas de f em derivadas cartesianas, pode-seescrever:

(6.46)

(6.47)

Onde a matriz 0 é uma matriz com valoresnulos e de dimensão 3x3. Sucintamente, (6.46) pode serreescrita como:

onde u,c é o vetor que contém as derivadascartesianas das componentes de deslocamentos u, υ ew, u,p o vetor que contém as derivadas paramétricasdas componentes de deslocamentos u, υ e w e Γu(ξ,ε,ζ)a matriz que transforma derivadas paramétricas dosdeslocamentos em derivadas cartesianas dosdeslocamentos.

É possível demonstrar também que o determinanteda matriz jacobiana é o fator de escala que transforma

(6.48)

(6.49)

o volume elementar dξ dε dζ no espaço paramétricoem volume elementar correspondente no espaçocartesiano dV = dx dy dz, como indicado a seguir.

As componentes de deformação em um problematridimensional expressas em (6.2) podem ser escritasalternativamente como:

ou, sucintamente,

(6.50)

(6.51)

(6.52)

(6.53)

Usando (6.47) e (6.45), a expressão (6.50) pode serreescrita como:

o que permite concluir que para o elementoem questão, a matriz B vale ,

A obtenção da matriz de rigidez se dá por integraçãonumérica no espaço paramétrico pelo método deGauss. Se forem usados ng pontos de Gauss comcoordenadas paramétricas pesos de integração amatriz de rigidez pode ser escrita como:

6.4 Exemplo de barra tracionadamodelada com sólidotridimensional, elementohexaedro

FIGURA 6.3 Malha de 1 elemento para o problema.

Dados:P = 10; υ = 0,2 ; E = 20000; a = 2(lado do cubo).Funções de interpolação trilineares e suas derivadas:

Matriz DNx(ξ,ε,ζ):

Matriz jacobiana J(ξ,ε,ζ) sua inversa Γ(ξ,ε,ζ):

Γ(ξ,ε,ζ) = DNx(ξ,ε,ζ)−1

Determinante da matriz jacobiana:detJ(ξ,ε,ζ) = |J(ξ,ε,ζ)|Matriz H e matriz 0:

Matriz Γu(ξ,ε,ζ):

Matriz DNdd(ξ,ε,ζ):

Matriz de compatibilidade cinemática B(ξ,ε,ζ):B(ξ,ε,ζ) = H·Γu(ξ,ε,ζ):Matriz constitutiva C:

Matriz de rigidez no ponto P(ξ,ε,ζ):KP(ξ,ε,ζ) = B(ξ,ε,ζ)τ·C·B(ξ,ε,ζ)·detJ(ξ,ε,ζ)Pontos notáveis e pesos para integração de Gauss:npg = número de pontos de Gaussnpg = 8 (opção de 8 pontos de Gauss = 2 x 2 x 2)

Volume do elemento:Vol = 8Vetor das cargas nodais f:i = 1..24fi = 0Forças P na direção do eixo y nos nós 3, 4, 7 e 8:f8 = P; f11= P; f20 = P; f23 = PRestrições:i= 1..6Ki,i = 106 · Ki,i

K14,14 = 106 · K14,14; K17,17 = 106 · K17,17;Observe que no loop, os 3 deslocamentos nas

direções x,y e z dos nós 1 e 2 e os deslocamentos nadireção y dos nós 5 e 6 estão sendo restringidos.

Cálculo do vetor de dos deslocamentos d:d = K−1 · fCálculo das deformações e tensões no ponto P(ξ,ε,ζ):

Resultados no ponto P(0,0,0):

Solução da resistência dos materiais:

Solução do modelo em elementos finitos:σ(0,0,0)2 = 10; ε(0,0,0)2 = 4,979 ×10−4; dε = 9,883×10−4

Demo version limitation, this page not show up.

C A P Í T U L O 8

Análise de estabilidade

8.1 IntroduçãoA seleção de pilares é muitas vezes a parte crucial deum projeto de uma estrutura porque qualquer falhapode ocasionar efeitos catastróficos.

Pilares “esbeltos” podem falhar por “flambagem”elástica, isto é , por deslocamento lateral excessivo comcomportamento linear do material.

Esforços axiais influenciam significativamente osdeslocamentos laterais em pilares assim como forçasde compressão podem produzir deslocamentostransversais indesejáveis em chapas e cascas. Forçasde tração podem diminuir esses deslocamentos, eforças de compressão tendem a aumentá-los oumesmo induzi-los.

Para que se possa avaliar o efeito das cargas axiaisem pilares é preciso realizar uma análise não lineargeométrica com equações de equilíbrio escritas naconfiguração deformada.

Nessa análise , supõe-se que os deslocamentos

laterais são pequenos o suficiente para validar aobtenção das equações de equilíbrio na configuraçãoindeformada.

A hipótese de grandes deslocamentos e pequenasdeformações tem sido suficiente para avaliaçõesprecisas da carga crítica por flambagem.

A flambagem de chapas pode ocorrer em almas oumesas de perfis metálicos. Cilindros de seçãotransversal circular vazada de parede finacomprimidos axialmente também devem seranalisados quanto à instabilidade por flambagemelástica.

8.2 Obtenção da carga críticaem pilares via solução dasequações diferenciais8.2.1 Carga Crítica No Pilar Ideal(Engaste – Extremidade Livre OuPilar Em Balanço Com CargaCentrada)Um pilar engastado na base e com a extremidade dotopo livre está representado na Figura 8.1. A carga decompressão aplicada no topo do pilar é P, a rigidez àflexão da seção transversal é EI, sendo E o módulo deelasticidade longitudinal e I o momento inércia à flexão,

e o seu comprimento L. A Figura 8.1.a representa aconfiguração deformada do pilar quando a cargavertical P é menor do que a carga crítica Pcr. A Figura8.1.b representa a sua configuração deformada quandoa carga P atinge a carga crítica. Quando isso acontecediz-se que ocorreu a flambagem do pilar. A deformadado pilar é representada por v(x). O deslocamentohorizontal na extremidade livre é d. Para calcular acarga crítica em pilares, a equação de equilíbrio deveser escrita na configuração deformada. Da resistênciados materiais sabe-se que, para uma seção transversaldistando x da base do pilar o momento interno, Mint(x) édado por:

FIGURA 8.1 Pilar engastado na base e livre no topo(pilar ideal).

(8.5)

(8.1)

(8.2)

(8.3)

(8.4)

onde υ″(x) representa a curvatura da seçãocalculada pela derivada segunda de υ(x) em relação ax.

O momento externo na mesma seção x é dado nessecaso por:

Para que haja equilíbrio em todas as seções:

logo,

Fazendo-se:

a expressão (8.4) pode ser reescrita como:

(8.8)

(8.6)

(8.7)

(8.9)

(8.10)

A solução da equação diferencial ordinária (8.6) édada pela soma da solução homogênea vH(x) com asolução particular υp(x), ou se ja:

onde,

Assim, a solução total vale :

A primeira derivada de v(x) é expressa por:

Cujas condições de contorno são:

(8.11)

(8.12)

(8.13)

(8.14)

(8.15)

(8.16)

Aplicando as condições de contorno em (8.10), chega-se a:

Introduzindo as constantes C1 e C2 em (8.9),obtém-se:

A expressão (8.15) representa o modo de flambagemda coluna, ou se ja, a forma com que ela flamba. Para x= L, υ(L) = δ, logo,

A equação (8.16) permite duas soluções como

(8.17)

(8.18)

indicado a seguir:

A primeira das duas possibilidades mostradascaracteriza uma situação de repouso ou estabilidade,pois δ ≠ 0. Essa solução não fornece nenhumainformação quanto à carga crítica de flambagem. Asegunda produz uma situação de flambagem ouinstabilidade, já que δ ≠ 0 e , portanto, há umdeslocamento lateral indeterminado da extremidadelivre. Essa solução informa sobre a carga crítica deflambagem, pois:

Quando P atinge o valor dado em (8.18), acarga é denominada de carga crítica Pcr por ser a cargaque produz a flambagem ou instabilidade da coluna.Como o valor de δ é indeterminado para P = Pcr, a curvaP x d é dada em azul na Figura 8.2. A interpretaçãofísica dessa curva é que nenhum deslocamento lateral

ocorre com P ≤ Pcr , ou se ja, δ = 0 nesse caso. Todavia,quando P = Pcr, d se torna indeterminado e a colunaflamba.

FIGURA 8.2 Relação P x δ para a curvatura dada por(8.19).

Se a solução desse problema tivesse sido obtida pelaexpressão mais precisa da curvatura da seção κ.

(8.19)

(8.20)

(8.21)

e não

como foi usado anteriormente em (8.1), a relação P xδ seria representada pela curva vermelha da Figura8.2.

8.2.2 Fórmula Geral Para Carga CríticaEm PilaresProcedendo de modo análogo ao apresentado no itemanterior, cargas críticas em pilares podem ser obtidaspara diversos tipos de condições de contorno. Umafórmula geral interessante que pode ser aplicada auma variedade de pilares é dada a seguir:

onde K é o fator de comprimento efetivo e KL = Le o

comprimento efetivo (ou de flambagem) do pilar. ATabela 8.1 apresenta vários valores de K para diversostipos de condições de contorno em pilares. A Figura 8.3mostra pilares com diferentes condições de contorno,seus respectivos modos de flambagem ecomprimentos efetivos Le. É interessante observar queos comprimentos efetivos representam distânciasentre seções de curvatura ou momento nulo do modode flambagem. No pilar ideal, a figura foi espelhadapara mostrar a distância entre as seções real e virtualde curvatura nula.

Tabela 8.1Valores de K para diversos tipos de condições decontorno

FIGURA 8.3 Comprimentos de flambagem Le parapilares com diferentes condições de contorno.

Observação: as soluções obtidas por meio dasequações diferenciais são importantes por váriosaspectos:

a) Permitem uma compreensão conceitual doproblema.

b) São úteis nos cursos de engenharia comoprimeiro contato com o problema.

c) Fornecem soluções que são benchmarks a serematingidos por outros métodos.

(8.23)

(8.22)

Todavia, a restrição a esse método reside na suacapacidade limitada de resolver problemas maiscomplexos em termos de cargas e condições decontorno.

8.2.3 Tensões Críticas Em PilaresUma vez obtida a carga crítica de um pilar, é possívelcalcular a tensão crítica definida como:

sendo r o raio de giração da seçãotransversal e l a esbeltez do pilar, dados por:

O conceito de tensão crítica introduz o parâmetro deesbeltez λ, tão importante como medida dasensibilidade do pilar à carga crítica.

A fórmula mencionada é chamada função de Euler egraficamente representa a curva de Euler como

indicado na Figura 8.4.

FIGURA 8.4 Tensões críticas em pilares em função daesbeltez l.

A Figura 8.4 é bastante esclarecedora quanto aospossíveis modos de colapso de um pilar. Pilares comíndice de esbeltez elevados, λ ≥ λlim, atingem o colapsopor flambagem elástica quando a tensão atuanteatinge a tensão crítica da curva de Euler antes datensão resistente. Por outro lado, pilares curtos, λ < λlim,têm colapso plástico, pois a tensão atuante atinge a

tensão resistente ao escoamento ou esmagamentoantes da tensão crítica de Euler.

8.3 Método aproximado deRayleigh-Ritz para cálculo dacarga crítica em pilaresComo visto no Capítulo 3, o Método de Rayleigh-Ritzusa funções aproximadoras para as deformadas paraobter soluções que se aproximam das soluçõesanalíticas quando refinadas, ou se ja, quandopolinômios de grau mais elevado ou sériestrigonométricas com mais termos são usados comofunções aproximadoras. Foi visto também que ométodo pode ser formulado a partir de princípios deenergia. Inicialmente, será deduzida a expressão dodeslocamento axial e lementar dΔ (na direção do eixo x)relativo a um comprimento dx da coluna para umdeslocamento lateral dv (na direção do eixo y) daextremidade superior do trecho dx como representadona Figura 8.5. A deformação axial da coluna devida àcarga P será desprezada.

(8.24)

(8.25)

FIGURA 8.5 Relação entre dx, dv, dΔ e v,x.

ou

A parcela dΔ2 pode ser desprezada em (8.25) porqueo incremento dΔ tem uma ordem de grandeza muitoinferior a dx e dυ, o que permite escrever:

(8.26)

(8.27)

(8.28)

(8.29)

O deslocamento vertical Δ na extremidade livre dopilar devido aos deslocamentos verticais incrementaisdΔ pode ser obtido por integração, ou se ja:

A energia potencial total do pilar é dada por PE = U +Wp, isto é ,

ou

(8.31)

(8.30)O princípio da mínima energia potencialtotal estabelece que se a estrutura estiver emequilíbrio estável, υ(x) minimiza o funcional PE(υ(x)).Posto dessa forma, o problema é de cálculo variacional.Usando uma função aproximadora para representarυ(x), o problema passa a ser como encontrar o mínimode uma função (método de Rayleigh-Ritz).

A função aproximadora deve satisfazer as condiçõesde contorno em deslocamento até a ordem dederivação “n-1”, sendo n a maior ordem de derivaçãoque aparece em PE(υ(x)).

8.3.1 Exemplo 1 Do Método De Rayleigh-RitzSeja obter uma estimativa da carga crítica para o pilarideal representado na Figura 8.1 pelo método deRayleigh-Ritz. A função aproximadora adotada pararepresentar o modo de flambagem da coluna é :

A função satisfaz as condições de contorno até a

(8.32)

(8.33)

(8.34)

(8.35)

ordem n-1 = 1, ou se ja, deslocamento transversal erotação. d é o parâmetro incógnito.

As derivadas primeira e segunda de υ(x) sãorespectivamente:

Substituindo as derivada de υ(x) em PE(υ(x)) eintegrando-as, chega-se a:

Aplicando-se a condição de mínimo, obtém-se:

ou

(8.36)

(8.37)

A equação (8.36) tem duas soluções possíveis:

A segunda solução corresponde, fisicamente, a umasituação de flambagem da coluna, pois produzdeslocamento lateral. Logo, a estimativa para Pcr pelométodo de Rayleigh-Ritz para a função aproximadoradada em (8.33) é :

Vale observar que, como a coluna flamba com P = Pcr,a possibilidade de ter

(8.38)

(8.39)

(8.40)

Não tem interesse físico.Como vimos no item 8.2.1., a solução exata desse

problema é:

A aproximação obtida representa um erro de 20%, oque é considerado muito alto.

8.3.2 Exemplo 2 Do Método De Rayleigh-RitzA nova função aproximadora adotada é :

Sendo δ, de novo, o deslocamento lateral daextremidade livre. A função satisfaz as condições decontorno do problema.

Repetindo o procedimento anterior, chega-se a:

(8.41)

(8.42)

Aplicando-se a condição de mínimo, obtém-se:

o que significa um erro de 1,2%, que pode agora serconsiderado satisfatório.

Observação: o método de Rayleigh-Ritz, além desimples, permite o tratamento de vários casos muitocomplexos de serem tratados via solução da equaçãodiferencial, tais como pilar com inércia variável, comdescontinuidades de inércia, com cargas diversas,dentre outros.

A grande limitação do método é , todavia, a escolhade uma função aproximadora adequada, capaz decobrir todo o domínio da estrutura. A solução para esseproblema foi obtida com o MEF, como visto no Capítulo3.

8.4 MEF para o cálculo da cargacrítica em pilaresPara um pilar de pórtico plano, a deformação ε(x),

(8.44)

(8.43)

considerando o alongamento axial devido à flexão dabarra, é dada por:

O termo do meio da expressão (8.43) tem osignificado da deformação produzida em um segmentode barra de comprimento dx devido a umdeslocamento transversal dv, como ilustrado na Figura8.5. Assim, a parcela de ε(x) em questão vale :

Na expressão (8.43), u(x) é a função que descreve odeslocamento axial e υ(x) a função que descreve odeslocamento transversal do elemento. Descrevendou(x) e υ(x) em função dos deslocamentos nodais doelemento, como representado na Figura 8.6, vem:

(8.47)

(8.45)

(8.46)

FIGURA 8.6 Elemento finito de um elemento depórtico plano.

sendo

(8.48)

(8.49)

(8.50)

(8.51)

(8.52)

onde ϕi (x) são as funções de interpolação para osdeslocamentos nodais.

A energia de deformação U de uma viga é dada por:

(8.53)

(8.56)

(8.55)

Substituindo a expressão dada em (8.43) para ε (x)em (8.53), observa-se que:

onde N é a força axial, positiva na tração, e ostermos υ,x

4 são desprezados por serem pequenos emcomparação com os demais. Assim, chega-se a:

Substituindo agora (8.45) e (8.46) naexpressão (8.55) e manipulando-se as equações,obtém-se:

onde Ke é a matriz de rigidez elástica convencionaldo elemento de, formada a partir do primeiro e terceiro

(8.57)

termos de U, e Kg é a chamada matriz de rigidezgeométrica formada a partir do segundo termo de U,ou se ja:

(8.59)

(8.58)

(8.60)

Com a aplicação do teorema deCastigliano, vem

sendo di um deslocamento nodal e fi a força externarelativa à direção de di, chega-se a:

A expressão (8.60) fornece o sistema de equações de

equilíbrio para uma barra. Para a solução de um pórticoqualquer, a matriz de rigidez global do pórtico deve serformada a partir da contribuição apropriada dasmatrizes de cada barra.

O sistema de equações de equilíbrio obtido para opórtico é não linear, pois a matriz de rigidez geométricaKg depende do esforço axial na barra N que, por suavez, depende dos deslocamentos axiais naextremidade da barra, ou se ja, Kg(d). A solução dosistema deve ser obtida por métodos apropriados paraa solução de sistemas de equações não lineares comoo método de substituições sucessivas, o método deNewton-Raphson, o método quase-Newton como oBFGS, dentre outros.

O primeiro passo dessa análise , em qualquer dosmétodos, consiste em uma análise linear elástica dopórtico para se determinar a força normal em cadabarra. Essa análise é executada com a matriz derigidez global representada somente pela matriz derigidez elástica Ke. Para cada barra, pode-se calcular aforça normal atuante N e formar a matriz de rigidezgeométrica Kg. Em uma segunda iteração, a matriz derigidez total da estrutura K seria representada pelasoma das matrizes de rigidez elástica Ke e da matriz derigidez geométrica Kg obtida da primeira iteração. Comas novas matrizes de rigidez, um novo vetor dedeslocamentos é calculado. Esse processo é repetidoiterativamente até a convergência do vetor dosdeslocamentos d.

(8.61)

A expressão (8.60) também pode ser usada para adeterminação do fator de carga crítica l. Esse fatorrepresenta a majoração das cargas nodais f necessáriapara produzir flambagem elástica na estrutura. Para sedeterminar l, é conveniente reescrever a expressão(8.60) como:

Apesar de o vetor de cargas nodais f não estarpresente em (8.61) e le não é dispensado do cálculo de λ> O vetor f é usado numa primeira etapa da análisepara se determinar os esforços normais N em cadabarra. Os esforços normais N serão necessários parase formar as matrizes de rigidez geométrica Kg decada barra e , a partir dessas, a matriz Kgdaestrutura.Se o vetor das cargas nodais f for majorado do fator l,os esforços normais N e conseqüentemente a matrizde rigidez geométrica Kg das barras também devemser majorados proporcionalmente de l. Isso aconteceporque nessa primeira etapa, os esforços normais Nsão determinados por uma análise linear com a matrizde rigidez da estrutura representada somente pelamatriz de rigidez elástica Ke.

Fisicamente, a expressão (8.65) pode serinterpretada da seguinte maneira. Uma formaaproximada de se realizar uma análise não linear éatravés de uma análise linear incremental explícita.

Nesse processo, a matriz de rigidez é atualizada para acarga Δf e , em seguida, um novo incremento de cargaΔf é aplicado à estrutura para o cálculo do novoincremento de deslocamentos Δd. A Figura 8.7esclarece o procedimento para um sistema de um graude liberdade.

FIGURA 8.7 Representação a um grau de liberdade damatriz de rigidez total K para um incremento de carga Δfa partir de uma carga λf.

Para um sistema de n graus de liberdade, o cálculode Δd para um incremento de carga Δf a partir de umacarga Δf pode ser obtido com a expressão (8.60),reescrita como:

(8.62)

(8.63)

Na situação de carga crítica, a Figura 8.7 deve sersubstituída pela Figura 8.8.

FIGURA 8.8 Situação para λ = λcrit.

Nesse caso, a expressão (8.62) seria atualizada para:

A expressão (8.63) coincide com a expressão (8.61).Ela responde à seguinte pergunta: Qual o fator de

carga l que precisa ser aplicado às cargas nodais f paraque a estrutura produza deslocamentos não triviais, Δd≠ 0, mesmo sem incremento nas cargas atuantes, Δf =0? O problema expresso em (8.67) recai em umproblema geral de autovalor cuja solução fornece nautovalores l e n autovetores f que representam Δd),sendo n a dimensão das matrizes Ke e Kg. O menorautovalor calculado é o fator de carga crítica e oautovetor associado ao menor autovalor é o modo deflambagem que representa o modo ou a forma deflambagem da estrutura.

8.5 MEF para cálculo da cargacrítica em placa à flexãoSeja a placa à flexão submetida às forças demembrana (forças que atuam no plano da placa), comoilustrado na Figura 8.9.

FIGURA 8.9 Esforços no plano da placa à flexão.

Para se obter uma matriz de rigidez geométrica paraum elemento de placa à flexão, deve-se proceder deforma análoga ao que foi fe ito para um elemento depórtico plano. Isso significa incorporar na expressão daenergia de deformação U da placa à flexão o trabalhofeito pelas forças de membrana nos deslocamentosproduzidos no plano da placa pelos deslocamentostransversais ao plano médio da placa (deslocamentosverticais na direção do eixo z). Nesse item, apenas aexpressão da matriz geométrica da placa à flexão Kgserá deduzida já que a matriz elástica Ke já foiapresentada no Capítulo 7. O elemento estudado seráo elemento da família Serendipity de 4 nós para a teoriade Mindlin.

Em analogia à expressão (8.48), as deformações no

(8.64)

(8.65)

(8.66)

plano médio da placa, associadas às pequenasrotações w,x e w,y , são dadas por:

sendo w(x,y) o deslocamento transversal nadireção do eixo z. A Figura (8.10) ilustra as diferentescomponentes de deslocamento.

(8.67)

FIGURA 8.10 Movimentos horizontais devidos a w,x ew,y.

Para se obter a expressão da matriz geométrica,será considerada apenas a parcela referente aotrabalho das forças de membrana Nx, Ny e Nxy nasdeformações associadas εx, εy e gxy. Assim:

ou, alternativamente,

(8.69)

(8.68)

(8.70)

No elemento de 4 nós da família Serendipity,o campo de deslocamentos transversais w(x,y) érepresentado no plano paramétrico por:

As derivadas paramétricas de w(ξ, h) são dadas por:

ou, sucintamente,

(8.71)

(8.72)

(8.73)

(8.74)

As derivadas cartesianas de w(ξ, η),x e w(ξ, η),ypodem ser obtidas das derivadas paramétricas de w(ξ,η),x e w(ξ, η),η por meio da pré-multiplicação pela matrizΓ (ξ, η) obtida como indicado a seguir.

Assim,

(8.75)

(8.77)

(8.76)

(8.78)

Logo,

ou,

onde,

e

(8.79)

(8.80)

(8.81)

(8.82)

Substituindo-se (8.77) em (8.68) chega-se a:

sendo,

A integração em (8.81) pode ser feita no planoparamétrico.

A integração em (8.82) pode ser feita

numericamente pelo método de Gauss. Observa-seque a matriz geométrica em (8.82) independe daspropriedades do material.

8.6 Exemplos de análise deestabilidade por elementosfinitos8.6.1 Carga Crítica Em Pilar Ideal

FIGURA 8.11 Pilar ideal estudado com um elemento.

Dados:

Vínculos para pilar ideal.

Solução por elementos finitos com um elemento.Cálculo do vetor dos autovalores l do problema de

autovalor generalizado.

Solução pela resistência dos materiais:

Comparação:Carga crítica por elementos finitos:

Carga crítica pela resistência dos materiais:

Por que Pcrit1 é maior que Pcrit2? Porque o MEF, aoaproximar a deformada, fornece modelos mais rígidosdo que os “exatos”.

Solução com L = 3.L = 3

nnodes = 4 gdl = 2 · nnodesInicialização:

Incidência:

FIGURA 8.12 Malha do pilar ideal com 3 elementos.

Ponteiros:

Matriz de rigidez global:

Vínculos pilar ideal:

Cálculo do vetor dos autovalores l do problema de

autovalor generalizado:

Carga crítica por elementos finitos:

Carga crítica pela resistência dos materiais:

erro = 1,028 × 10−4

Percebe-se que com mais elementos, ou se ja, com orefinamento da malha o erro relativo diminui.

8.6.2 Estudo De Um Pilar BiarticuladoDados:

FIGURA 8.13 Malha do pilar biarticulado com 3

elementos.

dois elementos.

Matriz de rigidez elástica e geométrica do pilar.L = 5;Montagem direta das matrizes KGe e KGg.

Vínculos do pilar ideal.

Cálculo do vetor dos autovalores l do problema deautovalor generalizado.

Carga crítica por elementos finitos.

Carga crítica pela resistência dos materiais.

Erro relativo.

Deslocamento lateral em pilar biarticulado com cargaexcêntrica P = 50

Solução da resistência dos materiais

exc = 0,1 (excentricidade da carga);

Solução por elementos finitos

Matriz de Rigidez Tangente:KGt = KGE + (-P)KGg

Comparação:Solução da Resistência dos Materiais δ = 0,042Solução em elementos finitos: d13 = −0,042Análise linear:

Erro da análise linear:

Deslocamento lateral do pilar com carga muitopróxima à carga critica.

P = 197KGt = KGE + (-P)KGg

8.6.3 Estudo Da Flambagem De UmaPlaca À Flexão, Elemento SerendipityIsoparamétrico BilinearA flambagem de uma placa à flexão será estudada comum modelo de 1 elemento finito. A matriz de rigidezelástica Keda placa à flexão é dada pelo elemento dateoria de Mindlin e a matriz geométrica Kg é a doelemento isoparamétrico bilinear.

FIGURA 8.14 Flambagem de placa à flexão commodelo de 1 elemento.

Dados: Coordenadas paramétricas dos nós do elemento:

Pontos notáveis e pesos de Gauss: w1 = 2 (peso para 1 ponto de integração). w2 = 1 (peso para 2 pontos de integração). Coordenadas nodais e espessura t :

Forças por unidade de comprimento aplicadas noplano médio:

Matriz das coordenadas nodais X:

Matrizes constitutivas Db e Ds:

onde

Funções de interpolação bilineares e suasderivadas.

Matriz DNx(ξ,η).

Matriz jacobiana: J(ξ,η)

Cálculo do determinante de J(ξ,η) no ponto P(ξ,η):

Cálculo da inversa da matriz jacobiana Γ(ξ,η) noponto P(ξ,η):

Cálculo da matriz G(ξ,η) no ponto P(ξ,η):

Usando dois pontos de integração de Gauss para Kbe um ponto para Ks (integração seletiva).

Pontos notáveis e pesos de Gauss: Para integração 1 × 1.

Para integração 2 × 2.

Matriz de rigidez total.

Introdução dos vínculos com a técnica dos númerosgrandes.

Matriz N dos esforços no plano.

Matriz de rigidez geométrica Kgp(ξ, η) no ponto P(ξ,η).

Cálculo das matrizes de rigidez geométrica doelemento.

Cálculo dos dois menores autovalores e doautovetor associado ao menor autovalor do problemade autovalor generalizado.

Coeficiente de Rayleigh.

Demo version limitation, this page not show up.

Demo version limitation, this page not show up.

Demo version limitation, this page not show up.

Referências Bibliográficas1. ABAQUS. Analysis User’s Manual Version 6.4.

Hibbitt, Karlsson; Sorensen. Civil Engineeringfor Practicing and Design Engineers, 4, 2002.

2. ARGYRIS JH, KELSEY S. Energy Theorems andStructural Analysis. London: Butterworths; 1960;(collection of papers published in AircraftEngineering in 1954 and 1955).

3. ASSAN AE. Método dos elementos finitos –primeiros passos. Campinas: Unicamp; 1999.

4. BATHE KJ. Finite Element Procedures in EngineeringAnalysis. New Jersey: Prentice-Hall; 1982.

5. BREBBIA CA, FERRANTE AJ. The Finite ElementTechnique. Edições Urgs 1975.

6. COOK RD, MALKUS DS. Concepts andApplications of Finite Element Analysis, Plesha. 3John Wiley & Sons 1989.

7. COURANT R. Variational Methods for theSolution of Problems of Equilibrium andVibration. Bulletin of the American MathematicalSociety. 4973;49:1–23.

8. ERGATOUDIS I, IRONS BM, ZIENKIEWICZ OC.Curved Isoparametric, Quadrilateral Elements

for Finite Element Analysis. Int Journal SolidsStructures. 1968;4(n. 1):31–42.

9. GERE JM, Gere W, Weaver Jr. W. Análise deestruturas reticuladas. Rio de Janeiro: GuanabaraDois; 1981.

10. HINTON E, OWEN DRJ. Finite ElementProgramming . Academic Press 1980.

11. HUGHES TJR. The Finite Element Method, LinearStatic and Dynamic Element Analysis. Dover 2000.

12. MARTHA LF. Análise de estruturas: conceitos emétodos básicos. São Paulo: Elsevier; 2010.

13. Martin HC, Carey GF. Introduction to FiniteElement Analysis, Theory and Application. McGraw-Hill 1973.

14. MELOSH RJ. Structural Engineering Analysis byFinite Elements. New Jersey: Prentice Hall; 1990.

15. MSC Nastran, Accurate , Efficient & AffordableFinite Element Analysis. MTool – BidimensionalMesh Tool – Manual do Usuário, PUC-Rio, 1992 eMG – Mesh Generation, Manual do Usuário,versão 3.0, PUC-Rio, 1996.

16. NITZ M, GALHA R. Mathcad 12 – guia prático. SãoPaulo: Érica; 2005.

17. Notas de aula, prof. Agustin Ferrante,COPPE/UFRJ, 1974.

18. PRZMIENIECKI JS. Theory of Matrix StructuralAnalysis. New York: McGraw-Hill Book Company;1968.

19. RUBINSTEIN MF. Structural Systems: Statics,Dynamics and Stability. Prentice-Hall 1970.

20. SAP2000 manual v12.21. Sistema de análises de estrutura (SALT) – DME –

UFRJ.22. SOBRINHO AS. Introdução ao método dos elementos

finitos. Rio de Janeiro: Ciência Moderna; 2006.23. SORIANO HL, SOUSA LIMA S. Método de

elementos finitos em análise de estruturas. SãoPaulo: Edusp; 2003.

24. SÜSSEKIND JC. Curso de análise estrutural. SãoPaulo: Globo; 1991; v. I e II.

25. TIMOSHENKO SP, GOODIER JN. Teoria daelasticidade. 3 Rio de Janeiro: Guanabara Dois;1980.

26. TIMOSHENKO S, GERE SP. Mecânica dos sólidos.Rio de Janeiro: LTC; 1994; v. I e II.

27. TIMOSHENKO S, WOINOWSKY-KRIEGER S.Theory of Plates and Shells. 2 New York: McGrawHill; 1959.

28. TURNER MJ, CLOUGH RW, MARTIN HC, ToppLJ. Stiffness and Deflection Analysis of ComplexStructures. Journal of Aeronautical Sciences.1956;23(n. 9):805–823.

29. VENÂNCIO FILHO F. Análise matricial deestruturas: Estática. Rio de Janeiro: AlmeidaNeves; 1975.

30. ZIENKIEWICZ OC. The Finite Element Method inEngineering Science. 2 McGraw Hill 1971.