Esta é uma pré-visualização de arquivo. Entre para ver o arquivo original
UNIVERSIDADE FEDERAL DE UBERLÂNDIA FACULDADE DE ENGENHARIA MECÂNICA Programa de Pós-Graduação em Engenharia Mecânica INCT de Estruturas Inteligentes em Engenharia Laboratório de Mecânica de Estruturas Prof. José Eduardo Tannús Reis Prof. Domingos Alves Rade 2011 MÉTODO DOS ELEMENTOS FINITOS APLICADOS À ENGENHARIA MECÂNICA Introdução ao Método dos Elementos Finitos D.A. Rade 1 CAPÍTULO 1 INTRODUÇÃO AO MÉTODO DOS ELEMENTOS FINITOS 1.1 - Fundamentos do método dos elementos finitos O método dos elementos finitos (MEF) é uma técnica de análise numérica destinada à obtenção de soluções aproximadas de problemas regidos por equações diferenciais. Embora o método tenha sido originalmente desenvolvido para a análise estática de sistemas estruturais, ele tem sido utilizado no estudo de uma grande variedade de problemas de Engenharia, nos domínios da Mecânica dos Sólidos, Mecânica dos Fluidos, Transmissão de Calor e Eletromagnetismo. Devido à sua eficiência e flexibilidade, além de sua adequação à implementação em computadores digitais, o MEF tem hoje uma grande difusão tanto no meio acadêmico como no industrial, estando disponível em grande número de “pacotes” comerciais existentes no mercado (ANSYS®, NASTRAN®, ABAQUS®, SYSTUS®, COMSOL®, etc.). Contudo, deve ser lembrado que a utilização eficaz destes programas e a correta interpretação dos resultados requerem o amplo conhecimento, por parte do Engenheiro, dos fundamentos do MEF. A principal motivação para o uso do MEF reside no fato que, devido à complexidade dos problemas práticos de Engenharia, soluções analíticas em forma fechada tornam-se inviáveis ou mesmo impossíveis. Assim, devemos recorrer a técnicas capazes de fornecer soluções numéricas aproximadas. A título de exemplo, consideremos os problemas de determinação da capacidade de carga de uma placa contendo enrijecedores e entalhes de formas complexas, ou de determinação da concentração de poluentes sob condições atmosféricas não uniformes, ou ainda de caracterização do perfil de velocidades em torno de um aerofólio. Para cada um destes problemas podemos obter, sem grande esforço, as equações governantes e as condições de contorno, utilizando princípios elementares da Física. Contudo, nenhuma solução analítica simples poderá ser obtida quando o problema exibir geometria e/ou condições de contorno complicadas, o que quase sempre ocorre em situações práticas. Para contornar esta dificuldade, uma estratégia possível é a simplificação do problema (em termos de sua geometria e/ou condições de contorno) de modo a viabilizar a construção de um modelo matemático cuja resolução analítica seja possível. Contudo, em grande número de casos (talvez na maioria das vezes), este procedimento tem como conseqüência graves imprecisões nas previsões do modelo. Uma segunda alternativa consiste em preservar a complexidade do modelo e empregar técnicas aproximadas de resolução. Esta segunda estratégia, na qual está inserido o MEF, tem sido cada vez mais viabilizada pela crescente capacidade de processamento dos computadores digitais. Em todo problema formulado em domínios contínuos, as incógnitas do problema, denominadas variáveis de campo (que podem ser grandezas escalares, como temperaturas ou vetoriais, como deslocamentos) podem assumir valores independentes em cada ponto do domínio. Conseqüentemente, o problema tem número infinito de incógnitas, sendo caracterizado como um problema infinito-dimensional. Este tipo de problema é geralmente modelado por equações diferenciais parciais, cuja solução analítica é dada por funções que fornecem os Introdução ao Método dos Elementos Finitos D.A. Rade 2 valores das variáveis de campo em função das coordenadas espaciais para todos os pontos do domínio. O MEF é essencialmente um processo de discretização, que visa transformar um problema infinito-dimensional em um problema finito-dimensional, com número finito de incógnitas. O método consiste em dividir o domínio sobre o qual o problema é estudado em várias regiões interconectadas, denominadas elementos. Cada elemento dispõe de um certo número de pontos (interiores e/ou limítrofes), denominados nós ou pontos nodais. O conjunto de elementos utilizados na discretização é denominado malha. Um exemplo é apresentado na Figura 1.1, que mostra a seção transversal de uma palheta de turbina de geometria complexa, discretizada em elementos de forma triangular, tendo, em cada vértice, um nó. Neste exemplo, o problema em questão poderia ser a determinação da distribuição de temperaturas sobre a seção da palheta, conhecidos o fluxo de calor e as condições de contorno. Uma vez definidos os elementos e seus respectivos nós, no interior de cada elemento são admitidas soluções aproximadas para as variáveis de campo, expressas como funções arbitrárias dos valores que as incógnitas assumem nos nós (valores nodais). Estas funções são denominadas funções de interpolação ou funções de forma. São também impostas condições garantindo a continuidade da solução no nós compartilhados por vários elementos. As incógnitas do problema, denominadas graus de liberdade (g.d.l.), passam a ser os valores das variáveis de campo nos pontos nodais, sendo o número destas incógnitas (agora finito), denominado número de graus de liberdade do modelo. Dependendo da natureza do problema, após a discretização, o modelo matemático regente resulta representado por um número finito de equações diferenciais ordinárias ou de equações algébricas, cuja resolução numérica conduz aos valores das incógnitas nodais. Uma vez determinadas estas incógnitas, os valores das variáveis de campo no interior dos elementos podem ser avaliados empregando as funções de interpolação. Conforme será visto mais adiante, a precisão da solução obtida depende essencialmente do número de elementos e do tipo de funções de forma empregadas na discretização. Sendo satisfeitas algumas condições, admite-se que a solução do problema discretizado convirja para a solução exata do problema contínuo à medida que se aumenta o número de incógnitas nodais. Figura 1.1 – Ilustração da malha de um modelo de elementos finitos Introdução ao Método dos Elementos Finitos D.A. Rade 3 Em comparação com outras técnicas numéricas, as principais vantagens o método dos elementos finitos são as seguintes: • elementos de diferentes formas e tamanhos podem ser associados para discretizar domínios de geometria complexa. • a divisão do contínuo em regiões facilita a modelagem de problemas envolvendo domínios não homogêneos, onde as propriedades físicas variam em função das coordenadas espaciais. • o método pode ser todo formulado matricialmente, facilitando sua implementação computacional. A implementação do MEF pode sempre ser efetuada em etapas sucessivas, de forma estruturada. As principais etapas são as seguintes: 1ª) Discretização do domínio. O primeiro passo é a divisão do domínio em elementos. O tipo e número de elementos a serem utilizados devem ser escolhidos de modo a representar adequadamente a geometria do problema e caracterizar convenientemente as variações da solução ao longo do domínio. Alguns tipos de elementos freqüentemente empregados para a discretização de domínios unidimensionais, bidimensionais e tridimensionais são ilustrados na Figura 1.2. Neste aspecto, deve-se observar que problemas unidimensionais são aqueles definidos em domínios representados por apenas uma coordenada espacial (linhas), ao passo que problemas bidimensionais e tridimensionais são aqueles definidos em domínios representados por duas coordenadas espaciais (superfícies) e três coordenadas espaciais (volumes), respectivamente. Os elementos axissimétricos, mostrados na Figura 1.2, são elementos utilizados para a discretização de problemas tridimensionais caracterizados pela existência de simetria geométrica e de carregamento em relação a um dado eixo. Neste caso, o problema tridimensional pode ser formulado como um problema bidimensional. (a) elementos unidimensionais (b) elementos bidimensionais Introdução ao Método dos Elementos Finitos D.A. Rade 4 (c) elementos tridimensionais (d) elemento axissimétrico Figura 1.2 – Ilustração de diferentes tipos de elementos 2ª) Escolha das funções de interpolação. Nesta etapa são escolhidas as funções de interpolação que representam as variáveis de campo no interior de cada elemento. Freqüentemente, mas nem sempre, funções polinomiais são escolhidas como funções de interpolação, devido à facilidade que oferecem para derivação e integração. Os graus dos polinômios utilizados estão relacionados ao número de incógnitas nodais de cada elemento, devendo também atender a certos requisitos de continuidade das variáveis de campo a serem satisfeitos nos nós e nas fronteiras entre elementos imediatamente vizinhos. 3ª) Construção das matrizes elementares. Uma vez escolhidos o tipo e número de elementos e as funções de interpolação, devemos estabelecer as relações matriciais expressando o comportamento (relações de causa-efeito), em termos de propriedades físicas e geométricas, para cada elemento, individualmente. Em outras palavras, procede-se à formulação em nível elementar. Para tanto, podem ser utilizados os seguintes processos: • processo direto, que é baseado no método da rigidez da análise estrutural, através do qual são obtidas as relações matriciais entre forças e deslocamentos nodais, a partir das relações de equilíbrio de forças e compatibilidade de deslocamentos. Procedimento similar pode ser utilizado na modelagem de problemas unidimensionais de transmissão de calor. Embora o processo direto só seja conveniente no tratamento de problemas mais simples, ele tem a grande vantagem de ser de fácil entendimento, permitindo clara interpretação física do significado das relações matriciais obtidas, interpretação esta que pode ser estendida a problemas mais complexos. • processo variacional, que é baseado no Cálculo Variacional e envolve a busca dos pontos críticos – geralmente pontos de mínimo – de um funcional associado ao problema estudado. De acordo com este processo, as relações matriciais em nível Introdução ao Método dos Elementos Finitos D.A. Rade 5 elementar resultam da imposição da condição de estacionaridade do funcional associado ao problema. Em problemas de Mecânica dos Sólidos, por exemplo, este funcional pode representar a energia de deformação ou a energia potencial complementar. O processo variacional, embora mais complicado que o processo direto sob o ponto de vista teórico, permite estender o MEF a problemas mais complexos. Contudo, sua aplicabilidade é limitada aos problemas regidos por princípios variacionais, que estabelecem a existência de funcionais. • processo dos resíduos ponderados. Este é um procedimento ainda mais versátil que os dois procedimentos anteriores, sendo baseado integralmente em operações matemáticas. O processo opera diretamente sobre as equações diferenciais que governam o problema e prescinde da existência de um funcional ou de um princípio variacional. 4ª) Montagem das matrizes elementares para obtenção das matrizes globais. Para caracterizar o comportamento do sistema completo, resultante da associação dos vários elementos, devemos agrupar as matrizes de cada um dos elementos de uma forma adequada. Em outras palavras, devemos combinar as equações matriciais expressando o comportamento dos elementos individuais para formar as equações matriciais que descrevem o comportamento do sistema em todo o domínio. Este processo é conhecido como montagem das matrizes globais. No processo de montagem, impõe-se a condição que em cada nó onde vários elementos estão interconectados, os valores das variáveis de campo são os mesmos para cada elemento compartilhando aquele nó. No final deste processo, as equações matriciais globais devem ser modificadas para satisfazer as condições de contorno do problema. A ordem das matrizes globais coincide com o número total de incógnitas nodais. Este número é chamado número de graus de liberdade do modelo. 5ª) Imposição dos carregamentos externos e das condições de contorno. As equações matriciais globais devem ser modificadas para satisfazer as condições de contorno do problema, que expressam o fato que alguns valores das incógnitas nodais são prescritos. Assim, por exemplo, em problemas de transferência de calor, os valores da temperatura em alguns pontos do contorno podem ser previamente conhecidos. Da mesma forma, deve-se alterar as equações globais para leva em conta que, em alguns nós, cargas externas conhecidas (forças, fluxos de calor, etc.) são aplicadas. Ao final deste processo, o número total de incógnitas nodais remanescentes define o chamado número de graus de liberdade do modelo. 6ª) Resolução do sistema de equações. Ao final do processo de montagem das matrizes globais, o modelo matemático do problema estará representado por um conjunto de equações, que podem ser lineares ou não lineares, algébricas ou diferenciais, dependendo da natureza do problema enfocado. Estas equações devem então ser resolvidas numericamente para a determinação dos valores das variáveis de campo nos pontos nodais. Neste processo de resolução, procedimentos numéricos apropriados, implementados sob a forma de rotinas computacionais, devem ser utilizados. Introdução ao Método dos Elementos Finitos D.A. Rade 6 7ª) Realização de cálculos complementares. Em várias situações, cálculos complementares devem ser realizados para a determinação de grandezas dependentes das variáveis de campo, determinadas na etapa precedente. Assim, por exemplo, nos problemas de Mecânica dos Sólidos, uma vez determinados os deslocamentos, cálculos adicionais são necessários para a determinação das deformações (utilizando as relações deformação-deslocamento) e das tensões (utilizando as relações tensão-deformação). 1.2 - Domínios de aplicação do MEF Os problemas passíveis de tratamento pelo MEF podem ser divididos em três categorias: • problemas de equilíbrio. Esta é a classe de problemas cuja solução é independente do tempo. São exemplos os problemas da Mecânica dos Sólidos envolvendo a determinação de tensões e deformações em elementos estruturais submetidos a carregamentos estáticos e os problemas da Mecânica dos Fluidos tratando da determinação de distribuições de pressão, velocidade em regime permanente e os problemas de Transferência de Calor em regime permanente. Para este tipo de problema, o processo de discretização através do MEF conduz a um modelo matemático representado por um conjunto de equações algébricas, que podem ser lineares ou não lineares. • problemas de autovalor. Nesta classe de problemas, o modelo matemático obtido é representado por um conjunto de equações lineares homogêneas, caracterizado pela dependência em relação a um parâmetro, cuja resolução conduz a um conjunto de autovalores e autovetores. São exemplos os problemas que tratam da determinação de freqüências naturais e modos de vibração de meios sólidos e fluidos, além de cargas de flambagem de elementos estruturais. No primeiro caso os autovalores correspondem às freqüências naturais e os autovetores associam-se aos modos naturais de vibração; no segundo, os autovalores correspondem às cargas de flambagem e os autovetores dizem respeito aos campos de deslocamentos correspondentes. • problemas de propagação. Os problemas de propagação são aqueles em que se busca caracterizar a evolução das variáveis de campo em função do tempo. É o caso típico de fenômenos que se desenvolvem em regime transitório. Os seguintes exemplos podem ser mencionados: determinação do movimento de sistemas estruturais submetidos a cargas de impacto e determinação de distribuições de temperatura geradas por fluxos de calor variáveis. Alguns exemplos de aplicações práticas do MEF são ilustrados a seguir. Introdução ao Método dos Elementos Finitos D.A. Rade 7 Figura 1.3 - Modelo de EF para análise aeroelástica de um Airbus A-320 (100.000 g.d.l.) (extraído de J.F. Imbert, "Analyse des Structures par Eléments Finis) Introdução ao Método dos Elementos Finitos D.A. Rade 8 (b) Figura 1.4 - Análise térmica por EF de uma tubulação em forma de estrela (de Desai C.S. e Abel J.F., "Introduction to the Finite Element Method") Introdução ao Método dos Elementos Finitos D.A. Rade 9 (c) Figura 1.5 - Análise por EF do escoamento de um fluido ideal em torno de um cilindro posicionado entre duas placas paralelas (de Desai C.S. e Abel J.F., "Introduction to the Finite Element Method") Introdução ao Método dos Elementos Finitos D.A. Rade 10 Figura 1.6 - Análise por EF do resfriamento de um lingote metálico (fonte: http://www.comsol.com (Comsol Group)) Figura 1.7 - Análise dinâmico de freios automotivos pelo MEF (fonte: http://www.simulia.com (Dassault Systems)) Introdução ao Método dos Elementos Finitos D.A. Rade 11 Figura 1.8 – Simulação numérica do processo de soldagem pelo MEF (fonte: http://www.simulia.com (Dassault Systems)) Figura 1.9 – Simulação do comportamento da bexiga humana pelo MEF (fonte: http://www.comsol.com (Comsol Group)) Introdução ao Método dos Elementos Finitos D.A. Rade 12 Figura 1.10 – Caracterização do fluxo magnético em um gerador pelo MEF (fonte: http://www.comsol.com (Comsol Group)) Figura 1.11 – Simulação de impacto em palhetas de uma turbina pelo MEF (fonte: http://www.ansys.com (Ansys, Inc.)) Introdução ao Método dos Elementos Finitos D.A. Rade 13 Figura 1.12 – Simulação de “crash-test” pelo MEF (fonte: http://www.ansys.com (Ansys, Inc.)) 1.3 - Sobre as limitações do método dos elementos finitos Embora o MEF seja reconhecidamente uma ferramenta útil e eficiente para a análise de diversos tipos de problemas de Engenharia, é importante estar atento às limitações do método, lembrando que ele fornece modelos matemáticos aproximados para representar o comportamento de sistemas físicos. Conforme será evidenciado no desenvolvimento do curso, em geral, a elaboração de um modelo de EF envolve a admissão de uma série de simplificações que terão efeito direto sobre a precisão das previsões do modelo. Deve também ser levado em conta que, na maioria das vezes, sob considerações de custo e limitações de recursos computacionais, busca-se estabelecer um compromisso entre a complexidade do modelo e a precisão considerada satisfatória. Algumas fontes de incerteza inerentes à modelagem por EF são: • a não consideração de certos tipos de efeitos físicos, tais como não linearidades, histerese, amortecimento, etc. • erros de discretização, devidos à impossibilidade de se obter uma perfeita representação de domínios de geometria complexa utilizando os tipos de elementos disponíveis. • conhecimento impreciso dos valores de alguns parâmetros físicos e/ou geométricos que são utilizados na elaboração do modelo (ex.: módulo de elasticidade, densidade, condutividade térmica, viscosidade, etc.) • dificuldade de modelar efeitos localizados, tais como junções parafusadas e rebitadas. • erros oriundos do processo de resolução numérica. Introdução ao Método dos Elementos Finitos D.A. Rade 14 Visando a validação de modelos de elementos finitos, um procedimento recomendável é a confrontação de suas previsões com dados provenientes de outros tipos de análises, em particular, de análises experimentais. Várias técnicas foram desenvolvidas, principalmente no âmbito da Dinâmica Estrutural, objetivando a correção sistemática de modelos de elementos finitos a partir da confrontação com dados experimentais (consultar, por exemplo, o livro de Friswell e Mottershead). Outro fator a ser considerado é que a elaboração de modelos de EF de problemas complexos é, na maioria das vezes, um processo interativo, fazendo apelo ao conhecimento do Engenheiro acerca do próprio método e do problema em estudo. 1.4 - Bibliografia • Bathe, K. J., Finite Element Procedures in Engineering Analysis, Prentice-Hall, 1982. • Cook, R.D., Concepts and Applications of Finite Element Analysis, John Wiley & Sons, 1974. • Cook, R., Finite Element Modeling for Analysis, John Wiley & Sons, 1995. • Desai, C.S., Abel, J. F., Introduction to the Finite Element Method, Van Nostrand Reinhold Company, 1972. • Imbert, J. F., Analyse des Structures par Eléments Finis, 3ième édition, Cépauès Editions, 1991. • Friswell, M.I., Mottershead, J.E., Finite Element Model Updating in Structural Dynamics, Kluwer Academic Publishers, 1995. • Moaveni, S., Finite Element Analysis. Theory and Application with ANSYS, Prentice-Hall, 1999. •Zienkiewicz, O.C., The Finite Element Method, 3rd edition, McGraw-Hill Book Co., 1977. •Zienkiewicz, O.C., Morgan, K, Finite Elements and Approximation, John Wiley & Sons, 1983 •Huebner, K.H., Thornton, E.A., The Finite Element Method for Engineers, John Wiley & Sons, 1982 D.A.Rade Formulação do MEF pelo Processo Direto 15 CAPÍTULO 2 FORMULAÇÃO DO MÉTODO DOS ELEMENTOS FINITOS PELO PROCESSO DIRETO No Capítulo 1 havíamos mencionado que as relações matriciais traduzindo o comportamento de cada subdomínio (elemento) podiam ser obtidas empregando um dos três métodos: Método Direto, Método Variacional e Método dos Resíduos Ponderados. Neste capítulo vamos ilustrar, com o auxílio de exemplos simples, a utilização do Processo Direto, que tem a vantagem de proporcionar facilitar a interpretação física para as diversas etapas da elaboração de um modelo de elementos finitos. Embora sua utilização seja viável apenas no tratamento de problemas unidimensionais simples, os conceitos derivados deste método podem ser estendidos, com vantagem, a problemas mais complexos. 2.1 – Análise estática de sistemas compostos por associações de molas lineares. Um dos problemas mais elementares que podem ser examinados sob o enfoque do MEF é a análise estática do sistema constituído por um conjunto de molas lineares, com constantes de rigidez ik , i=1, 2, ..., como aquele exemplificado na Figura 2.1. Nesta figura, if e iu designam, respectivamente, as forças externas aplicadas e os deslocamentos dos nós. O problema consiste em determinar os deslocamentos nodais, dados os valores das forças aplicadas. Figura 2.1 1k 2k 3k nó 1 nó 2 nó 3 nó 4 2f 3f 4f 1f 1u 2u 3u 4u D.A.Rade Formulação do MEF pelo Processo Direto 16 2.1.1 – Obtenção das equações de equilibro em nível elementar Cada mola é identificada com um elemento finito. Para obter a matriz de rigidez para um elemento genérico, este elemento é considerado isoladamente, conforme mostrado na Figura 2.2, onde os índices E e D designam as grandezas associadas aos nós das extremidades esquerda e direita, respectivamente. Figura 2.2 Admitindo que todo o conjunto e, portando, cada elemento, encontra-se em equilíbrio, podemos escrever: EiDi ff −= Sabemos ainda que para as molas lineares, o alongamento é proporcional à força aplicada em suas extremidades, sendo a constante de proporcionalidade o coeficiente de rigidez. Assim, escrevemos: ( )EiDiiEi uukf −−= ( )EiDiiDi uukf −= As duas equações acima podem ser postas na seguinte forma matricial: ⎪⎭ ⎪⎬⎫⎪⎩ ⎪⎨⎧=⎪⎭ ⎪⎬⎫⎪⎩ ⎪⎨⎧⎥⎦ ⎤⎢⎣ ⎡ − − D i E i D i E i ii ii f f u u kk kk (2.1) A equação acima pode ainda ser escrita na forma compacta: E iu Diu ik E if D if D.A.Rade Formulação do MEF pelo Processo Direto 17 ( )[ ] ( ){ } ( ){ }eieiei FUK = i=1,2,... (2.2) onde ( ){ } [ ] ( ){ } [ ]TDiEieiTDiEiei ffFeuuU == são, respectivamente, os vetores de deslocamentos nodais e forças nodais em nível elementar e ( )[ ]eiK é denominada matriz de rigidez elementar. Sobre esta matriz, podemos observar: • o elemento genérico ( )( )mneiK representa a força de reação aplicada no nó m quando provocamos um deslocamento unitário no nó n, mantendo bloqueado o nó m. • é uma matriz simétrica ( )[ ] ( )[ ]Teiei KK = . Levando em conta a interpretação dada acima, a simetria da matriz de rigidez traduz o Princípio de Reciprocidade de Maxwell-Betti, aplicável a sistemas mecânicos lineares. Isto significa que a força de reação que surge no nó m quando provocamos um deslocamento unitário no nó n, mantendo bloqueado o nó m é idêntica à força de reação que surge no nó n quando provocamos um deslocamento unitário no nó m, mantendo bloqueado o nó n. • é uma matriz singular (não inversível). Isto porque, como não foram introduzidas as condições de contorno em nível elementar, a relação (2.2) deve contemplar a existência de uma configuração de equilíbrio sem deformação da mola ( )1+= ii uu e sem o aparecimento de forças nos nós ( )01 == +ii ff . • é uma matriz semi-definida positiva { } ( )[ ] { } { } { }⎟⎠⎞⎜⎝⎛ ≠∀≥ 00 x,xKx eiT 2.1.2 – Obtenção das equações de equilíbrio em nível global Após a obtenção das equações que descrevem o comportamento de cada elemento, isoladamente, devemos considerar o fato que os elementos estão, na realidade, interconectados nos pontos nodais. Fisicamente, a interconexão significa que deve haver, nos nós compartilhados por mais de uma mola, equilíbrio das forças e compatibilidade de deslocamentos. Considerando dois elementos vizinhos, i e i+1, ilustrados na Figura 2.3, estas condições são expressas segundo: 11 ++ =+ iEiDi fff (equilíbrio do nó i+1) (2.3) 11 ++ == iEiDi uuu (compatibilidade de deslocamentos no nó i+1) (2.4) Para imposição destas condições, vamos primeiramente desenvolver a equação matricial (2.1) para os elementos i e i+1: D.A.Rade Formulação do MEF pelo Processo Direto 18 E i D ii E ii fukuk =− (2.5) D i D ii E ii fukuk =+− (2.6) E i D ii E ii fukuk 11111 +++++ =− (2.7) D i D ii E ii fukuk 11111 +++++ =+− (2.8) Somando as equações (2.7) e (2.8) e introduzindo nas três equações resultantes as equações (2.3) e (2.4), obtemos as relações. E iii E ii fukuk =− +1 (2.9) ( ) iDiiiiiEii fukukkuk =−++− ++++ 1111 (2.10) D i D iiii fukuk 1111 ++++ =+− (2.11) Retornando à notação matricial, as equações (2.9) a (2.11) são dispostas sob a forma: ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ = ⎪⎭ ⎪⎬ ⎫ ⎪⎩ ⎪⎨ ⎧ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ − −+− − + + + + ++ ++ D i i E i D i i E i ii iiii ii f f f u u u kk kkkk kk 1 1 1 1 11 11 0 0 (2.12) Figura 2.3 E iu Diu ik E if E iu 1+ Diu 1+ 1+ik E if 1+ Dif 1+ D if nó i nó i+1 nó i+1 nó i+2 D.A.Rade Formulação do MEF pelo Processo Direto 19 Considerando o exemplo ilustrado na Figura 2.1, aplicando sucessivamente o mesmo procedimento aos dois pares de molas que compartilham os nós 2 e 3, obtemos a seguinte relação matricial: ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ = ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − −+− −+− − 4 3 2 1 4 3 2 1 33 3322 2211 11 00 0 0 00 f f f f u u u u kk kkkk kkkk kk (2.13) A equação acima pode ainda ser escrita na forma compacta: ( ) ( ){ } ( ){ }g g gK U F⎡ ⎤ =⎣ ⎦ (2.14) onde ( ){ } [ ] ( ){ } [ ]1 2 3 4 1 2 3 4eT Tg gU u u u u F f f f f= = são, respectivamente, os vetores de deslocamentos e forças em nível global e ( )gK⎡ ⎤⎣ ⎦ é denominada matriz de rigidez global. Sobre esta matriz, podemos fazer as mesmas observações anteriormente apresentadas para as matrizes de rigidez elementares: • o elemento genérico ( )g mn K⎡ ⎤⎣ ⎦ representa a força de reação aplicada no nó m quando provocamos um deslocamento unitário no nó n, mantendo bloqueado o nó m. • é uma matriz simétrica ( ) ( ) Tg gK K⎡ ⎤ ⎡ ⎤=⎣ ⎦ ⎣ ⎦ . • é uma matriz singular. • é uma matriz semi-definida positiva { } ( ) { } { } { }( )0, 0T gx K x x⎡ ⎤ ≥ ∀ ≠⎣ ⎦ O processo de obtenção da matriz de rigidez global a partir das matrizes elementares e denominado montagem da matriz global. Um procedimento mais geral de montagem, bem adaptado à implementação computacional, pode ser formulado através da introdução de matrizes de transformação, compreendendo as seguintes etapas: 1ª) Para cada elemento, estabelecemos relações matriciais expressando as transformações dos deslocamentos nodais em nível elementar com os deslocamentos nodais em nível global.. Assim, no exemplo considerado: D.A.Rade Formulação do MEF pelo Processo Direto 20 [ ]( ) ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ =⎪⎭ ⎪⎬⎫⎪⎩ ⎪⎨⎧ × 4 3 2 1 42 u u u u T u u iD i E i i=1,2,3 2ª) Pré-multiplicamos as matrizes de rigidez elementares por [ ]TiT e a pós- multiplicamos por [ ]iT , para obter as matrizes elementares expandidas, com a mesma dimensão da matriz de rigidez global: ( )[ ]( ) [ ]( ) ( )[ ]( ) [ ]( )42222444 ×××× = ieiTiei TKTK i=1,2,3 (2.15) 3ª) Adicionamos as matrizes elementares expandidas para obter finalmente a matriz de rigidez global. No exemplo: ( ) ( )3 1 g e i i K K = ⎡ ⎤ ⎡ ⎤=⎣ ⎦ ⎣ ⎦∑ (2.16) Detalhamos, a seguir, o procedimento de construção da matriz de rigidez global para o exemplo considerado: • elemento nº 1: ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ ⎥⎦ ⎤⎢⎣ ⎡=⎪⎭ ⎪⎬⎫⎪⎩ ⎪⎨⎧ 4 3 2 1 1 1 0010 0001 u u u u u u D E ( )[ ] [ ] ( )[ ] [ ] ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − − =⎥⎦ ⎤⎢⎣ ⎡⎥⎦ ⎤⎢⎣ ⎡ − − ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ == 0000 0000 00 00 0010 0001 00 00 10 01 11 11 11 11 1111 kk kk kk kk TKTK eTe • elemento nº 2: ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ ⎥⎦ ⎤⎢⎣ ⎡=⎪⎭ ⎪⎬⎫⎪⎩ ⎪⎨⎧ 4 3 2 1 2 2 0100 0010 u u u u u u D E D.A.Rade Formulação do MEF pelo Processo Direto 21 ( )[ ] [ ] ( )[ ] [ ] ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − −=⎥⎦ ⎤⎢⎣ ⎡⎥⎦ ⎤⎢⎣ ⎡ − − ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ == 0000 00 00 0000 0100 0010 00 10 01 00 22 22 22 22 2222 kk kk kk kk TKTK eTe • elemento nº 3: ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ ⎥⎦ ⎤⎢⎣ ⎡=⎪⎭ ⎪⎬⎫⎪⎩ ⎪⎨⎧ 4 3 2 1 3 3 1000 0100 u u u u u u D E ( )[ ] [ ] ( )[ ] [ ] ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − −=⎥⎦ ⎤⎢⎣ ⎡⎥⎦ ⎤⎢⎣ ⎡ − − ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ == 33 3333 33 3333 00 00 0000 0000 1000 0100 10 01 00 00 kk kkkk kk TKTK eTe ( ) ( ) ( ) ( ) 1 2 3 g e e eK K K K⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤= + +⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − −+− −+− − = 33 3322 2211 11 00 0 0 00 kk kkkk kkkk kk (2.17) 2.1.3 – Imposição das condições de contorno As equações (2.13) devem ser modificadas para levar em conta que os valores dos deslocamentos nos nós extremos (1 e 4) podem ser prescritos. Admitamos que as condições de contorno sejam as seguintes: 11 uu = (2.18.a) 44 uu = (2.18.b) onde 1u e 4u são os valores conhecidos. Desenvolvendo (2.13) considerando as relações (2.18), obtemos o seguinte conjunto de equações: D.A.Rade Formulação do MEF pelo Processo Direto 22 ( ) ( ) ( ) ( ) ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ =+− +=++− −=−+ =− ⇒ ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ =+− =−++− =−++ =− 44333 43333222 11232221 12111 44333 34333222 23222111 12111 fukuk ukfukkuk ukfukukk fukuk fukuk fukukkuk fukukkuk fukuk (2.19) Determinamos os deslocamentos nodais incógnitos 2u e 3u resolvendo simultaneamente a segunda e a terceira equações do sistema (2.19): ⎭⎬ ⎫ ⎩⎨ ⎧ − − ⎥⎦ ⎤⎢⎣ ⎡ +− −+= ⎭⎬ ⎫ ⎩⎨ ⎧⇒ ⎭⎬ ⎫ ⎩⎨ ⎧ − −= ⎭⎬ ⎫ ⎩⎨ ⎧⎥⎦ ⎤⎢⎣ ⎡ +− −+ − 433 112 1 322 221 3 2 433 112 3 2 322 221 ukf ukf kkk kkk u u ukf ukf u u kkk kkk (2.20) Os valores das forças de reação nos nós 1 e 4 são finalmente determinados através da primeira e quarta equações do sistema (2.19): ⎪⎩ ⎪⎨ ⎧ +−= −= 43334 21111 ukukf ukukf (2.21) É importante observar que, após a imposição das condições de contorno, a matriz de rigidez modificada, deve ser inversível para que a operação indicada em (2.20) possa ser efetuada. Este será sempre o caso quando as condições de contorno impostas forem suficientes para impedir que o sistema mecânico seja cinematicamente variável ou, em outras palavras, que possa se movimentar sem que haja deformação de pelo menos uma das molas. Um procedimento sistemático para imposição das condições de contorno e cálculo das forças de reação, bem adaptado à implementação computacional, consiste em particionar os graus de liberdade em dois conjuntos: graus de liberdade livres e graus de liberdade impostos. No nosso exemplo, estes dois conjuntos seriam: • graus de liberdade livres: ( ){ } [ ]2 3 TgU u u=A ; ( ){ } [ ]2 3 TgF f f=A • graus de liberdade impostos: ( ){ } [ ]1 4 TgU u u=A ; ( ){ } [ ]1 4 TgF f f=A Após reordenação das equações, o sistema (2.13) pode ser expresso sob a forma: ( ) ( ) ( ) ( ) ( ){ } ( ){ } ( ){ } ( ){ } g gg g i g gg g i ii ii U FK K U FK K ⎡ ⎤ ⎧ ⎫ ⎧ ⎫⎡ ⎤ ⎡ ⎤⎣ ⎦ ⎣ ⎦ ⎪ ⎪ ⎪ ⎪⎢ ⎥ =⎨ ⎬ ⎨ ⎬⎢ ⎥⎡ ⎤ ⎡ ⎤ ⎪ ⎪ ⎪ ⎪⎢ ⎥ ⎩ ⎭ ⎩ ⎭⎣ ⎦ ⎣ ⎦⎣ ⎦ A AAA A A D.A.Rade Formulação do MEF pelo Processo Direto 23 Do sistema acima tiramos duas equações matriciais: ( ) ( ){ } ( ) ( ){ } ( ){ }g g g g gi iK U K U F⎡ ⎤ ⎡ ⎤+ =⎣ ⎦ ⎣ ⎦AA A A A (2.22) ( ) ( ){ } ( ) ( ){ } ( ){ }g g g g gi ii i iK U K U F⎡ ⎤ ⎡ ⎤+ =⎣ ⎦ ⎣ ⎦A A (2.23) De (2.22) obtemos os valores dos deslocamentos nodais correspondentes aos graus de liberdade livres e, em seguida, a partir de (2.23) calculamos as forças de reação correspondentes aos graus de liberdade impostos: ( ){ } ( ) ( ){ } ( ) ( ){ }( )1g g g g gi iU K F K U−⎡ ⎤ ⎡ ⎤= −⎣ ⎦ ⎣ ⎦A AA A A (2.24) ( ){ } ( ) ( ){ } ( ) ( ){ }g g g g gi i ii iF K U K U⎡ ⎤ ⎡ ⎤= +⎣ ⎦ ⎣ ⎦A A (2.25) 2.2 – Exemplo resolvido utilizando o MATLAB® Nesta seção apresentamos a resolução do problema mostrado na Figura 2.1, particularizado para a seguinte situação: 41 20 10k = × N/m 42 30 10k = × N/m 43 50 10k = × N/m 10002 =f N 1 3 0f f= = N 041 == uu A Figura 2.4 mostra o código implementado em ambiente MATLAB® . A solução obtida para o problema é a seguinte: 32 3,64 10 mu −= × 33 2,73 10 mu −= × N37271 ,f −= N72724 ,f −= D.A.Rade Formulação do MEF pelo Processo Direto 24 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Programa para análise estática por EF de sistemas constituídos por molas lineares % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Domingos Alves Rade %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % ENTRADA DE DADOS % valores dos coeficientes de rigidez das molas k_values = 1e4*[ 20 30 10 ]; % valores em N/m % número de nós nb_node = 4; % matriz de conectividade mat_conect=[1 2 ; 2 3 ; 3 4 ] % condições de contorno nos gdl impostos cond_cont=[1 0 ; 4 0]; % forças externas aplicadas nos gld livres forcas_aplic= [ 2 1000 ; 3 0]; % valores em Newtons nb_ele=length(k_values); % CONSTRUÇÃO DAS MAT. ELEMENTARES E MONTAGEM DA MATRIZ DE % RIGIDEZ GLOBAL K_global=zeros(nb_node); for ii=1:nb_ele K_elementar=k_values(ii)*[1 -1 ; -1 1]; mat_ident=eye(nb_node); mat_transf=[mat_ident(mat_conect(ii,1),:); mat_ident(mat_conect(ii,2),:)]; K_global=K_global+mat_transf'*K_elementar*mat_transf; end %IMPOSIÇÃO DAS CONDIÇÕES DE CONTORNO PELO MÉTODO DO %PARTICIONAMENTO DA MATRIZ DE RIGIDEZ % identificação dos gdl livres e gdl impostos gdl_livres=forcas_aplic(:,1); gdl_impostos = cond_cont(:,1); % construção das submatrizes de rigidez K_ll=K_global(gdl_livres,gdl_livres); K_li=K_global(gdl_livres,gdl_impostos); K_ii=K_global(gdl_impostos,gdl_impostos); % construção dos vetores de forças nos gdl livres e de deslocamentos nos gdl % impostos f_liv=forcas_aplic(:,2); d_imp=cond_cont(:,2); % CÁLCULO DOS DESLOCAMENTOS NOS GDL LIVRES E FORÇAS DE REAÇÃO % NOS GDL IMPOSTOS d_liv=inv(K_ll)*(f_liv-K_li*d_imp) f_imp=K_li*d_liv+K_ii*d_imp Figura 2.4 D.A.Rade Formulação do MEF pelo Processo Direto 25 2.3 – Análise estática de treliças planas. A idéia de modelar estruturas como uma série de elementos finitos surgiu como uma extensão dos métodos matriciais utilizados para análise de treliças. Estas estruturas consistem de barras interconectadas em certos pontos considerados como rótulas perfeitas, de modo que apenas forças axiais podem ser transmitidas às barras, não sendo considerados efeitos de flexão e cisalhamento. Um exemplo de treliça é ilustrado na Figura 2.5. Figura 2.5 Na modelagem, cada barra da treliça é tratada com um elemento finito, e os pontos de conexão das barras são associados aos nós do modelo de EF. Para determinar as relações forças-deslocamentos para a treliça completa, devemos, primeiramente, determinar estas relações para um elemento genérico, considerado isoladamente, levando em conta sua orientação, e proceder, em seguida, à montagem das equações matriciais que levam em conta a compatibilidade de deslocamentos e o equilíbrio de forças nos nós. Considerando a Figura 2.5, que mostra um elemento genérico i, orientado segundo um ângulo iα , cujas propriedades relevantes são: comprimento iL , área de seção transversal iA e módulo de elasticidade iE . Vamos utilizar os conceitos elementares da Resistência dos Materiais para obter as relações forças- deslocamentos em nível elementar. Para cada elemento, definimos dois sistemas de referência ilustrados na Figura 2.6: - o sistema de referência global X-Y, comum a todos os elementos. - um sistema de referência local x-y, cuja orientação é determinada pela orientação de cada barra. D.A.Rade Formulação do MEF pelo Processo Direto 26 Figura 2.6 Deve ser notado, na Figura 2.6, que cada um dos nós apresenta duas componentes independentes de deslocamento nas direções X e Y, havendo, portanto, dois graus de liberdade por nó. Além disso, as forças aplicadas nos nós são orientadas na direção da barra. Da condição de equilíbrio do elemento, temos: ( ) ( )E Di ix xF F= − (2.26) A relação força-deslocamento, derivada da lei de Hooke para o elemento considerado permite escrever: ( ) ( ) ( )Ei iD E xi i ix x i i F L L u u E A Δ = − = − (2.27) ( ) ( ) ( )Di iD E xi i ix x i i F L L u u E A Δ = − = (2.28) As duas últimas equações acima podem ser dispostas na seguinte forma matricial. Y X y x iL i α E D ii A,E ( ) ( ),D Di ix xu F ( ) ( ),E Ei ix xu F ( ) ( )XEiXEi F,u ( ) ( )YEiYEi F,u ( ) ( ),D Di iY Yu F ( ) ( ),D Di iX Xu F D.A.Rade Formulação do MEF pelo Processo Direto 27 ( ) ( ) ( ) ( ) i i i i E E i ii i x x D Di i i i i ix x i i E A E A u FL L E A E A u F L L ⎡ ⎤− ⎧ ⎫ ⎧ ⎫⎢ ⎥ ⎪ ⎪ ⎪ ⎪⎢ ⎥ =⎨ ⎬ ⎨ ⎬⎢ ⎥ ⎪ ⎪ ⎪ ⎪−⎢ ⎥ ⎩ ⎭ ⎩ ⎭⎣ ⎦ (2.29) Devemos observar que, na equação acima os deslocamentos e forças nodais estão expressas em termos de suas componentes nas direções dos eixos locais. Como cada um dos elementos possui, em geral, orientação diferente da dos demais, antes da montagem das equações matriciais globais é necessário expressar as relações forças-deslocamentos em termos das componentes no sistema de referência global. Isto é feito introduzindo a seguinte transformação para os deslocamentos, que pode ser facilmente extraída da Figura 2.4: ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) cos sen cos sen 0 0 0 0 cos sencos sen E i X E E E E E i i i i i i ix X Y x Yi i D D D D Di ii i i i i i ix X Y x X D i Y u u u u u u u u u u u u α α α α α αα α ⎧ ⎫⎪ ⎪⎪ ⎪⎧ ⎫= + ⎪ ⎪⎡ ⎤⎪ ⎪ ⎪ ⎪⇒ =⎨ ⎬ ⎨ ⎬⎢ ⎥⎣ ⎦= + ⎪ ⎪ ⎪ ⎪⎩ ⎭ ⎪ ⎪⎪ ⎪⎪ ⎪⎩ ⎭ ou, em forma compacta: ( ) ( ) ( ) ( ) ( ) ( ) i E i X E E i ix Y D D i ix X D i Y u u u T u u u α ⎧ ⎫⎪ ⎪⎪ ⎪⎧ ⎫ ⎪ ⎪⎪ ⎪ ⎪ ⎪⎡ ⎤=⎨ ⎬ ⎨ ⎬⎣ ⎦⎪ ⎪ ⎪ ⎪⎩ ⎭ ⎪ ⎪⎪ ⎪⎪ ⎪⎩ ⎭ (2.30) onde: [ ] ⎥⎦⎤⎢⎣⎡= iiii sencossencosT i ααααα 00 00 (2.31) Introduzindo (2.30) em (2.29) e pré-multiplicando a equação resultante por [ ]T i Tα , obtemos: D.A.Rade Formulação do MEF pelo Processo Direto 28 [ ]T i Tα ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ − − i ii i ii i ii i ii L AE L AE L AE L AE [ ] i Tα ( )( )( )( ) ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ Y D i X D i Y E i X E i u u u u = [ ]T i Tα ( ) ( ) E i x D i x F F ⎧ ⎫⎪ ⎪⎨ ⎬⎪ ⎪⎩ ⎭ (2.32) Com base na Figura 2.4, podemos desenvolver o membro do lado direito de (2.32) sob a forma: [ ]T i Tα ( ) ( ) E i x D i x F F ⎧ ⎫⎪ ⎪⎨ ⎬⎪ ⎪⎩ ⎭ = ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ i i i i cos cos sen cos α α α α 0 0 0 0 ( ) ( ) E i x D i x F F ⎧ ⎫⎪ ⎪⎨ ⎬⎪ ⎪⎩ ⎭ = ( )( )( )( ) ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ = ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ Y D i X D i Y E i X E i i D i i D i i E i i E i F F F F senF cosF senF cosF α α α α (2.33) Assim, a equação (2.32) assume a seguinte forma, representando a relação forças-deslocamentos expressa no sistema de referência global: ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ − −− −− −− iiiiii iiiiii iiiii iiiiii i ii sensencossensencos sencoscossencoscos sensencossensencos sencoscossencoscos L AE αααααα αααααα αααααα αααααα 22 22 22 22 ( )( )( )( ) ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ Y D i X D i Y E i X E i u u u u = ( )( )( )( ) ⎪⎪⎭ ⎪⎪⎬ ⎫ ⎪⎪⎩ ⎪⎪⎨ ⎧ Y D i X D i Y E i X E i F F F F (2.34) ou: ( )[ ] ( ){ } ( ){ }eieiei FUK = i=1,2,... (2.35) onde: ( )[ ] =eiK ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ − −− −− −− iiiiii iiiiii iiiii iiiiii i ii sensencossensencos sencoscossencoscos sensencossensencos sencoscossencoscos L AE αααααα αααααα αααααα αααααα 22 22 22 22 (2.36) é a matriz de rigidez elementar do elemento de treliça. Devemos observar que cada elemento possui quatro graus de liberdade (2 graus de liberdade por nó, correspondendo aos deslocamentos nas direções X e Y. Uma vez desenvolvida a expressão geral da matriz de rigidez para um elemento de barra genérico, podemos utilizar os procedimentos já detalhados na Seção 2.1 para a montagem da matriz de rigidez global – traduzindo a compatibilidade de deslocamentos e equilíbrio de forças nos nós – e a imposição das D.A.Rade Formulação do MEF pelo Processo Direto 29 condições de contorno. No que segue, estas etapas serão detalhadas a partir de um exemplo numérico. 2.4 – Exemplo resolvido usando o MATLAB®. Consideremos a treliça plana constituída de três barras, ilustrada na Figura 2.7(a), cujas características físicas e geométricas são dadas na Tabela 2.1. Interessa-nos determinas os deslocamentos nodais, as reações nos apoios e as tensões em cada uma das barras. A Figura 2.7(b) mostra a numeração escolhida para os nós e elementos, bem como a indicação das componentes das forças e deslocamentos nodais no sistema de referência global adotado. (a) (b) Figura 2.7 45o 100 KN Y X 3 2 1 1 2 3 ( )Xu3 ( )Yu3 ( )Xu2 ( )Yu2 ( )XF3 ( )YF3 ( )YF2 ( )XF2 ( )Xu1 ( )Yu1 ( )YF1 ( )XF1 D.A.Rade Formulação do MEF pelo Processo Direto 30 Tabela 2.1 Elemento A(cm2) E(N/m2) L(m) Conectividade (nós) α (graus) 1 32,3 6,9×1010 2,54 2,3 90 2 38,7 20,7×1010 2,54 1,2 0 3 25,8 20,7×1010 3,59 1,3 135 O programa MATLAB® utilizado e os resultados obtidos são apresentados na Figura 2.5. Deve ser observado que, após o cálculo das forças e deslocamentos nodais, as tensões nas barras foram calculadas de acordo com o seguinte procedimento: Utilizando a equação (2.28), escrevemos: ( ) ( ) ( )Di D Ex ii i ix xi i F E u u A L σ ⎡ ⎤= = −⎢ ⎥⎣ ⎦ Os valores de ( )Di xu e ( )Ei xu são calculados a partir dos deslocamentos nodais expresso no sistema de referência global utilizando a equação (2.30). Os resultados obtidos para este exemplo são os seguintes: • Deslocamentos ( )1 0,6858mmXu = − ( )2 0,9100 mmXu = − ( )2 0,8059 mmYu = − • Reações de apoio ( )2 70,711KNYF = ( )3 70,711KNXF = ( )3 0YF = • Tensões nas barras 1 21,819MPaσ = 2 18,271MPaσ = 3 38,760MPaσ = − D.A.Rade Formulação do MEF pelo Processo Direto 31 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Programa para análise estática por EF de treliças planas %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Domingos Alves Rade Última modificação: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % ENTRADA DE DADOS % propriedades físicas e geométricas das barras [A(m^2) E(N/m^2) L(m) alfa(grau)] propriedades = [ 32.3e-4 6.9e10 2.54 90 38.7e-4 20.7e10 2.54 180 25.8e-4 20.7e10 2.59 135]; % número total de graus de liberdade nb_gdl = 6; % matriz de conectividade (graus de liberdade ordenados segundo: % 1: nó 1, direção x % 2: nó 1, direção y % 3: nó 2, direção x % 4: nó 2, direção y % 5: nó 3, direção x % 6: nó 3, direção y mat_conect=[ 3 4 5 6 1 2 3 4 1 2 5 6]; % condições de contorno nos gdl impostos cond_cont=[2 0 5 0 6 0]; % forças externas aplicadas nos gld livres(valores em Newtons) forcas_aplic= [ 1 0 3 -100000*cos(pi/4) 4 -100000*sin(pi/4)]; [nb_ele,dummy]=size(propriedades); % CONSTRUÇÃO DAS MATRIZES ELEMENTARES E MONTAGEM DA MATRIZ DE RIGIDEZ GLOBAL K_global=zeros(nb_gdl); for ii=1:nb_ele coef=propriedades(ii,1)*propriedades(ii,2)/propriedades(ii,3); alfa=propriedades(ii,4)*pi/180; K_elementar=coef*[1 -1 ; -1 1]; mat_transf=[cos(alfa) sin(alfa) 0 0 0 0 cos(alfa) sin(alfa)]; K_elementar=mat_transf'*K_elementar*mat_transf; mat_ident=eye(nb_gdl); mat_transf1=[mat_ident(mat_conect(ii,1),:); mat_ident(mat_conect(ii,2),:); ... mat_ident(mat_conect(ii,3),:); mat_ident(mat_conect(ii,4),:)]; K_global=K_global+mat_transf1'*K_elementar*mat_transf1; end % IMPOSIÇÃO DAS CONDIÇÕES DE CONTORNO PELO MÉTODO DO PARTICIONAMENTO %identificacao dos gdl livres e gdl impostos gdl_livres=forcas_aplic(:,1); gdl_impostos = cond_cont(:,1); % construção das submatrizes de rigidez K_ll=K_global(gdl_livres,gdl_livres); K_li=K_global(gdl_livres,gdl_impostos); K_ii=K_global(gdl_impostos,gdl_impostos); % construção dos vetores de forças nos gdl livres e de deslocamentos nos gdl impostos f_liv=forcas_aplic(:,2); d_imp=cond_cont(:,2); % CÁLCULO DOS DESLOCAMENTOS NOS GDL LIVRES E FORÇAS DE REAÇÃO NOS GDL IMPOSTOS disp('Deslocamentos nos gdl livres (em milímetros)'); d_liv=inv(K_ll)*(f_liv-K_li*d_imp)*1000 disp('Reacoes de apoio (em kN) '); f_imp=(K_li'*d_liv/1000+K_ii*d_imp)/1000 des_global=zeros(nb_gdl,1); des_global(cond_cont(:,1))=cond_cont(:,2); des_global(forcas_aplic(:,1))=d_liv/1000; % CÁLCULO DAS TENSÕES NAS BARRAS for ii=1:nb_ele alfa=propriedades(ii,4)*pi/180; D.A.Rade Formulação do MEF pelo Processo Direto 32 mat_transf=[cos(alfa) sin(alfa) 0 0 0 0 cos(alfa) sin(alfa)]; des_local=mat_transf*des_global(mat_conect(ii,:)); sigma(ii)=propriedades(ii,2)/propriedades(ii,3)*(des_local(2)-des_local(1)); end disp('Tensões nas barras (em MN/m^2) '); sigma/1e6 Figura 2.8 D.A.Rade Formulação do MEF pelo Processo Direto 33 2.5 - Aplicação do processo direto a outros tipos de problemas unidimensionais A formulação do MEF pelo processo direto, ilustrado na seção 2.1 para sistemas formados por associações de molas lineares pode ser utilizado para a resolução de vários outros tipos de problemas unidimensionais em regime permanente nas áreas de Mecânica dos Sólidos, Mecânica dos Fluidos e Transferência de Calor. Para cada um destes tipos de problemas o processo consiste em determinar, a partir da aplicação das leis Físicas pertinentes, as equações equivalentes às equações de equilíbrio das molas. A partir daí, o procedimento de montagem e resolução das equações globais segue o mesmo procedimento estudado anteriormente. A obtenção das equações equivalentes é ilustrada a seguir. 2.5.1 - Barras em solicitação axial Consideramos aqui problemas formados por barras solicitadas por forças axiais, como exemplificado na Figura 2.9, na qual, para cada segmento da barra, , ,i i iE A l indicam, respectivamente, o módulo de elasticidade, área de seção transversal e comprimento. ei iN U designam, respectivamente, as forças axiais e os deslocamentos aplicados nos nós que delimitam cada segmento. Figura 2.9 1 (1) (2) (3) N1 U1 N2 U2 N3 U3 N4 U4 E1, A1, l1 E2, A2, l2 E3, A3, l3 2 3 4 (i) NiE Ui E Ni D Ui D Ei, Ai, li i i+1 ki= EiAi/li fiE ui E fi D ui D D.A.Rade Formulação do MEF pelo Processo Direto 34 Da Resistência dos Materiais sabemos que para cada segmento em equilíbrio, considerado isoladamente, admitindo comportamento linear elástico, valem as seguintes relações: D D Ei ii i i i i i i i i E AN A E A U U l E D E Di ii i i i i E AN N U U l As duas equações acima por se postas na seguinte forma matricial: 1 1 1 1 E E i ii i D D i i i U NE A l U N (2.38) Comparando as equações (2.1) e (2.38) concluímos que cada segmento da barra solicitada axialmente pode ser considerado como uma mola linear cuja constante de rigidez equivalente é dada por i i i i E Ak l Desta forma, o procedimento de montagem das equações globais e imposição das condições de contorno e carregamentos é feito da mesma forma adotada para os sistemas de molas considerados na Seção 2.1. 2.5.2 Eixos sujeitos a torção Consideramos aqui sistemas formados por eixos solicitados por momentos torsores (torques), conforme mostrado na Figura 2.10, na qual, para cada segmento, , ,i i iG J l indicam, respectivamente, o módulo de cisalhamento, momento polar de inércia da seção transversal e comprimento. Além disso, ei i designam, respectivamente, os torques aplicados e ângulos de torção das seções transversais nos nós que delimitam cada segmento. D.A.Rade Formulação do MEF pelo Processo Direto 35 Figura 2.10 Da Resistência dos Materiais sabemos que para cada segmento em equilíbrio, considerado isoladamente, admitindo comportamento linear elástico, valem as seguintes relações: D D Ei ii i i i G J l E D E Di ii i i i i G J l As duas equações acima por se postas na seguinte forma matricial: 1 1 1 1 E E i ii i D D i i i G J l (2.39) Comparando as equações (2.1) e (2.39) concluímos que cada segmento da barra solicitada em torção pode ser considerado como uma mola linear cuja constante de rigidez equivalente é dada por i i i i G Jk l Para o caso de eixos de seção circular de raio iR , o momento polar de inércia é dado por: 4 2 i i RJ G1, J1, l1 G2, J2, l2 T1 ,1 T2 ,2 T3 , 3 (2) (1) ki = GiJi,/li TiE iE TiD iD Gi, Ji, li TiE , i E TiD (1) , i D D.A.Rade Formulação do MEF pelo Processo Direto 36 O procedimento de montagem das equações globais e imposição das condições de contorno e carregamentos é feito da mesma forma adotada para os sistemas de molas considerados na Seção 2.1. 2.5.3 Problemas de escoamento em regime permanente Nesta seção tratamos o problema de escoamento unidimensional em regime permanente em tubulações que podem ser consideradas constituída pela associação de vários segmentos, como ilustrado na Figura 2.11(a) que mostra a vazão total Q e as vazões que percorrem os segmentos, indicadas por qi , i=1,2, ... O problema consiste em determinar os valores destas vazões, bem como das pressões nos nós que delimitam os segmentos, os quais são também indicados na Figura 2.11. Figura 2.11(a) A Figura 2.11(b) mostra um segmento arbitrário, delimitado por dois nós E e D, nos quais estão indicadas as respectivas pressões Eip e Dip e vazões Eiq e Diq . Figura 2.11(b) Em virtude da hipótese de escoamento em regime permanente, para cada segmento i deve-se ter: A E H D Q Q q1 q2 q2 q3 q4 F G B Cq1 q1 q3 1/Ri E iq D iq E ip D ip i E iq E ip D ip D iq D.A.Rade Formulação do MEF pelo Processo Direto 37 D E i i iq q q (2.41) Além disso, sabe-se da Mecânica dos Fluidos que: D E i i i ip p R q (2.42) onde Ri é o coeficiente de resistência, que depende essencialmente do tipo de fluido, da geometria da seção transversal e do acabamento da superfície interna do duto. Combinando as equações (2.41) e (2,42) escrevemos: 1E E Di i i i q p p R 1D E Di i i i q p p R ou: 1 11 1 1 E E i i D D ii i q p Rq p (2.43) Como nos casos anteriores, a comparação de (2.1) e (2.43) permite concluir que o problema de escoamento unidimensional em regime permanente pode ser modelado da mesma forma que o sistema de molas lineares considerado na Seção 2.1, sendo a rigidez de mola equivalente dada por 1/ Ri. 2.5.4 Transferência de calor unidimensional em regime permanente Tratamos aqui problema de transferência de calor unidimensional em regime permanente em uma barra constituída pela associação de vários segmentos, como ilustrado na Figura 2.12(a), onde Ki, Ai e li designam, respectivamente, o coeficiente de condutividade térmica, a área da seção transversal e o comprimento dos segmentos que compõem a barra. O problema a ser resolvido consiste em determinar a distribuição de temperatura T(x) para um conjunto de condições de contorno impostas nas extremidades 1 e 4. A Figura 2.12(b) mostra um segmento genérico i da barra e as grandezas a ele associadas. D.A.Rade Formulação do MEF pelo Processo Direto 38 Figura 2.12(a) Figura 2.12(b) Recordamos abaixo as equações pertinentes aos dois mecanismos de transferência de calor (condução e convecção) que consideramos estar envolvidos no problema em questão. Condução (Lei de Fourier): i i i dTq K A dx (2.44.a) Convecção: i i i Pq h A T T (2.44.b) qi D Ki, Ai, li EiT E i D iT E iq D iq E D i iT T Ki TiE qiE TiD qiD Q1 T1 Q4 K1, A1, l1 K2, A2, l2 K3, A3, l3 1 3 4 T2 T3 T4 2 x T TP qi D.A.Rade Formulação do MEF pelo Processo Direto 39 Análise em nível elementar E i iq q Di iq q E Di ii i i i K Aq T T l Di iq q Assim, para o elemento i: 1 1 1 1 E E i ii i D D ii i q TK A lq T (2.45) Como nos casos anteriores, a comparação de (2.1) e (2.44) mostra que o problema de transferência de calor unidimensional em regime permanente pode ser modelado da mesma forma que o sistema de molas lineares considerado na Seção 2.1, sendo a rigidez de mola equivalente dada por KiAi/li. Uma especificidade a ser considerada neste problema é que são possíveis combinações de três tipos de condição de contorno aplicadas nas extremidades da barra: 1) valor do calor imposto; 2) valor da temperatura imposta; 3) condição de contorno do tipo convecção com o valor da temperatura T imposto; A condição de contorno do tipo temperatura imposta deve ser tratada da mesma forma que os deslocamentos impostos, apresentada na Seção 2.1.3 e a condição de valor imposto é análoga à aplicação de uma força de valor conhecido nos problemas estruturais. O tratamento das condições de contorno do tipo convecção requer algumas transformações adicionais das equações, uma vez que este tipo de condição de contorno envolve as temperaturas das extremidades da barra, conforme pode ser visto na Equação (2.44.b). Considerando, para exemplificação, a barra discretizada com três elementos finitos, as equações globais resultantes do procedimento de montagem, tomam a forma: D.A.Rade Formulação do MEF pelo Processo Direto 40 1 1 1 1 1 1 1 11 1 1 1 2 2 2 2 1 1 2 2 2 33 3 3 32 2 2 2 4 42 2 3 3 3 3 3 3 3 3 0 0 0 0 0 0 0 0 K A K A l l q TK A K A K A K A l l l l T TK A K AK A K A q Tl l l l K A K A l l (2.46) Admitamos que na extremidade 1 tenhamos uma condição de contorno do tipo temperatura imposta e na extremidade 4 tenhamos uma condição de contorno de convecção, as quais são expressas segundo : 1 1T T 4 4 3 4q h A T T Introduzindo estas expressões em (2.46), obtemos: 1 1 1 1 1 1 1 11 1 1 1 2 2 2 2 1 1 2 2 2 33 3 3 32 2 2 2 4 3 2 2 3 3 4 3 3 3 3 4 3 3 3 0 0 0 0 0 0 0 0 K A K A l l q TK A K A K A K A l l l l T TK A K AK A K A h A T l l l l T K A K A h A l l O sistema de equações acima pode ser decomposto em dois conjuntos de equações cuja resolução fornece o valor das temperaturas desconhecidas T2, T3 e T4 e o calor transferido através da seção transversal no nó 1, segundo: 1 1 2 2 2 2 1 2 2 2 3 3 3 32 2 2 2 3 2 2 3 3 4 3 4 3 3 3 3 4 3 3 3 0 0 0 0 K A K A K A l l l T K A K AK A K A T l l l l h A T T K A K A h A l l D.A.Rade Formulação do MEF pelo Processo Direto 41 1 21 1 1 1 1 1 1 3 4 0 0 T TK A K Aq l l T T D.A. Rade Fundamentos de Cálculo Variacional 42 CAPÍTULO 3 FUNDAMENTOS DE CÁLCULO VARIACIONAL É possível, e muito conveniente, considerar o procedimento de obtenção de modelos de elementos finitos como um processo variacional. Através deste processo procura-se estabelecer as equações do movimento ou de equilíbrio em nível elementar e em nível global mediante a busca dos pontos estacionários (ou pontos críticos, ou pontos extremos) de certas funções, denominadas funcionais, que são funções escalares de funções que definem a distribuição da variável de campo sobre o domínio considerado. Neste contexto o Cálculo Variacional constitui o conjunto de técnicas matemáticas destinadas à determinação dos pontos estacionários de funcionais. Os chamados Princípios Variacionais são leis que estabelecem que a configuração de equilíbrio estático ou de equilíbrio dinâmico de um dado sistema mecânico devem corresponder a pontos estacionários de determinados funcionais associados ao problema. Em diversos problemas da Mecânica do Contínuo, os funcionais representam a energia do sistema. As equações de equilíbrio ou do movimento podem então ser obtidas desenvolvendo, matematicamente, as condições de estacionaridade destes funcionais. Neste Capítulo, apresentaremos os fundamentos do Cálculo Variacional e ilustraremos os desenvolvimentos apresentados através de exemplos baseados no uso de Princípios Variacionais da Mecânica, notadamente o Princípio da Energia Potencial Estacionária e o Princípio de Hamilton. Tais princípios serão apresentados detalhadamente no Capítulo 5. 3.1 Funcionais envolvendo uma variável dependente e uma variável independente Consideremos inicialmente o problema de determinação dos pontos estacionários do seguinte funcional: dxxu,xu,xFuI x x 21 (3.1) Admitiremos ainda que a função u deva satisfazer as seguintes condições de contorno: 11 uxu ; 22 uxu (3.2) Nas equações acima, x designa a variável independente, xu é a variável dependente, que representa a variável de campo associada ao problema, u indica a derivada de xu em relação à variável independente, xu,xu,xF indica genericamente uma relação funcional envolvendo x, xu e xu e 1u e 2u D.A. Rade Fundamentos de Cálculo Variacional 43 designam valores prescritos, assumidos pela variável dependente nos extremos do intervalo em que o problema é definido. Designando doravante por u=u(x) a função que corresponde à solução estacionária exata do funcional definido em (3.1), o procedimento variacional consiste em definir as chamadas funções de comparação, representadas sob a forma: xuxuxu~ 21 xxx (3.3.a) com: xxu (3.3.b) onde 1 é um parâmetro escalar arbitrário e x é uma função contínua arbitrária definida no intervalo 21 xxx , que assume valor nulo nos extremos deste: 01 x 02 x (3.4) Em decorrência de (3.4), as funções de comparação satisfazem as condições de contorno do problema: 11 uxu~ ; 22 uxu~ (3.5) Em (3.3.b), x é entendida como uma função arbitrária que representa uma perturbação infinitesimal acrescentada à solução exata u(x) para cada valor da variável independente x, conforme ilustra a Figura 3.1, onde ficam evidenciadas as condições (3.4) e (3.5). Figura 3.1 x 2x 1x u, u xy xu xu~ 1u 2u D.A. Rade Fundamentos de Cálculo Variacional 44 Na seqüência do processo variacional, devemos introduzir a solução variada (3.3) no funcional dado por (3.1) e impor a condição de estacionariedade do funcional, que determina que sua variação I , associada a qualquer variação arbitrária u deve ser nula. Esquematicamente, este processo é ilustrado como segue: 0 IIIFFuu Buscaremos, primeiramente, expressar a variação da função xu,xu,xF , designada por xu,xu,xF , resultante da substituição de xu por xuxu . Assim, definimos: xu,xu,xFuu,uu,xFxu,xu,xF (3.6) Em decorrência de (3.3.b), temos: xu (3.7) Introduzindo (3.3.b) e (3.7) em (3.6), obtemos: xu,xu,xFxu,xu,xFxu,xu,xF (3.8) Negligenciando os termos de ordem igual e superior à segunda em , procedemos à expansão do primeiro termo do lado direito de (3.8) em série de Taylor, obtendo: x u Fx u Fxu,xu,xFxu,xu,xF (3.9) Introduzindo (3.9) em (3.8), temos: F FF x,u x ,u x x x u u (3.10) A variação do funcional I correspondente à variação da função F é definida segundo: FIFFII = dxxu,xu,xFxu,xu,xFx x 21 dxxu,xu,xFx x 21 dxxu,xu,xFxx 21 (3.11) Introduzindo (3.10) em (3.11), temos: D.A. Rade Fundamentos de Cálculo Variacional 45 2 1 x x F FI x,u x ,u x x x dx u u (3.12) A condição necessária para a estacionaridade do funcional I é que sua variação deve anular-se, ou seja, 0 xu,xu,xI . Assim, com base em (3.12), devemos ter: 2 1 0 x x F Fx x dx u u (3.13) Procedemos, em seguida, à integração por partes da segunda parcela do integrando: 2 1 2 1 2 1 x x x x x x dxx u F dx dx u Fdxx u F (3.14) Introduzindo (3.14) em (3.13), temos: 22 1 1 0 xx x x F F d Fx x dx u u dx u (3.15) Em decorrência de (3.4), a primeira parcela de (3.15) se anula: 02 1 x x x u F , (3.16) de modo que: 2 1 0 x x F d F x dx u dx u (3.17) Neste ponto, evocamos o seguinte lema: Lema: Se 1x e 2x ( 1x < 2x ) são constantes e G(x) é uma função contínua de x no intervalo 21; xx e se 02 1 dxxGx x x para toda e qualquer função continuamente diferenciável x que satisfaz D.A. Rade Fundamentos de Cálculo Variacional 46 01 x 02 x , conclui-se que 0xG identicamente no intervalo 21; xx . Com base neste lema podemos afirmar que a equação (3.17) se verifica se: 0F d F u dx u (3.18) A equação (3.18) é denominada Equação de Euler-Lagrange associada ao funcional dado por (3.1). Por outro lado, a equação (3.16) indica todas as possíveis combinações de condições de contorno impostas nos extremos do intervalo 21; xx . No caso do funcional analisado, as condições de contorno que satisfazem (3.16) são: • Em x= 1x : 01 x (3.18.a) 0 1 xxu F (3.18.b) • Em x= 2x : 02 x (3.18.c) 0 2 xxu F (3.18.d) Neste ponto, introduzimos a seguinte classificação para as condições de contorno: a) condições de contorno geométricas ou de Dirichlet: são aquelas que envolvem diretamente a função de comparação e suas derivadas. b) condições de contorno naturais ou de Neumann: são aquelas representadas por expressões que envolvem as derivadas da função xu,xu,xF em relação a um de seus argumentos ou derivadas de x . De acordo com estas definições, as condições de contorno indicadas pelas equações (3.18.a) e (3.18.c) são condições de contorno geométricas, ao passo que aquelas indicadas pelas equações equações (3.18.b) e (3.18.d) são condições de contorno naturais. D.A. Rade Fundamentos de Cálculo Variacional 47 Observações: • a equação de Euler-Lagrange (3.18) apresenta-se sob a forma de uma equação diferencial, cuja resolução, complementada pelas condições de contorno (3.2), fornece a função u(x) que corresponde ao ponto estacionário do funcional (3.1). • para um determinado ponto estacionário do funcional, pode-se determinar sua natureza (máximo ou mínimo), a partir da análise do sinal da segunda variação do funcional I. A segunda variação é obtida pela expansão de uuI em torno de uI utilizando uma série de Taylor limitada aos termos de segunda ordem em u . Para o funcional em questão, escrevemos: uu uu Fu u Fu u Fu u Fu u Fu,u,xFuu,uu,xF 2 2 2 2 2 2 2 2 1 2 1 ou: u,u,xFu,u,xFu,u,xFuu,uu,xF 2 (3.19) com: u u Fu u Fu,u,xF (3.20) uu uu Fu u Fu u Fu,u,xF 2 2 2 2 2 2 2 2 2 1 2 1 (3.21) Introduzimos agora a segunda variação do funcional I: dxuu,uu,xFuuI x x 21 (3.22) Introduzindo (3.19) em (3.22), temos: 2 1 2 1 2 1 22 x x x x x x dxu,u,xFdxu,u,xFdxu,u,xFuuI (3.23) Levando em conta (3.1) e (3.11), reescrevemos a equação acima sob a forma: IIuIuuI 2 (3.24) onde: D.A. Rade Fundamentos de Cálculo Variacional 48 2 1 222 x x dxu,u,xFI (3.25) é a segunda variação do funcional I. Introduzindo (3.21) em (3.25) e levando em conta (3.3.b), escrevemos: 2 1 2 2 2 2 2 2 2 22 2 2 1 x x dx uu F u F u FI (3.26) Pode-se mostrar que um dado ponto estacionário do funcional I é um ponto de mínimo se 02 I e um ponto de máximo se 02 I para qualquer função arbitrária . Exemplo 1: O problema de equilíbrio de uma corda tensionada Consideremos a corda elástica de comprimento L ilustrada na Figura 3.2, sujeita a uma tração axial constante T e a um carregamento transversal distribuído (força por unidade de comprimento) q(x). Nesta mesma figura u=u(x) indica o campo de deslocamento transversais correspondente à posição de equilíbrio da corda. Pode-se mostrar que a energia potencial total do sistema, incluindo a energia de deformação da corda e o trabalho da carga distribuída é dada pelo seguinte funcional: dxquuTuI L 0 2 2 (a) Figura 3.2 Confrontando a equação (a) com a equação (3.1), reconhecemos: L q(x) x y u(x) D.A. Rade Fundamentos de Cálculo Variacional 49 quuTu,uF 2 2 (b) Após o cálculo das derivadas parciais necessárias, as equações (3.18) e (3.16), nos fornecem, respectivamente, a equação de Euler-Lagrange e os termos de contorno do problema: T xqu (c) 00 LuT (d) Admitindo constanteqxq , integrando duas vezes (c) e introduzindo as condições de contorno determinadas pelos vínculos mostrados na Figura 3.1: 00 u 0Lu , obtemos a seguinte expressão para o campo de deslocamentos transversais da corda em equilíbrio: xLx T qxu 2 Investiguemos agora a natureza do ponto estacionário, mediante a análise do sinal da segunda variação do funcional (a), dada pela equação (3.26). Para tanto, a partir de (b) avaliamos as seguintes derivadas: 02 2 u F ; T u F 2 2 ; 0 2 uu F Introduzindo as equações acima em (3.26), obtemos: 22 0 1 2 2 L I T dx Observamos que a segunda variação do funcional é sempre positiva, independentemente da função x , o que indica que o ponto estacionário determinado corresponde a um ponto de mínimo do funcional (a). D.A. Rade Fundamentos de Cálculo Variacional 50 3.2 Funcionais envolvendo derivadas de ordem superior. O procedimento desenvolvido na seção precedente pode ser estendido para a determinação dos pontos extremos de um funcionais envolvendo derivadas de ordem superior à primeira. Consideremos o funcional expresso dado por: dxxu,,xu,xu,xFuI x x n 2 1 (3.27) onde a função u x satisfaz às seguintes condições de contorno: 11 uxu ; 22 uxu 11 uxu ; 22 uxu (3.28) 1111 nn uxu ; 1221 nn uxu Por procedimento similar àquele apresentado na seção anterior, introduzimos a seguinte função de comparação: xuxuxu~ 21 xxx com: xxu onde 1 é um parâmetro escalar arbitrário e x é uma função contínua arbitrária definida no intervalo 21 xxx , que satisfaz às seguintes condições: 01 x ; 02 x 01 x ; 02 x (3.29) 011 xn ; 021 xn Em decorrência de (3.4), as funções de comparação satisfazem às condições de contorno do problema dadas por (3.28). A variação da função F resulta: nnn uFuFuFu,,u,u,xF (3.30) D.A. Rade Fundamentos de Cálculo Variacional 51 e a condição de estacionariedade fica: 02 1 dxu,u,u,xFI xx n (3.31) Introduzindo (3.30) em (3.31), obtemos: 2 1 0 x x )n( )n( dxu F u F u F (3.32) Desenvolvemos (3.32) fazendo integrações por partes sucessivas até a eliminação, no integrando, das derivadas n,, . Para tanto, utilizamos o seguinte desenvolvimento para a integração por partes: 2 2 1 1 1 x x k kk kk k x x F d Fdx dx dxu u 2 1 2 1 0 1 2 11 2 3 2 11 1 1 1 xk kk k k kk k k k x F d F d F d F dx dx dxu u u u ou, de forma mais compacta: 22 2 1 1 1 1 1 0 1 1 xx x k jk k jk k j k jk k k jx x x F d F d Fdx dx dx dxu u u (3.33) Introduzindo (3.33) em (3.32), a condição de estacionariedade do funcional toma a forma: 2 1 0 1 x kn k k k kx d F dx dx u 01 0 1 0 1 2 1 n k k j x x jk kj j j u F dx d (3.34) Em (3.34), a segunda parcela do lado esquerdo se anula em virtude das condições (3.29), o que conduz a: 2 1 0 1 0 x kn k k k kx d F dx dx u (3.35) D.A. Rade Fundamentos de Cálculo Variacional 52 Utilizando mais uma vez o Lema fundamental do Cálculo Variacional apresentado na seção precedente, concluímos que a equação (3.35) se verifica para qualquer função arbitrária x à condição que: 0 1 0 kn k k k k d F dx u (3.36) A equação (3.36) é a Equação de Euler-Lagrange associada ao funcional dado por (3.27). Os termos de contorno que figuram em (3.34): 01 0 1 0 1 2 1 n k k j x x jk kj j j u F dx d (3.37) representam o conjunto de todas as possíveis condições de contorno do problema. Exemplo 2: O problema de equilíbrio de uma viga de Euler-Bernoulli Consideremos a viga de uniforme ilustrada na Figura 3.3, de módulo de rigidez à flexão EI, submetida a um carregamento transversal distribuído q(x). Neste caso, com base na teoria de vigas de Euler-Bernoulli, a energia potencial total do sistema, incluindo a energia de deformação elástica e o trabalho da carga distribuída, associada a um campo de deslocamentos transversais representado por u(x), é dada pelo seguinte funcional: dxxyxq dx udEII L 0 2 2 2 2 1 (a) Figura 3.3 L q(x) x u u(x) D.A. Rade Fundamentos de Cálculo Variacional 53 Confrontando as equações (a) e (3.27), identificamos n=2 e 21 2 F u,u EI u q x u x (b) Particularizando (3.36) e (3.37) para n=2, obtemos: equação de Euler-Lagrange: 02 2 u F dx d u F dx d u F (c) termos de contorno: 0 0 L u F u F dx d u F (d) Do desenvolvimento de (c) a partir de (b), resulta a seguinte equação de Euler-Lagrange: EI xqu 4 , (e) que é a bem conhecida equação diferencial da linha elástica de vigas de Euler- Bernoulli. Desenvolvendo os termos de contorno dados por (d), obtemos a expressão: 00 LuEIuEI (f) Em (f), associamos a função ao deslocamento transversal da viga, à rotação da seção transversal. Relembramos ainda as seguintes definições: • uEIM : momento fletor • uEIV : esforço cortante Observamos que qualquer que seja o tipo de vínculo existente nas extremidades da viga, em x=0 e x=L, a equação (f) será satisfeita. Com efeito: • Para uma extremidade livre: M = 0 ; V = 0. • Para uma extremidade engastada: = 0 ; = 0 D.A. Rade Fundamentos de Cálculo Variacional 54 • Para uma extremidade simplesmente apoiada: M = 0 ; = 0 3.3 Funcionais envolvendo várias variáveis independentes Em diversos tipos de problemas de Engenharia, duas ou mais variáveis independentes podem intervir. É o caso, por exemplo, de problemas de equilíbrio bidimensionais nos quais a variável de campo são funções de duas coordenadas espaciais. Consideraremos, no desenvolvimento que segue, funcionais do tipo: D yx dydxu,u,y,x,uFuI (3.38) onde x e y são as variáveis independentes, D é um domínio bidimensional definido no plano x-y (ver Figura 3.4), y,xuu é a variável dependente e yx uu e designam, respectivamente, as derivadas da função u, em relação a x e y. Admitiremos, ainda, que a função u deva assumir um conjunto de valores prescritos ao longo de uma curva fechada C que constitui a fronteira do domínio D: y,xuy,xu CC (3.39) Figura 3.4 x y D C D.A. Rade Fundamentos de Cálculo Variacional 55 Para o estabelecimento das condições de estacionariedade do funcional (3.38), procedemos de forma análoga àquela adotada na seção anterior, introduzindo a função de comparação: y,xuy,xuy,xu , com: y,xy,xu (3.40) A função de comparação deve satisfazer as seguintes condições de fronteira: y,xuy,xu~ CC (3.41) De (3.39) e (3.41) decorre que a função y,x deve anular-se identicamente na fronteira C: 0Cy,x (3.42) A variação da função yx u,u,y,x,uF associada a y,xu é dada por: y y x x y y x x u F u F u Fu u Fu u Fu u FF (3.43) de modo que a condição necessária para que o funcional I tenha um ponto estacionário é: D dydxFI 0 D y y x x u F u F u F (3.44) o que conduz a: 0 D y y x x u F u F u F (3.45) Para efetuar as integrações por partes das duas últimas parcelas do integrando em (3.45), empregaremos o Teorema de Green em duas dimensões, expresso sob a forma: dxy,xFdyy,xGy,xdydx v y,xy,xF x y,xv,xG CD D.A. Rade Fundamentos de Cálculo Variacional 56 dydx y y,xF x y,xGy,x D (3.46) Utilizando (3.46), procedemos à integração dos dois últimos termos de (3.45) da seguinte forma: D y y x x dydx u F u F = dyuFdxuF yxC dydxu F dy d u F dx d D yx (3.47) Introduzindo (3.47) em (3.45), obtemos: dyuFdxuF yxC dydxu F dy d u F dx d u F D yx = 0 (3.48) Em virtude da condição (3.42) e do Lema 1, de (3.48) resulta: 0 yx u F dy d u F dx d u F (3.49) dyuFdxuF yxC = 0 (3.50) A equação (3.49) é a Equação de Euler-Lagrange associada ao funcional dado por (3.38) e os termos de contorno (3.50) fornecem todas as combinações possíveis de condições de contorno naturais e geométricas a serem satisfeitas na fronteira do domínio de definição do problema. D.A. Rade Fundamentos de Cálculo Variacional 57 Exemplo 3: O problema de equilíbrio de uma membrana tensionada Consideremos a membrana mostrada na Figura 3.5, a qual se encontra em equilíbrio estático, submetida a um carregamento transversal externo q(x,y) e a uma tensão constante T em seu plano. Figura 3.5 A energia potencial total, incluindo o trabalho da carga distribuída q(x,y) é dada pelo funcional: D yx dxdyquuuTuI 2221 (a) onde u=u(x,y) designa o campo de deslocamentos transversais da membrana. Confrontando a equação (a) com (3.38), identificamos: yx u,u,u,y,xF quuuT yx 2221 (b) Aplicando a equação (3.49), obtemos a seguinte equação de Euler- Lagrange: 022 quuT yx (c) Em particular, se D for um domínio circular de raio R, convém expressar (c) em coordenadas polares, introduzindo as transformações (ver Figura 3.6): y x z q(x,y) T T D.A. Rade Fundamentos de Cálculo Variacional 58 cosrx senrx 222 yxr Figura 3.6 Assim procedendo, (c) assume a forma: 0 q dr dur dr d r T (d) Admitindo que a carga distribuída seja constante, a integração da equação (d), subsidiada pela condição de contorno: 0,ru para r = R, conduz à seguinte expressão para o campo de deslocamentos transversais da membrana: 22 4 rR T qru Observemos que, devido à axissimetria geométrica e de carregamento, o deslocamento transversal é independente do ângulo . x y r P(x,y) R D.A. Rade Fundamentos de Cálculo Variacional 59 3.4 Funcionais envolvendo várias variáveis dependentes Nesta seção, consideraremos funcionais envolvendo mais de uma variável dependente e uma única variável independente. No que segue, trataremos de funcionais envolvendo apenas duas variáveis dependentes, sendo imediata a extensão dos desenvolvimentos ao caso de funcionais envolvendo um número maior de variáveis dependentes. Buscaremos estabelecer as condições de estacionariedade do seguinte funcional: dxv,u,v,u,xFv,uI x x xx 21 (3.51) onde x é a variável independente, u=(x) e v=(x) são as variáveis dependentes e xx vu e são, respectivamente, as derivadas de u e v em relação a x. Admitiremos que as funções u e v satisfaçam as seguintes condições de contorno: 11 uxu ; 22 uxu (3.52) 11 vxv ; 22 vxv Adotando procedimento análogo àquele empregado nas seções anteriores, introduzimos as seguintes funções de comparação: xuxuxu~ 21 xxx (3.53.a) xvxvxv~ 21 xxx (3.53.b) com: xxu (3.53.a) xxv (3.53.b) onde 1 é um parâmetro escalar arbitrário e x e x são funções contínuas arbitrárias definidas no intervalo 21 xxx , que devem assumir valores nulos nos extremos deste intervalo: 01 x 02 x (3.55.a) 01 x 02 x (3.55.b) D.A. Rade Fundamentos de Cálculo Variacional 60 Em decorrência de (3.55), as funções de comparação satisfazem as condições de contorno do problema: 11 uxu~ ; 22 uxu~ (3.56.a) 11 uxv~ ; 22 uxv~ (3.56.b) Empregando mais uma vez o procedimento baseado na expansão em série de Taylor limitada aos termos de primeira ordem, obtemos a seguinte expressão para a variação F : x x x x v F u F v F u FF (3.57) e a variação do funcional I correspondente à variação da função F resulta: 2 1 x x x x x x dx v F u F v F u FI (3.58) De (3.58) resulta a seguinte condição de estacionariedade do funcional I: 00 2 1 x x x x x x dx v F u F v F u FI (3.59) Procedemos, em seguida, à integração por partes das duas últimas parcelas do integrando, conforme detalhado abaixo: 2 1 2 1 2 1 x x x x xx x x x x dx u F dx d u Fdx u F (3.60.a) 2 1 2 1 2 1 x x x x xx x x x x dx v F dx d v Fdx v F (3.60.b) Introduzindo (3.60) em (3.50), obtemos: 2 2 2 1 1 1 0 x x x x x x xx x x F F F d F F d Fdx dx u v u dx u v dx v (3.61) Em virtude de (3.55), os termos de contorno em (3.61) se anulam: D.A. Rade Fundamentos de Cálculo Variacional 61 0 2 1 x xxx v F u F (3.62) de modo que, de (3.61), resulta a condição: 2 2 1 1 0 x x x xx x F d F F d Fdx dx u dx u v dx v Mais uma vez, evocando o Lema Fundamental do Cálculo Variacional empregado anteriormente, concluímos que a condição acima se verifica para quaisquer funções arbitrárias x e x se, e somente se: 0 x F d F u dx u (3.63.a) 0 x F d F v dx v (3.63.b) Os termos de contorno expressos por (3.62) representam todas as possíveis condições de contorno geométricas e naturais do problema, enquanto (3.63) constituem as Equações de Euler-Lagrange associadas ao funcional I. O procedimento desenvolvido acima pode ser imediatamente estendido ao caso de funcionais envolvendo um número qualquer de variáveis dependentes, representados genericamente sob a forma: dxz,,v,u,z,v,u,xFz,,v,uI x x xxx 21 (3.64) Neste caso, os termos de contorno resultam: 0 2 1 x xxxx z F v F u F (3.65) e haverá tantas equações de Euler-Lagrange quantas forem as variáveis dependentes: 0 x F d F u dx u (3.66.a) 0 x F d F v dx v (3.66.b) D.A. Rade Fundamentos de Cálculo Variacional 62 0 x F d F z dx z (3.66.c) Exemplo 4: O problema de vibrações de um sistema mecânico discreto de dois graus de liberdade. Para o sistema mecânico ilustrado na Figura 3.7, dado um conjunto de deslocamentos 1u e 2u e as respectivas velocidades 1u e 2u , são as seguintes as expressões para a energia cinética e a energia de deformação do sistema: 22221121 umumT (a) 212122220211021 uukukukV (b) Figura 3.6 Formemos o seguinte funcional: 2 1 t t dtVTI dtuukukukumumt t 2 1 2 1221 2 220 2 110 2 22 2 112 1 (c) Em (c), identificamos: 21221222021102222112121 21 uukukukumumu,u,u,uF (d) Empregando as equações (3.66), obtemos as seguintes equações de Euler- Lagrange: 10k 20k 12k 1m 2m tu1 tu2 D.A. Rade Fundamentos de Cálculo Variacional 63 0121211011 uukukum (e.1) 0121222022 uukukum (e.2) Podemos colocar as equações (e) sob a seguinte forma matricial: 10 12 121 1 1 12 20 122 2 2 0 0 0 0 k k km u u k k km u u (f.1) ou ainda: 0 uKuM (f.2) onde: 2 1 u u u (g.1) 2 1 0 0 m m M (matriz de inércia) (g.2) 10 12 12 12 20 12 k k k K k k k (matriz de rigidez) (g.3) As equações (f) constituem as equações do movimento em regime livre (sem forças externas) do sistema mecânico. Sua resolução, subsidiada por um conjunto de condições iniciais, conduzem à determinação das funções tu1 e tu2 , que descrevem o movimento das massas 1m e 2m em função do tempo. D.A. Rade Princípios Variacionais da Mecânica 64 CAPÍTULO 4 PRINCÍPIOS VARIACIONAIS DA MECÂNICA Neste Capítulo apresentaremos dois dos principais princípios variacionais da Mecânica: o Princípio da Energia Potencial Estacionária, aplicável a corpos deformáveis em equilíbrio estático e o Princípio de Hamilton, que rege os problemas de dinâmica de sistemas rígidos ou flexíveis. A partir do princípio de Hamilton, aplicando as técnicas do Cálculo Variacional, apresentadas no Capítulo 3, obteremos as Equações de Lagrange do Movimento. Será também apresentado o Método de Rayleigh-Ritz destinado à obtenção de soluções aproximadas de problemas regidos por princípios variacionais. Nos capítulos subseqüentes, o método de Rayleigh-Ritz, combinado com a noção de aproximação por subdomínios, será utilizado para a formulação do método dos elementos finitos, considerando alguns tipos particulares de elementos. 4.1 Fundamentos da Teoria da Elasticidade em pequenos deslocamentos Nesta seção recapitularemos os fundamentos da Teoria da Elasticidade com base na hipótese de pequenos deslocamentos. Esta hipótese admite que as componentes do deslocamento de um ponto qualquer P do corpo, denotadas doravante por u, v e w (ver Figura 4.1), sejam pequenas o bastante para justificar a linearização das equações que governam o problema. Figura 4.1 y x z v u w P’ P D.A. Rade Princípios Variacionais da Mecânica 65 A Figura 4.2 ilustra o estado mais geral de tensões a que pode estar submetido um elemento diferencial de volume de um corpo elástico, onde zyx ,, são as componentes de tensões normais, zyyzzxxzyxxy ,,,,, são as componentes de tensões cisalhantes e X , Y e Z são as componentes da força externa de volume (força por unidade de volume) atuante sobre o elemento. Figura 4.2 • Tensor das tensões. O estado de tensões em um ponto qualquer do corpo deformável é determinado por nove componentes formando o chamado tensor das tensões, definido da seguinte forma: zzyzx yzyyx xzxyx (4.1) y x z x xz xy z zx zy y yx yz Z Y X dx dy dz D.A. Rade Princípios Variacionais da Mecânica 66 • Equações de equilíbrio. Impondo as o equilíbrio de forças e momentos para o elemento diferencial de volume mostrado na Figura 4.1, pode-se mostrar que as componentes de tensão devem satisfazer as seguintes equações de equilíbrio: yxxy (4.2a) zxxz (4.2b) zyyz (4.2c) 0 X zyx xzxyx (4.2d) 0 Y zyx yzyxy (4.2e) 0 Z zyx zyzxz (4.2f) As equações (4.2a)-(4.2c) implicam a simetria do tensor das tensões que, doravante, será expresso em termos de apenas seis componentes independentes de tensões: zyzxz yzyxy xzxyx (4.3) • Tensor das deformações. O estado de deformações em um ponto qualquer do corpo deformável é determinado por nove componentes formando o chamado tensor das deformações, definido da seguinte forma: zzyzx yzyyx xzxyx (4.4) onde zyx ,, são as deformações normais e yzxzxy ,, são as deformações cisalhantes. D.A. Rade Princípios Variacionais da Mecânica 67 • Relações deformações-deslocamentos. Pode-se demonstrar que, admitindo a hipótese de pequenos deslocamentos, são válidas as seguintes relações envolvendo as componentes de deslocamento e as deformações: x u x (4.5a) y v y (4.5b) z w z (4.5c) x v y u xy (4.5d) x w z u xz (4.5e) y w z v yz (4.5f) • Relações tensões-deformações. Para materiais obedecendo à Lei de Hooke (proporcionalidade entre tensões e deformações), as chamadas equações constitutivas, relacionando as componentes de tensões e as componentes de deformações são expressas através da seguinte relação: E (4.6) onde: Tyzxzxyzyx (4.7) é o vetor contendo as componentes de tensões Tyzxzxyzyx (4.8) é o vetor contendo as componentes de deformações e D.A. Rade Princípios Variacionais da Mecânica 68 666564636261 565554535251 464544434241 363534333231 262524222221 161514131211 eeeeee eeeeee eeeeee eeeeee eeeeee eeeeee E (4.9) é o chamado tensor de rigidez ou tensor das constantes elásticas. As constantes elásticas satisfazem às seguintes relações, que conferem simetria ao tensor de rigidez. jiij ee para i,j = 1,2, ..., 6 (4.10) A relação (4.6) pode ser invertida, fornecendo: S (4.11) onde: SE 1 666564636261 565554535251 464544434241 363534333231 262524232221 161514131211 ssssss ssssss ssssss ssssss ssssss ssssss (4.12) é o chamado tensor de flexibilidade. Este é um tensor simétrico cujas componentes satisfazem as relações: jiij ss para i,j = 1,2, ..., 6 (4.13) Para materiais isotrópicos, cujas características mecânicas independem da direção, o número de constantes elásticas independentes fica reduzido a apenas dois. Neste caso, as relações tensões-deformações assumem a forma: zyxxx G 212 (4.14a) zyxyy G 212 (4.14b) zyxzz G 212 (4.14c) D.A. Rade Princípios Variacionais da Mecânica 69 xyxy G (4.14d) xzxz G (4.14e) yzyz G (4.14f) e, inversamente: zyxx E 1 (4.15a) zxyy E 1 (4.15b) yxzz E 1 (4.15c) com: 12 EG (4.16) • Condições de contorno. No que diz respeito às condições de contorno, podemos imaginar a superfície externa do corpo deformável dividida em duas regiões: a região 1S , sobre a qual as condições de contorno são estabelecidas em termos de forças externas aplicadas e a região 2S , sobre a qual as condições de contorno são expressas em termos de deslocamentos impostos. Denotando por Z,Y,X as componentes cartesianas as forças impostas por unidade de área em 1S , as condições de contorno que traduzem o equilíbrio de forças na fronteira são expressas segundo: XX YY ZZ em 2S (4.17) onde: nmX xzxyx (4.18a) nmY yzyxy (4.18b) nm zyzxz (4.18c) D.A. Rade Princípios Variacionais da Mecânica 70 As equações (4.18) são obtidas impondo o equilíbrio do tetraedro infinitesimal mostrado na Figura 4.3, onde knjmin é o vetor unitário na direção da normal exterior à superfície 1S . Figura 4.3 Denotando as componentes de deslocamentos impostos por w,v,u , as condições de contorno ditas geométricas, aplicáveis à região 2S , são expressas segundo: uu (4.19a) vv em 2S (4.19b) ww (4.19c) No problema geral de Elasticidade temos um total de 15 incógnitas: 3 componentes de deslocamento, 6 componentes de deformação e 6 componentes de tensão. Estas incógnitas podem ser determinadas mediante a resolução de 15 equações diferenciais: 3 equações de equilíbrio (4.2), 6 relações deformações- deslocamentos (4.5) e 6 relações tensões-deformações (4.6). A solução deve também satisfazer as condições de contorno (4.17) e (4.19). 4.2 O Princípio do Trabalho Virtual em Elasticidade Considerando um corpo deformável em equilíbrio sob a ação de forças externas e um conjunto de restrições cinemáticas impostas, introduzimos o conceito de deslocamentos virtuais: kwjviur (4.20) x y z n z zx zy y yx yz x xy xz dS1 D.A. Rade Princípios Variacionais da Mecânica 71 entendidos como um campo de deslocamentos infinitesimais arbitrários aplicados aos pontos do corpo. Os deslocamentos virtuais são deslocamentos imaginários, diferindo dos deslocamentos reais pelo fato que a eles não associamos nenhuma variação do tempo. Deste modo, as forças atuantes sobre o corpo permanecem inalteradas durante a ocorrência dos deslocamentos virtuais. Admite-se ainda que os deslocamentos virtuais não violem as condições de contorno geométricas. Desta forma, devemos ter: 0u 0v 0w em 2S (4.21) Define-se o trabalho virtual de uma força F como sendo o trabalho realizado por esta força quando seu ponto de aplicação sobre um deslocamento virtual: rFW (4.22) Em termos de componentes cartesianas: wFvFuFkwjviukFjFiFW zyxzyx (4.23) Com base nestas definições, é a seguinte a expressão para o trabalho virtual das forças externas impostas sobre o contorno do corpo: dSwZvYuXW S S 1 (4.24) Considerando as condições de contorno na fronteira 1S , dadas por (4.17), a equação (4.24) se reescreve: dSwnmvnmunmW S zyzxzyzyxyxzxyx S 1 (4.25) Aplicando o Teorema de Gauss, expresso sob a forma: dV z N y M x LdSNnMmL S V (4.26) a equação (4.25) fica: dVwvu z wvu y wvu x W V yzyzxzyzyxyxzxyx S Desenvolvendo as derivadas indicadas na equação acima, temos: D.A. Rade Princípios Variacionais da Mecânica 72 dVw zyx v zyx u zyx W V zxyxzyzyxyxzxyxS + dV z w z v y w y v z u x w y u x v x u V zyzyxzxyx (4.27) Introduzimos agora as seguintes deformações associadas aos deslocamentos virtuais: z w y v x u zyx (4.28) z v y w z u x w y u x v yzxzxy Utilizando as definições (4.28) e levando em conta as equações de equilíbrio (5.2), reescrevemos (4.27) sob a forma: dVwZvYuXW V S + dV V yzyzxzxzxyxyzzyyxx (4.29) Por fim, confrontando (4.24) e (4.29), escrevemos: dSwZvYuX S 1 dVwZvYuX V = dV V yzyzxzxzxyxyzzyyxx (4.30) ou ainda: IVS WWW (4.31) onde: D.A. Rade Princípios Variacionais da Mecânica 73 • dSwZvYuXW S S 1 (4.32) é o trabalho virtual das forças externas impostas na fronteira S1; • dVwZvYuXW V V (4.33) é o trabalho virtual das forças externas de volume; • dVW V yzyzxzxzxyxyzzyyxx I (4.34) é o trabalho virtual das forças internas, associadas às tensões distribuídas sobre o volume do corpo. As equações (4.30) e (4.31) expressam o chamado Princípio do Trabalho Virtual para corpos deformáveis em equilíbrio, sob a hipótese de pequenos deslocamentos. 4.3 O Princípio da Energia Potencial Estacionária Consideremos os casos em que seja possível definir uma função energia potencial de deformação ,,,,,,UU yzxzxyzyx de modo que: x x U y y U z z U (4.35) xy xy U xz xz U yz yz U No caso de materiais obedecendo as relações constitutivas (4.6), podemos facilmente verificar que a energia potencial de deformação é dada por: EU TT 2 1 2 1 (4.36) Mais particularmente no caso particular de materiais isotrópicos, regidos pelas equações constitutivas (4.14), o desenvolvimento de (4.36) conduz à expressão: 2222222 212121 yzxzxyzyxzyxU (4.37) Com base em considerações físicas, pode-se mostrar que, em (4.36) U é uma função definida positiva das componentes de deformação. Isto significa que U>0 para qualquer campo de deformações correspondente a um campo de D.A. Rade Princípios Variacionais da Mecânica 74 deslocamentos cinematicamente admissível. Tal fato fica evidenciado na forma particular (4.37). Introduzindo (4.35) em (4.34), podemos expressar o trabalho virtual das forças internas em termos da energia potencial de deformação da seguinte forma: dVUUUUUUW V yz yz xz xz xy xy z z y y x x I Na última equação acima observamos que o integrando representa o diferencial total da função ,,,,,,UU yzxzxyzyx . Assim, podemos escrever: dVUW V I (4.38) com: yz yz xz xz xy xy z z y y x x UUUUUUU (4.39) Podemos ainda permutar a ordem dos operadores integração e diferenciação em (4.38), escrevendo: dVUW V I (4.40) Consideremos agora o trabalho virtual das forças externas (forças de volume e forças impostas na fronteira S1. Admitindo que estas forças sejam conservativas (o trabalho por elas realizado independe do caminho percorrido por seu ponto de aplicação), suas componentes podem ser expressas em termos de derivadas parciais de funções potenciais da seguinte forma: u X v Y w Z (5.41) u X vY wZ onde w,v,u e w,v,u são, respectivamente, as funções potenciais associadas às forças externas de volume e às forças externas impostas na fronteira S1. Nestas condições, partindo de (4.32) e (4.33), o trabalho virtual das forças externas pode ser expresso sob a forma: D.A. Rade Princípios Variacionais da Mecânica 75 VS IS dVw w v v u u dSw w v v u u WW 1 Os integrandos na equação acima correspondem aos diferenciais totais das funções e , fato que nos permite escrever: dVwvudSwvuWW SS IS 11 ,,,, (4.42) Assim, combinando as equações (4.38) e (4.42), o Princípio do Trabalho Virtual, expresso por (5.31), assume a forma: 0 1 V V S dSdVdVU (4.43a) ou ainda: 0 (4.43b) onde: V V S dSdVdVU 1 (4.44) é a energia potencial total. As equações (4.43) expressam o chamado Princípio da Energia Potencial Estacionária, o qual estabelece que, dentre todos os campos de deslocamentos (u,v,w) que satisfazem as condições de contorno geométricas, aquele que corresponde à real configuração de equilíbrio é um ponto estacionário da energia potencial total. Do ponto de vista do Cálculo Variacional, podemos observar que as equações de equilíbrio e as condições de contorno naturais em S1 são, respectivamente, as equações de Euler-Lagrange e os termos de contorno do problema de estacionariedade do funcional (4.44). Resta mostrar que o extremo associado ao campo de deslocamentos da configuração de equilíbrio estático corresponde a um ponto de mínimo do funcional (4.44). Para tanto, denotaremos, a seguir, por (u,v,w) as componentes de deslocamento correspondentes à configuração de equilíbrio e por wvu ,, um campo de deslocamentos arbitrário, próximo a (u,v,w), que satisfazem as condições de contorno geométricas. Relacionaremos ambos os campos de deslocamentos através das relações: uuu vvv www (4.45) D.A. Rade Princípios Variacionais da Mecânica 76 Introduzindo (4.45) em (4.44) e empregando uma expansão em série de Taylor limitada à segunda ordem, escrevemos: wvuwvuwvuwvu ,,,,,,,, 2 Como (u,v,w) é um ponto estacionário do funcional , temos 0,, wvu e a equação acima fica: wvuwvuwvu ,,,,,, 2 (4.46) A equação (5.46) nos mostra que a natureza do ponto crítico (mínimo ou máximo) é determinada pelo sinal da segunda variação wvu ,,2 , de modo que: • 0,,2 wvu ponto de mínimo • 0,,2 wvu ponto de máximo Podemos mostrar que: dVwvuU V ,,2 (4.47) onde U é a energia potencial de deformação, dada por (4.36). Do fato de U ser uma função definida positiva resulta que 02 para quaisquer valores arbitrárias de wvu e, , o que demonstra que o campo de deslocamentos correspondente à posição de equilíbrio do corpo deformável corresponde a um ponto de mínimo do funcional que representa a energia potencial total. Por esta razão, o Princípio da Energia Potencial Estacionária é às vezes denominado Princípio da Energia Potencial Mínima. 4.4 – Princípio do Trabalho Virtual Aplicado à Dinâmica de Sistemas de Partículas Visando deduzir o Princípio Variacional de Hamilton, aplicável a sistemas mecânicos em movimento, estabeleceremos primeiramente o Princípio do Trabalho Virtual para um conjunto arbitrário de n partículas, ilustrado na Figura 4.4, ficando subentendido que qualquer sistema mecânico discreto ou contínuo poderá ser considerado como formas particulares deste tipo de conjunto. Na Figura 4.4, sobre uma partícula genérica iP , de massa im , estão indicadas: • iF : a resultante das forças impostas sobre a partícula, representando as ações mecânicas dos corpos a ela vizinhos, incluindo as forças exercidas pelas outras partículas do sistema; D.A. Rade Princípios Variacionais da Mecânica 77 • if : a força de restrição exercida sobre pela superfície de restrição iS sobre a qual iP é forçada a se mover. Negligenciaremos, por enquanto, o atrito, de modo que a força de restrição terá sempre a direção normal à superfície de restrição, indicada na Figura 7.1 por in . • ir : vetor posição da partícula em relação ao sistema de referência, suposto fixo, OXYZ . Figura 4.4 Aplicando a Segunda Lei de Newton para cada partícula do sistema, considerada isoladamente, escrevemos: iiii rmfF i=1, 2, ..., n (4.48) O Princípio de d’Alembert, apresentado na Seção 3.6, nos permite escrever (4.49) sob a forma: 0 iiii rmfF i=1, 2, ..., n (4.49) com a interpretação de que, para um observador movimentando-se juntamente com a partícula, esta é observada em estado de equilíbrio sob a ação das forças Y nP 1r X Z ir nr iP 1P O 1F iF nF nf if 1f 1r ir nr it in 1S iS nS D.A. Rade Princípios Variacionais da Mecânica 78 iF , if e da força de inércia ii rm . Esta situação é denominada equilíbrio dinâmico. Introduzimos agora o conceito dos chamados deslocamentos virtuais, que serão denotados por 1r , ..., ir , ..., nr . Tais deslocamentos virtuais são concebidos com as seguintes propriedades: • são perturbações imaginárias, infinitesimais e arbitrárias das posições das partículas, que não violam as restrições cinemáticas impostas pelas superfícies de restrição. Esta última condição implica que os deslocamentos virtuais devem ter a direção tangente à superfície de restrição, na posição instantaneamente ocupada pela partícula. Assim, conforme ilustra a Figura 7.1, os vetores if e ir são mutuamente perpendiculares. • são admitidos ocorrer instantânea e simultaneamente, de modo que aos deslocamentos virtuais não associamos nenhum lapso de tempo finito. Desta forma, as forças aplicadas às partículas do sistema não variam em decorrência dos deslocamentos virtuais. Em termos de coordenadas cartesianas, os vetores posição das partículas são dados por: kZjYiXr iiii , i=1, 2, ..., n de modo que os deslocamentos virtuais podem ser expressos segundo: kZjYiXr iiii , i=1, 2, ..., n onde iii ZYX ,, são entendidos como variações dadas às coordenadas que representam a posição das partículas em relação ao sistema de referência empregado. Em virtude de (4.49), podemos afirmar que, na situação de equilíbrio dinâmico, o trabalho realizado por todas as forças aplicadas sobre a partícula iP , incluindo a força de inércia, associado ao deslocamento virtual ir é nulo, ou seja: 0 iiiii rrmfF i=1, 2, ..., n (4.50) Do fato dos vetores if e ir serem mutuamente perpendiculares resulta que as forças de restrição produzem trabalho virtual nulo ( 0 ii rf ) e (4.50) fica: 0 iiii rrmF i=1, 2, ..., n (4.51) Adicionando as n equações (4.51), obtemos: D.A. Rade Princípios Variacionais da Mecânica 79 0 1 n i iiii rrmF (4.52) ou ainda: 0 IF WW (4.53) onde: n i ii F rFW 1 e IW n i iii rrm 1 (4.54) desingam o trabalho virtual das forças impostas e o trabalho virtual das forças de inércia, respectivamente. A equação (4.53) expressa o Princípio do Trabalho Virtual, o qual estabelece que para qualquer posição de um sistema de partículas, o trabalho virtual total realizado por todas as forças impostas e todas as forças de inércia resulta nulo para todo e qualquer conjunto arbitrário de deslocamentos virtuais introduzidos a partir daquela posição. 4.5 – Princípio de Hamilton No desenvolvimento que segue, ao operador , que designa uma variação virtual, serão atribuídas as mesmas propriedades do operador diferencial do Cálculo Diferencial. A fim de desenvolver (4.52), introduzimos a identidade: iiiiiiiiii rrrrrrrrdt rrd 2 1 donde: iiiiii rrdt rrdrr 2 1 (4.55) Introduzindo (4.55) em (4.52) e desenvolvendo a equação resultante, escrevemos: n i ii rF 1 0 2 1 11 n i iii n i ii i rrmdt rrdm (4.56) Em (4.56), reconhecemos: D.A. Rade Princípios Variacionais da Mecânica 80 n i iii rrmT 12 1 (energia cinética do sistema de partículas), (4.57) de modo que, considerando (4.54), podemos reescrever (4.56) sob a forma: TW F n i ii i dt rrd m 1 (4.57) Observamos que, no caso geral, a energia cinética pode ser expressa como uma função de várias variáveis. Em termos das componentes cartesianas dos vetores posição e velocidade das partículas do sistema, tal função assume a forma: nnnnnn zyxzyxzyxzyxzyxzyxTT ,,,,,,,,,,,,,,,,,,, 222111222111 , (4.58) de modo que T designa o diferencial total da função T, expresso segundo: n n n n n n z z Ty y Tx x Tz z Ty y Tx x TT 1 1 1 1 1 1 n n n n n n z z Ty y Tx x Tz z Ty y Tx x T 1 1 1 1 1 1 (4.59) Multiplicando (4.57) por dt e integrando entre dois instantes de tempo 1t e 2t , temos: dtTWtt F21 2 1 2 1 11 t t n i iii t t n i ii i rrmdtdt rrd m (4.60) Neste ponto, admitiremos que os deslocamentos virtuais, embora arbitrários, sejam nulos nos instantes 1t e 2t ., ou seja, 021 trtr ii , i=1,2,...,n, de modo que, da equação (4.60), resulta: 02 1 dtTWtt F (4.61) Admitindo que as forças impostas sejam todas conservativas, podemos escrever: VW F (4.62) onde V é a função potencial ou energia potencial associada às forças impostas, que depende exclusivamente das posições ocupadas pelas partículas do sistema. Em coordenadas cartesianas, esta função é expressa segundo: D.A. Rade Princípios Variacionais da Mecânica 81 nnn z,y,x,,z,y,x,z,y,xVV 222111 (4.63) de modo que V é interpretado como o diferencial total da função V, dado por: n n n n n n z z Vy y Vx x Vz z Vy y Vx x VV 1 1 1 1 1 1 (4.64) Assim, introduzindo (4.62) em (4.61), escrevemos: 2 1 0 t t dtVT ou ainda: 0 2 1 t t dtL (4.65) onde: VTL (4.66) é o chamado Lagrangeano do sistema. Levando em conta (4.59) e (4.64), observamos que o Lagrangeano assume a forma de uma função de várias variáveis, em termos das componentes dos vetores posição e velocidade da partículas. Assim, escrevemos: nnnnnn zyxzyxzyxzyxzyxzyxLL ,,,,,,,,,,,,,,,,,,, 222111222111 , (4.67) e, n n n n n n z z Ly y Lx x Lz z Ly y Lx x LL 1 1 1 1 1 1 n n n n n n z z Ly y Lx x Lz z Ly y Lx x L 1 1 1 1 1 1 (4.68) Podemos permutar os operadores integração e diferenciação, expressando então (4.65) sob a seguinte forma: 2 1 0 t t dtL (4.69) D.A. Rade Princípios Variacionais da Mecânica 82 A equação (4.69) expressam o Princípio Variacional de Hamilton, que se aplica aos casos em que todas as forças impostas são conservativas. Este princípio estabelece que dentre todos os movimentos possíveis de ser realizados pelo sistema mecânico, satisfazendo as restrições cinemáticas impostas, entre dois instantes de tempo subseqüentes quaisquer t1 e t2, o movimento realmente desenvolvido pelo sistema será aquele que corresponde a um ponto estacionário do funcional: 2 1 t t dtLI (4.70) Conforme ilustrado no exemplo a seguir, a imposição das condições necessárias para a ocorrência de um ponto crítico, estudadas no âmbito do Cálculo Variacional, abordado no Capítulo 3, conduz às equações do movimento do sistema mecânico, que são determinadas pelas equações de Euler-Lagrange associadas ao funcional expressso em (4.70). 4.6 O Princípio de Hamilton Estendido Quando, além das forças impostas conservativas, houver forças impostas não conservativas, podemos, em (4.61), separar os trabalhos realizados por estes dois tipos de forças, escrevendo: FncFcF WWW (4.71) onde FncFc WW e designam, os trabalhos realizados pelas forças impostas conservativas e não conservativas, respectivamente. Fazendo uso da relação (4.62), para as forças impostas conservativas, escrevemos: VW Fc (4.72) Associando (4.71) e (4.72) com (4.61), escrevemos: 02 1 2 1 dtWdtL tt Fnctt (4.73) A equação (7.27) traduz o chamado Princípio de Hamilton Estendido, aplicável aos casos em que existem forças impostas não conservativas. É importante observar que, no lado esquerdo de (4.73), o termo FncW , que representa o trabalho virtual das forças não conservativas não pode ser expresso sob a forma de diferencial total de uma função de várias variáveis. Em conseqüência (4.73) não pode ser interpretada como a condição de ocorrência de um ponto crítico de uma função de várias variáveis, o que significa que, estritamente, o Princípio de Hamilton Estendido não pode ser considerado como D.A. Rade Princípios Variacionais da Mecânica 83 um princípio variacional. Contudo, (4.73) pode ser utilizada para obter as equações do movimento, conforme mostraremos com o auxílio do exemplo apresentado a seguir. Uma outra observação importante a ser feita é que, embora tenham sido desenvolvidos para sistemas de partículas, o Princípio de Hamilton e o Princípio de Hamilton Estendido por ser imediatamente aplicáveis a sistemas mecânicos formados por corpos rígidos e deformáveis, desde que sejam corretamente computadas suas energias cinética e potencial e os trabalhos virtuais das forças e momentos aplicadas a estes corpos. Exemplo 4.1: Sabendo que no sistema ilustrado a mola está indeformada quando r=ro, obteremos as equações diferenciais do movimento. Na figura abaixo, t representa um torque externo, m e mo designam, respectivamente, a massa do cursor e da barra uniforme OA e k é a constante de rigidez da mola. Obtenhamos as equações do movimento do sistema em termos das grandezas r e definidas na figura. Considerando o torque externo t como um esforço não conservativo, utilizaremos o Princípio de Hamilton Estendido, dado pela equação (7.27), repetida abaixo: 02 1 2 1 dtWdtL tt Fnctt (a) mo m L O A r k r t Ref. D.A. Rade Princípios Variacionais da Mecânica 84 Em termos das coordenadas radial-transversal indicadas na figura, escrevemos as seguintes expressões para a energia cinética e potencial do sistema: barracursor TTT onde: 2222 2121 rrmvvmT rcursor 222 6 1 2 1 LmJT oobarra (rotação não baricêntrica) (elást.)(gravit.)(gravit.) molabarracursor VVVV onde: cosmgrVcursor cos 2 LgmV obarra 2 2 1 omola rrkV Com base nas equações acima, escrevemos o Lagrangeano sob a forma: 2221 rrmL 2261 Lmo cosmgr cos2Lgmo 221 orrk (b) O trabalho virtual do torque não conservativo é dado por: FncW (c) Com o entendimento que o operador funciona como o diferencial total, escrevemos: LLr r Lr r LL (d) D.A. Rade Princípios Variacionais da Mecânica 85 A partir de (b), avaliamos as derivadas parciais indicadas na equação acima: 2mr r L cosmgrrk o rmr L sen2sen LgmmgrL o 22 3 1 LmmrL o Introduzindo estas derivadas em (d), escrevemos: rrmrmgrrkmrL o cos2 222 3 1sen 2 sen LmmrLgmmgr oo (e) Introduzindo (c) e (e) em (a), obtemos: rrmrmgrrkmrt t o 2 1 cos2 0 3 1 2 222 dtLmmrsenLgmsenmgr oo (f) Desenvolvemos, em seguida, eliminar os diferenciais das variações temporais das coordenadas que aparecem em (f). Para tanto, efetuamos as seguintes integrações por partes: • 2 1 2 1 2 1 t t t t t t dtrrrrmdtrrm • 2 1 2 1 2222 3 1 3 1t t t t oo LmmrdtLmmr dtLmmrrmrt t o 12 22 3 12 Introduzindo as duas últimas equações em (f), após rearranjo, escrevemos: 2 1 2 t t o rrmcosmgrrkmr dtLmmrrmrsenLgmsenmgr oo 22 3122 D.A. Rade Princípios Variacionais da Mecânica 86 0 3 1 2 1 22 t t o Lmmrrrm (g) Relembrando que e r devem ser nulas em t=t1 e t=t2, o último termo do lado esquerdo da equação acima resulta nulo. Além disso, como a equação resultante deve ser satisfeita para quaisquer incrementos arbitrários e r , (g) é satisfeita se forem nulas as funções que multiplicam estes incrementos. Deste fato, resultam as equações diferenciais do movimento do sistema: 02 cosmgrrkrrm o (h) tgLmmrrmrLmmr oo sen 2 2 3 1 22 (i) Observemos que as equações do movimento são acopladas e não lineares. A partir de um conjunto de condições iniciais e dada a função que define o torque externo aplicado t , estas equações podem ser integradas numericamente (ver Seção 3.8) para obtenção das funções que tr e t que definem a configuração do sistema em função do tempo. 4.7 Número de graus de liberdade e coordenadas generalizadas Considerando o sistema de n partículas mostrado na Figura 4.4, podemos afirmar que, se todas as partículas puderem se mover livre e independentemente uma das outras, a posição de cada uma delas fica determinada por três coordenadas espaciais iii zyx ,, i=1,2,...,n, de modo que, para a completa determinação da configuração espacial do sistema são necessárias 3n coordenadas. Contudo, se houver restrições cinemáticas que limitam o movimento das partículas a uma dada região do espaço ou impedem algum tipo de movimento relativo entre as partículas, algumas destas 3n coordenadas tornam-se dependentes entre si. Neste caso, o número de coordenadas independentes necessárias e suficientes para determinar a configuração espacial do sistema é dado por: N = 3n – p (4.74) onde p é o número de equações de restrição impostas pelas restrições cinemáticas e N é o denominado número de graus de liberdade do sistema. No caso de sistemas formados por corpos rígidos, o número de coordenadas necessárias para definir a posição espacial de cada corpo é 6 (3 coordenadas lineares e 3 coordenadas angulares), de modo que o número de graus de liberdade de um sistema formado por n corpos rígidos é dado por: N = 6n – p (4.75) D.A. Rade Princípios Variacionais da Mecânica 87 A escolha do conjunto de coordenadas independentes é arbitrária, embora seu número permaneça invariável. Qualquer conjunto de N coordenadas independentes constitui um conjunto de coordenadas generalizadas, que serão aqui denotadas por Nqqq ,,, 21 . Ilustramos os conceitos de número de graus de liberdade e coordenadas generalizadas com o auxílio do seguinte exemplo. Exemplo 4.2: No sistema ilustrado na figura abaixo, as partículas P1 e P2 são forçadas a se mover sobre o plano x-y, estando ligadas pela barra rígida de comprimento L. Para este sistema, determinemos o número de graus de liberdade e um conjunto de coordenadas generalizadas. No caso em questão temos um sistema com duas partículas (n=2) e número toal de coordenadas é 3n=6. Em termos das componentes de seus vetores posição em relação ao sistema de eixos indicados na figura, estas seis coordenadas são denotas por 222111 z,y,x,z,y,x . No entanto, estas seis coordenadas não são todas indepentes, uma vez que as seguintes as equações de restrição devem ser satisfeitas. z1=0 z2=0 2212212 Lyyxx P1 P2 L y x z D.A. Rade Princípios Variacionais da Mecânica 88 Temos então p = 3 (três equações de restrição), de modo que, de acordo com a equação (4.74), o número de graus de liberdade do sistema resulta ser N=323=3. Alguns conjuntos possíveis de coordenadas generalizadas são: ( 11 xq , 12 yq , 23 xq ) ( 11 xq , 12 yq , 13 yq ) ( 11 xq , 12 yq , 3q ) 4.8 Equações de Lagrange Nesta seção, desenvolveremos a formulação que, partindo do Princípio de Hamilton Estendido, conduz às chamadas Equações de Lagrange, as quais constituem uma forma bastante elegante e expedita para a obtenção das equações do movimento de sistemas dinâmicos. As equações de Lagrange são formuladas em termos de coordenadas generalizadas, cujo conceito foi introduzido na Seção 4.7. Para o conjunto de n partículas mostrado na Figura 4.4, podemos sempre expressar os vetores posição das partículas em função de um conjunto previamente escolhido de N coordenadas generalizadas através de relações do tipo: tqqqrr Nii ,,,, 21 i=1,2, ..., n (4.76) Em termos de coordenadas cartesianas, a transformação (4.76) pode ser detalhada como segue: kzjyixr iiii (4.77) com: tqqqxx Nii ,,, 21 tqqqyy Nii ,,, 21 (4.78) tqqqzz Nii ,,, 21 i=1,2, ..., n Buscaremos, em seguida, representar as velocidades das partículas em termos das coordenadas generalizadas. Para tanto, empregando a regra da cadeia da derivação, derivamos (4.77) em relação ao tempo, escrevendo: kzjyixv iiii (4.79) D.A. Rade Princípios Variacionais da Mecânica 89 com: t x q q x t x q q x q q x q q x dt dx x i N j j j ii N N iiii i 1 2 2 1 1 t y q q y q q y q q y dt dy y iN N iiii i 2 2 1 1 t yq q y iN j j j i 1 (4.80) t z q q z q q z q q z dt dz z iN N iiii i 2 2 1 1 t z q q z iN j j j i 1 i=1,2, ..., n Levando em conta (4.79), a expressão da energia cinética do sistema de partículas, dada por (4.57), é desenvolvida da seguinte forma: n i iiii n i iii zyxmvvmT 1 222 1 2 1 2 1 (4.81) Introduzindo (4.80) em (4.81), escrevemos: 2 1 2 1 2 112 1 t zq q z t yq q y t xq q xmT i N j j j ii N j j j ii N j j j i n i i Na equação acima podemos perceber que, de modo geral, a energia cinética será função das coordenadas generalizadas nqqq ,, 21 , de suas derivadas temporais nqqq ,, 21 e do tempo t. Assim, escrevemos: tqqqqqqTT NN ,,,,,, 2121 (4.82) Da mesma forma, a energia potencial, que é função exclusiva dos vetores posição (e, eventualmente, do tempo), é expressa sob a forma: tqqqqqqVV NN ,,,,,, 2121 (4.83) Com base em (4.82) e (4.83), podemos expressar o Lagrangeano como uma função das coordenadas generalizadas, das coordenadas generalizadas nqqq ,, 21 , de suas derivadas temporais nqqq ,, 21 e do tempo t: tqqqqqqLVTL NN ,,,,,, 2121 (4.84) A partir de (4.84), podemos expressar da seguinte forma a variação do Lagrangeano associada a um conjunto arbitrário de deslocamentos virtuais: D.A. Rade Princípios Variacionais da Mecânica 90 N N N N q q Lq q Lq q Lq q Lq q Lq q LL 2 2 1 1 2 2 1 1 ou: N j j j j j q q Lq q LL 1 (4.85) Notemos que, de acordo com as propriedades atribuídas aos deslocamentos virtuais (ver Seção 4.4), a variação no tempo, representada pelo termo dttL não é incluída no desenvolvimento que conduz a (4.85). No desenvolvimento que segue, buscaremos expressar o trabalho virtual das forças não conservativas em termos das coordenadas generalizadas. Para tanto, escrevemos: n i ii F nc rFW 1 (4.86) Fazendo uso da equação (4.77), expressamos os deslocamentos virtuais sob a forma: ir N N iii q q r q q r q q r 2 2 1 1 N j j j i q q r 1 (4.87) onde os termos jq , j=1,2, ..., N, são entendidos como variações arbitrárias e independentes introduzidas nas coordenadas generalizadas. Notemos, mais uma vez que, de acordo com as propriedades atribuídas aos deslocamentos virtuais na Seção 4.4, a variação no tempo, representada pelo termo dttri não é incluída no desenvolvimento (4.87). Introduzindo (4.87) em (4.86), escrevemos: FncW N j j n i j i i qq rF 1 1 (4.88) ou ainda: FncW N j jj qQ 1 (4.89) onde: D.A. Rade Princípios Variacionais da Mecânica 91 n i j i ij q rFQ 1 (4.90) são os denominados esforços generalizados. Voltemos agora ao Princípio de Hamilton Estendido, expresso pela equação (4.73), repetida abaixo: 02 1 2 1 dtWdtL tt Fnctt (4.91) Introduzindo (7.85) e (4.89) em (4.91), escrevemos: dtq q Lq q Lt t N j j j j j 2 1 1 0 2 1 1 dtqQ t t N j jj Rearrajando: 0 2 1 1 1 dtq q LqQ q Lt t N j N j j j jj j (4.92) Efetuamos em seguida as seguintes integrações por partes: dtq q L dt dq q Ldtq q L t t N j j j t t j j t t N j j j 2 1 2 1 2 1 11 (4.93) Considerando a hipótese que os deslocamentos virtuais (e, portanto, as variações correspondentes das coordenadas generalizadas) devem anular-se nos instantes 1t e 2t , ou seja, 021 tqtq jj , o primeiro termo do lado direito de (4.93) resulta nulo. Introduzindo então a equação resultante em (4.92), escrevemos, após alguns rearranjos: 0 2 1 1 dtqQ q L q L dt dt t N j jj jj (4.94) Uma vez que a equação acima deve ser satisfeita para todo e qualquer conjunto de variações virtuais arbitrárias e independentes, jq , j=1,2, ..., N, resulta que as funções que multiplicam cada uma das variações jq em (4.94) devem anular-se identicamente. Assim, temos: j jj Q q L q L dt d j=1,2, ..., N (4.95) D.A. Rade Princípios Variacionais da Mecânica 92 As equações (4.95) são as denominadas Equações de Lagrange do movimento. Note-se que tais equações apresentam-se sob a forma de equações diferenciais de segunda ordem na variável tempo, que constituem as equações do movimento do sistema mecânico. Devemos também observar que, embora tenham sido deduzidas para sistemas de partículas, as Equações de Lagrange na forma (4.95) podem ser imediatamente aplicadas a sistemas constituídos de corpos rígidos e deformáveis, bastando que se usem as expressões adequadas para o cálculo das energias cinética e potencial que compõem o Lagrangeano. O uso das equações de Lagrange para a obtenção das equações do movimento é ilustrado no exemplo a seguir. Exemplo 4.3: Considerando o sistema mecânico apresentado do Exemplo 4.1, obtenhamos as equações diferenciais do movimento utilizando as Equações de Lagrange. Como as coordenadas ,r são independentes, podemos adotá-las como coordenadas generalizadas e as duas equações de Lagrange do movimento, expressas por (4.95) tomam as formas: rQr L r L dt d (i) Q LL dt d (ii) Sendo o trabalho virtual do torque não conservativo, dado por: FncW e desenvolvendo (4.89) considerando a escolha feita para as coordenadas generalizadas: FncW QrQr , constatamos que: 0rQ Q (iii) Efetuando as derivadas parciais indicadas em (i) e (ii), a partir da expressão do Lagrangeano dada pela Eq. (a) do Exemplo 4.1, obtemos: 2221 rrmL 2261 Lmo cosmgr cos2Lgmo 221 orrk (b) D.A. Rade Princípios Variacionais da Mecânica 93 r L rm rmr L dt d r L orrkmgmr cos2 L 22 3 1 Lmmr o L dt d 22 3 12 Lmmrrmr o L senLgmmgrsen o 2 (iv) Introduzindo as equações (iii) e (iv) em (i) e (ii), obtemos, as seguintes equações do movimento, que são idênticas às equações (h) e (i) obtidas no Exemplo 4.1: 0cos2 mgrrkrrm o (h) tgsenLmmrrmrLmmr oo 2 2 3 1 22 (i) 4.9 Método de Rayleigh-Ritz Conforme vimos nas seções anteriores deste capítulo, os métodos variacionais conduzem às equações diferenciais de equilíbrio, no caso de problemas estáticos, ou equações diferenciais do movimento, no caso de problemas dinâmicos. Em diversas situações, notadamente no caso de sistemas mecânicos complexos, as equações de equilíbrio ou do movimento podem ser obtidas de forma mais simples e elegante, em comparação com os procedimentos baseados nos princípios da Mecânica Newtoniana. Entretanto, a resolução do problema completo requer a resolução destas equações, a qual pode ser feita através de métodos analíticos apenas nos casos mais simples. Este fato justifica o uso de métodos que forneçam soluções aproximadas dos problemas estudados, conhecidos como métodos diretos do Cálculo Variacional. O método de Rayleigh-Ritz, também conhecido como Método de Ritz é um método destinado à obtenção de soluções aproximadas de problemas regidos por princípios variacionais. Ilustraremos seus princípios considerando um problema regido por um funcional do tipo: dxxu,,xu,xu,xFuI x x n 2 1 (4.96) com as condições de contorno: 11 uxu ; 22 uxu D.A. Rade Princípios Variacionais da Mecânica 94 11 uxu ; 22 uxu (4.97) 1111 nn uxu ; 1221 nn uxu Diante da impossibilidade (ou não conveniência) de se encontrar a solução exata da função u x que extremiza o funcional (4.96), propõe-se uma solução aproximada expressa da seguinte forma: 1 n i i i u x c x (4.98) As funções , 1, 2, ,i x i n são escolhidas arbitrariamente e devem ser tais que pelo menos as condições de contorno geométricas do problema sejam satisfeitas por u x . Estas funções são denominadas funções de comparação (trial functions). Note-se que, na aproximação (4.98), permanecem desconhecidas as constantes , 1, 2, ,ic i n . Introduzindo (4.98) em (4.96), e efetuando as operações de derivação e integração indicadas, obtém-se o funcional sob a forma: 1 2, , , nI u I c c c Observa-se que este funcional resulta expresso como uma função escalar de n variáveis, para o qual as condições de estacionariedade são: 1 2 0; 0; ; 0 n I u I u I u c c c (4.99) As equações (4.99) representam um sistema de n equações algébricas com n incógnitas que, uma vez resolvido, fornece os valores de , 1, 2, ,ic i n . Introduzindo estes valores em (4.98) tem-se a solução aproximada do problema. Conforme veremos no próximo capítulo, o Método de Rayleigh-Ritz é a base para a formulação do Método de Elementos Finitos pelo processo variacional. D.A. Rade Princípios Variacionais da Mecânica 95 Exemplo 4.4. Ilustraremos o método através do problema de flexão de uma viga de Euler-Bernoulli, uniforme, de módulo de rigidez à flexão EI, submetida a um carregamento transversal uniformemente distribuído q, conforme ilustrado na Figura 4.5. Desejamos obter uma aproximação para o campo de deslocamentos transversais da viga, indicado por y(x). Figura 4.5 O problema de equilíbrio é regido pelo Princípio da Energia Potencial Mínima, que, para o problema em questão, é formulado sob a seguinte forma variacional: 22 2 0 1 0 2 L d yI y EI qy dx dx (a) Vamos admitir uma solução aproximada do problema sob a forma: 1 2 2x xy x a sen a senL L (b) ou 1 1 2 2y x a x a x (c) onde: 1 xx sen L 2 2 xx sen L são denominadas funções de Ritz. Estas funções são escolhidas de modo que y x satisfaça às condições de contorno geométricas do problema. 0 0y y L y(x) L q EI x D.A. Rade Princípios Variacionais da Mecânica 96 Visando formular a condição de estacionariedade (a), vamos primeiramente avaliar: 22 2 0 2 L EI d yI y qy dx dx 22 2 1 2 1 2 0 2 2 2 2 L EI x x x xa sen a sen q a sen a sen dx L L L L L L Desenvolvendo a expressão acima e efetuando as integrações, temos: 4 42 2 2 21 2 0 2 2 2 L EI x xI y a sen a sen L L L L 2 2 1 2 1 2 2 2 22 x x x xa a sen sen qa sen qa sen dx L L L L L L Efetuando as integrações indicadas na equação acima, obtemos: 4 42 2 11 22 22 2 2 2 a LEI L EI LI y a a q L L (d) Devemos agora impor a condição de estacionariedade do funcional I. Observamos em (4.99), que ele é função de duas incógnitas a1 e a2. Assim, a condição de estacionariedade conduz a: 1 2 1 2 0I II a a a a 1 0 I a ; 2 0I a 1 0I a 4 1 5 4 qla EI 2 0I a 2 0a Assim, o campo de deflexões aproximado fica: 454 qL xy x senEI L (4.100) As soluções exata e aproximada para o deslocamento no ponto médio da viga são: D.A. Rade Princípios Variacionais da Mecânica 97 45 2 384 exato L qLy EI 4 . 5 4 2 aprox L qLy EI O erro da aproximação é de apenas 0,39 %, o que mostra a boa precisão. 4.9 Bibliografia FOWLES, G.R. CASSIDAY, G.L., Analytical Mechanics, 6th Edition, Harcourt Brace College Publishers, 1998. WASHIZU, K., Variational Methods in Elasticity and Plasticity (Monographs in Aeronautics & Astronautics), 2nd Ed., Elsevier, 1975. LEMOS, N., Mecânica Analítica, 2ª Edição, Editora Livraria da Física, 2007. LANCZOS, C., The Variational Principles of Mechanics, 4th Edition, University of Toronto Press, 1977. D.A. Rade Formulação do MEF pelo Processo Variacional 98 CAPÍTULO 5 FORMULAÇÃO DO MEF PELO PROCESSO VARIACIONAL. APLICAÇÃO À MODELAGEM DO COMPORTAMENTO DINÂMICO DE VIGAS DE EULER –BERNOULLI 5.1 – Introdução Neste capítulo, o método dos elementos finitos é formulado para obtenção das equações do movimento para vigas de Euler-Bernoulli, utilizando o processo variacional enfocado no Capítulo 3. São admitidas as seguintes hipóteses: 1ª. Consideram-se exclusivamente os efeitos da flexão transversal, não sendo considerados os efeitos associados ao movimento longitudinal; 2ª. São adotadas as simplificações da teoria de vigas de Euler-Bernoulli: desprezam-se os efeitos do cisalhamento transversal e de inércia de rotação das seções transversais; 3ª. Os deslocamentos e rotações são considerados pequenos e os materiais apresentam comportamento linear elástico linear isotrópico. 4ª. Desprezam-se os efeitos dissipativos, isto é, considera-se o sistema sem amortecimento. Será considerada a viga mostrada na Figura 5.1 de comprimento e com condições de contorno arbitrárias, para a qual se admitem as hipóteses da teoria de Euler-Bernoulli, sendo m x (kg/m) e E(x) (N.m-2) a densidade linear e módulo de elasticidade longitudinal do material que constitui a viga, respectivamente, I x (m4) é o momento de inércia de área da seção transversal da viga em relação ao eixo baricêntrico z, q(x,t) (N.m-1) é o carregamento transversal distribuído aplicado externamente e v(x,t) (m) indica o campo de deslocamentos transversais da viga. D.A. Rade Formulação do MEF pelo Processo Variacional 99 Figura 5.1 Serão a seguir detalhadas as etapas típicas de uma análise por elementos finitos, a saber: 1ª. Obtenção das equações do movimento em nível elementar, empregando as o Princípio Variaciona de Hamilton; 2ª. Obtenção das equações do movimento em nível global; 3ª. Imposição das condições de contorno; 4ª. Resolução das equações do movimento para obtenção de diferentes tipos de respostas dinâmicas. 5.2 – Obtenção das equações do movimento em nível elementar A viga mostrada na Figura 5.1 é admitida dividida em um certo número arbitrário de elementos cada um deles delimitado por dois nós. A Figura 5.2 mostra um destes elementos, identificado genericamente pelo índice i, contendo dois nós iL e iR , cujas propriedades relevantes são indicadas por i , im x , iE x , iI x . Além disso ,iq x t e ,iv x t indicam o carregamento externo e o campo de deslocamentos transversais no interior do elemento ( 0 ix ). As forças transversais ,L Ri iV t V t e os momentos ,L Ri iM t M t representam as ações mecânicas exercidas pelos elementos vizinhos, à esquerda e à direita, sobre o elemento considerado. Também são indicados os graus de liberdade escolhidos para o elemento, que são constituídos pelos valores dos deslocamentos transversais nos nós, indicados por Liv t e Riv t , e v(x,t) q(x,t) y x L D.A. Rade Formulação do MEF pelo Processo Variacional 100 pelas rotações das seções transversais nos nós, indicadas por Li t e Ri t . Estes graus de liberdade são agrupados no vetor: Te L L R Ri i i i iu t v t t v t t (5.1) onde o superscrito (e) designa os graus de liberdade em nível elementar. Deve-se observar que, de acordo com a teoria de Euler-Bernoulli, vale a seguinte relação entre as rotações das seções transversais e os deslocamentos transversais da viga: ,, ii v x tx t x (5.2) Figura 5.2 Os deslocamentos transversais são interpolados pelo seguinte polinômio cúbico: 2 31 2 3 4, 0i iv x t c t c t x c t x c t x x (5.3) Em conseqüência, levando em consideração a Eq. (5.2), as rotações das seções transversais ficam aproximadas pela seguinte função quadrática: 22 3 4, 2 3 0i ix t c t c t x c t x x (5.4) As equações (5.3) e (5.4) devem satisfazer as condições cinemáticas: 0, Ei iv t v t , Di i iv t v t (5.5.a) y x Ei Di Riv t Ri t Li t Liv t Ei(x), Ii(x), i(x) ,i x t ,iv x t qi(x) RiV t LiV t LiM t RiM t D.A. Rade Formulação do MEF pelo Processo Variacional 101 0 Ei i x v t x i Di i x v t x (5.5.b) Impondo as condições (5.5) às aproximações (5.3) e (5.4), são obtidas quatro equações algébricas lineares que, uma vez resolvidas, fornecem as expressões das constantes ic t em função dos graus de liberdade que formam o vetor dado por (5.1) e o campo de deslocamentos transversais resulta expresso da seguinte forma: 1 2 3 4, E E D Di i i i iv x t x v t x t x v t x t (5.6.a) ou: , ei iv x t x u t (5.6.b) com: 1 2 3 4x x x x x (5.7) 2 31 1 3 2x xx l l (5.8.a) 2 32 2 x xx x l ll l (5.8.b) 2 33 3 2x xx l l (5.8.c) 2 34 x xx l ll l (5.8.d) As funções , i = 1 a 4 são chamadas funções de forma ou funções de interpolação. São mostradas graficamente na Figura 5.3. D.A. Rade Formulação do MEF pelo Processo Variacional 102 Figura 5.3 As equações (5.6) mostram que o campo de deslocamentos transversais e rotações das seções transversais são expressos como combinações lineares dos deslocamentos e rotações nodais, sendo os coeficientes de combinação linear dados pelas funções de forma. Visando obter as equações do movimento em nível elementar através do Princípio Variacional de Hamilton, devemos formular as energias cinética e potencial, além do trabalho virtual do carregamento externo, conforme desenvolvimento a seguir. Energia Cinética A energia cinética do elemento é dada pela expressão: 2 0 ,1 2 i i i i v x t T m x dx t (5.9) 1(x) x 1 4(x) x 45o 45o 3(x) x 1 2(x) x i i i i D.A. Rade Formulação do MEF pelo Processo Variacional 103 Introduzindo (5.6.b) em (5.9) após derivação de iv t em relação ao tempo, escrevemos: 12 Te e ei i i iT u t M u t (5.10) onde: 4 4 0 i Te i iM m x x x dx (5.11) é a chamada matriz de massa ou matriz de inércia do elemento de viga de Euler- Bernoulli, cujo elemento geral é calculado segundo: 0 i e i i m nmn M m x x x dx , , 1a 4m n (5.12) Para o caso em que im x é constante, as equações (5.11) e (5.12) conduzem à seguinte matriz: 2 2 2 156 22 54 13 4 13 3 156 22420 4 i i e i i ii i i i i mM simétrica (5.13) Energia Potencial Para a viga de Euler-Bernoulli, a energia potencial, incluindo a energia de deformação e os trabalhos das forças e momentos externos aplicados, é expressa segundo: 2 2 2 0 ,1 2 i i i i i v x t V E x I x dx x 0 , ,i i iq x t v x t dx E E E E D D E Di i i i i i i iV t u t M t t V t u t M t t D.A. Rade Formulação do MEF pelo Processo Variacional 104 Introduzindo (5.6.b) em (5.14), após derivação dupla de ,iv x t em relação a x, escrevemos: 12 Te e ei i i iV u t K u t Te ei iu t Q t Te ei iu t Q t (5.15) onde: 4 4 0 i Te i i iK E x I x x x dx , , 1a 4m n (5.16) é a chamada matriz de rigidez do elemento de viga de Euler-Bernoulli, cujo elemento geral é calculado segundo: 0 i e i i i m nmn K E x I x x x dx , , 1a 4m n , (5.17) Além disso: 0 ,i Tei iQ t x q x t dx (5.18) é o vetor de esforços generalizados em nível elementar, cujo elemento geral é dado por: 0 ,iei m imQ t x q x t dx , 1a 4m , (5.19.a) e Te E E D Di i i i imQ t V t t V t t (5.19.b) Para o caso em que o módulo de rigidez à flexão e o carregamento distribuído são independentes de x no interior do elemento: i i i iE x I x E I , ,i iq x t q t , as equações (5.16) e (5.18) conduzem às seguintes expressões: D.A. Rade Formulação do MEF pelo Processo Variacional 105 2 2 3 2 12 6 12 6 4 6 2 12 6 4 i i e i i ii i i i i i E IK simétrica (5.20) 2 2 2 12 2 12 i i i i e i i i i i q t q t Q t q t q t (5.21) O Princípio Variacional de Hamilton estabelece: 2 1 0 t i it T V dt (5.22) Combinando as equações (5.10), (5.15) e (5.22), temos: 2 1 1 1 0 2 2 T T T Tt e e e e e e e e e e i i i i i i i i i it u t M u t u t K u t u t Q t u t Q t dt Na equação acima reconhecemos um funcional do tipo dado pela equação (3.51), no qual a variável independente é o tempo, as variáveis dependentes são os graus de liberdade elementares que formam o vetor eiu t e as derivadas de mais alta ordem são as derivadas de primeira ordem em relação ao tempo. As equações de Euler-Lagrange associadas ao funcional, dadas por (3.63), podem ser expressas, de forma alternativa, segundo: 0i ie ei iL Lddtu u (5.23) sendo i i iL T V o Lagrangeano. Efetuando as operações indicadas em (5.23), obtemos as seguintes equações do movimento em nível elementar: D.A. Rade Formulação do MEF pelo Processo Variacional 106 ( ) ( )e e e ee ei i i i i iM u t K u t Q t Q t (5.24) 5.2 Montagem das equações do movimento em nível global Na seção anterior, as equações do movimento foram obtidas para um elemento genérico considerado isolado. Entretanto, deve ser considerado que a viga é discretizada em vários elementos de modo que elementos vizinhos compartilham nós, conforme ilustrado na Figura 5.4 para dois elementos contíguos i e j. Te L L R Ri i i i iu t v t t v t t Te L L R Rj j j j ju t v t t v t t 1 1 2 2 3 3 Tgu t v t t v t t v t t Figura 5.4 graus de liberdade elementares Li Ri Riv t Ri t Li t Liv t Lj Rj Rjv t Rj t Lj t Ljv t 3 3v t 3 t 1 2 2v t 2 t 1 t 1v t graus de liberdade globais D.A. Rade Formulação do MEF pelo Processo Variacional 107 As equações do movimento dos dois elementos devem ser combinadas de modo que seja assegurado o equilíbrio dinâmico e a continuidade de deslocamentos e rotações nos nós que são compartilhados dois elementos vizinhos. Este procedimento será ilustrado a seguir para os dois elementos mostrados na Figura 5.4. Estabelecem-se, para cada elemento, transformações lineares relacionando os graus de liberdade elementares e os graus de liberdade globais, sob a forma: 4 64 1 6 1e gi iu t L u t (5.25.a) 4 64 1 6 1e gj ju t L u t (5.25.b) O esquema da Figura 5.4 permite facilmente construir as matrizes iL e jL , conforme mostrado abaixo. Nota-se que são matrizes booleanas (formadas por zeros e uns) nas quais as posições dos uns permitem localizar os graus de liberdade elementares dentre os graus de liberdade globais. Elemento i 1 1 2 2 3 3 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 L i L i R i R i v t v t t t v t v t t t v t t (5.26.a) Elemento j 1 1 2 2 3 3 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 L j L j R j R j v t v t t t v t v t t t v t t (5.26.b) Para o sistema formado pelos elementos i e j escreve-se a energia cinética total sob a forma: D.A. Rade Formulação do MEF pelo Processo Variacional 108 i jT T T (5.27) Utilizando a Eq. (5.10) para ambos os elementos e as equações (5.25), obtém- se a energia cinética total expressa sob a seguinte forma, em termos dos graus de liberdade globais: 12 Tg g gT u t M u t (5.28) onde a matriz de massa global é expressa segundo: 6 6 TTg e e i i i j j jM L M L L M L (5.39) De modo semelhante, opera-se sobre a energia potencial total e o trabalho virtual total do carregamento externo: i jV V V = 12 T Tg g g g gu t K u t u t Q t (5.30) com: 6 6 TTg e e i i i j j jK L K L L K L (5.31) e: 6 1 TTg e ei i j jQ L Q L Q (5.32) Vale observar que, de acordo com a terceira Lei de Newton, as forças e momentos exercidos pelo elemento 1 sobre o elemento 2 são iguais aos correspondentes exercidos pelo elemento 2 sobre o elemento 1. Isso faz com que as operações que conduzem a (5.30) eliminem tais esforços correspondentes aos nós intermediários da formulação, permanecendo apenas aqueles correspondentes aos nós das extremidades da viga. Estes últimos são tratados quando são impostas as condições de contorno, conforme detalhado na Seção 4.3. Aplicando novamente o princípio de Hamilton expresso segundo: 2 1 0 t i it T V dt (5.33) obtemos as seguintes equações de Euler-Lagrange associadas: D.A. Rade Formulação do MEF pelo Processo Variacional 109 gg gi L d L Q dtu u com L T V (5.37) Associando as equações (5.28), (5.30) e (5.37), obtêm-se as equações do movimento em nível global, formado pelos dois elementos finitos conectados: g g g g gM u t K u t Q t (5.38) As equações (5.39) constituem um sistema de 6 equações diferenciais ordinárias no tempo que possuem exatamente a mesma forma das equações do movimento de sistemas discretos não amortecidos. Desta forma, os métodos apresentados naquele capítulo podem ser utilizados para calcular os diferentes tipos de respostas dinâmicas do sistema sujeito a diferentes formas de excitação, que podem ser completadas por um conjunto de condições iniciais: 00g gu u 00g gu u (5.39) O procedimento aqui ilustrado para o acoplamento de apenas dois elementos pode ser aplicado sucessivamente para o acoplamento de um número qualquer de elementos. 4.3 Imposição das condições do contorno Ao contrário do que ocorre nos problemas de equilíbrio (problemas estáticos) os problemas de análise dinâmica podem ser resolvidos sem dificuldades particulares quando o sistema estrutural estiver livre no espaço, sem nenhuma vinculação ao sistema de referência fixo. Este é o caso de estruturas aeronáuticas ou espaciais. Nestes casos, o sistema de equações diferenciais (5.38) pode ser processado diretamente. Por outro lado, caso haja vinculações cinemáticas representadas por deslocamentos impostos em alguns pontos da estrutura, as equações do movimento devem ser devidamente modificadas para levar em conta tais vinculações. Apresenta-se, a seguir um procedimento para esta finalidade. As matrizes e vetores que figuram em (5.38) são particionados segundo os graus de liberdades ditos “livres”, que não são afetados pelos vínculos cinemáticos e os graus de liberdade denominados “impostos” que são aqueles cujos valores são prescritos. Os primeiros são doravante indicados por subscritos “L” e os últimos por subscritos “I” nas equações do movimento, que se expressam segundo: g g gg g g g L L LLL LI LL LI g g g gg g g IL II IL III I I u t u t Q tM M K K M M K Ku t u t Q t (5.40) D.A. Rade Formulação do MEF pelo Processo Variacional 110 Os dois blocos de equações resultantes do particionamento indicado em (5.40) podem ser postos sob a forma: g g g g g g g g gLL L LL L L LI I LI IM u t K u t Q t M u t K u t (5.41) g g g g g g g g gI IL L II I IL L II IQ t M u t M u t K u t K u t (5.42) A resolução das equações (5.41) e (5.42) permite determinar a resposta dinâmica nas coordenadas livres e, subseqüentemente, os esforços de reação associados aos graus de liberdade impostos. 4.4 Realização de análise dinâmicas A partir das equações do movimento globais resolução das equações (5.41) e (5.42), podem-se computar diferentes tipos resposta dinâmica, cuja formulação é apresentada a seguir: Análise Modal Para o problema homogêneo (respostas em regime livre): 0g g g gLL L LL LM u t K u t (5.43) procuram-se soluções harmônicas do tipo: g i tLu t U e (5.44) Introduzindo (5.44) em (5.43), chega-se ao seguinte problema de autovalor: 2 0g gLL i LL iK M U (5.45) cujas soluções ,i iU , i=1,2, ...são, respectivamente, as freqüências naturais e os modos naturais de vibração. Análise Harmônica Considera-se o problema não-homogêneo (respostas em regime forçado): g g g g gLL L LL L LM u t K u t Q t (5.43) com as forças excitadoras e respostas harmônicas dadas, respectivamente, por: D.A. Rade Formulação do MEF pelo Processo Variacional 111 g i tL LQ t Q e (5.44) g i tL Lu t U e (5.45) Introduzindo (5.44) e (5.45) em (5.43), obtêm-se a seguinte relação entre as amplitudes de excitação e amplitudes das respostas harmônicas: 2g gLL LL L LK M U Q (5.46) ou: 12g gL LL LL LU K M Q (5.47) Análise Transiente O problema de análise transiente consiste na determinação das respostas temporais da estrutura sob a ação de forças externas variáveis no tempo e/ou a um conjunto de condições iniciais (deslocamentos e/ou velocidades) impostas no instante inicial. Alternativamente, a excitação pode ser provocada por deslocamentos impostos variáveis no tempo. Nestes casos, a resolução é feita integrando as equações do movimento (5.41), utilizando algoritmos de integração numérica, tais como os métodos de Newmark, Runge-Kutta, Huboldt, dentre outros existentes para esta finalidade. Ao final da integração, são conhecidos os valores dos graus de liberdade do modelo de elementos finitos em um conjunto discreto de pontos dentro do intervalo de tempo escolhido. 4.5. Bibliografia CRAIG JR., R. R. e KURDILA, A. J. Fundamentals of Structural Dynamics. 2.ed., New Jersey – USA: John Wiley & Sons, 2006. 744p. MEIROVITCH, L. Methods of Analytical Dynamics. New York – USA: McGraw-Hill, 1970. 524p.