Elementos Finitos - Apostilas - Engenharia Mecânica Part3, Notas de estudo de Engenharia Mecânica

Baixe Elementos Finitos - Apostilas - Engenharia Mecânica Part3 e outras Notas de estudo em PDF para Engenharia Mecânica, somente na Docsity! 75 Apêndice E Revisão Conceitual de Equações Diferenciais Parciais O conjunto de equações a derivadas parciais que regem os fenômenos físicos é bem vasto, e compreende um número de casos particulares que seria ilusório querer descrever de maneira exaustiva, sem ir de encontro a uma redação enciclopédia. Pode-se, no entanto, classificar a maioria das equações em três grandes classes, cada uma ilustrada por um tipo de fenômeno bem particular. Seja a u(x, y) e equação diferencial de segunda ordem abaixo: 0 2 22 2 2 =++ ∂ ∂+ ∂ ∂+ ∂ ∂+ ∂∂ ∂+ ∂ ∂ GFu y u E x u D y u C yx u B x u A - Se B2 – 4 AC > 0, têm-se as Equações Hiperbólicas, como, por exemplo, as equações de ondas: . ),( 2 2 2 2 x u t u temostxuPara ∂ ∂= ∂ ∂ - Se B2 – 4 AC < 0, têm-se as Equações Elípticas, como nas equações de Laplace ou de Equilíbrio: .0 ),( 2 2 2 2 = ∂ ∂+ ∂ ∂ y u x u temosyxuPara - Se B2 – 4 AC = 0, têm-se as Equações Parabólicas, exemplificado pelas equações de condução de calor: . ),( 2 2 x u t u temostxuPara ∂ ∂= ∂ ∂ As equações do tipo elípticas são representativas dos problemas de potencial que aparecem nos estudos em regime permanente na eletricidade (eletrostática ou magnetostática), mecânica (deformação de um sólido, escoamento Laplaciano de um fluido) e térmica (distribuição de temperaturas). As condições de contorno normalmente associadas são do tipo: Dirichlet, ))()(( 00 sfusu == Neuman
     = ∂ ∂ )()( 0 sfsn u ou mista
     = ∂ ∂+ )()()( 0 sfsn u su As equações parabólicas são representativas dos problemas de difusão, cuja equação de difusão da temperatura em um corpo incompreensível é o caso típico. Esta equação é igual a de penetração das correntes induzidas em um corpo condutor de eletricidade. As condições de contorno associadas à equação parabólica são de dois tipos: a) Condições de Dirichlet, Neuman ou mista sobre a fronteira espacial do domínio. b) Uma condição inicial (t=0) em todo o domínio. 76 As equações hiperbólicas caracterizam o fenômeno da propagação das ondas, sejam elas vibratórias do tipo mecânicas ou eletromagnéticas. As condições de contorno relacionadas à equação de propagação são associadas às condições de Cauchy ao instante inicial (dada sua função u e sua derivada ∂u/∂t em relação ao tempo). Desde que o domínio de estudo seja formado de vários sub-domínios, sobre os quais os coeficientes da equação diferem por causa das propriedades físicas, diferentes materiais que compõem estes sub-domínios, há na fronteira de separações derivadas de ordens elevadas e formação de condições de transmissão que exprimem as condições atribuídas às diversas funções e suas derivadas. Como exemplo, na eletrostática é bem conhecido que na passagem de um meio isolante a outro, há refração das linhas de campo elétrico. Chamando de V o potencial elétrico definido no interior de um domínio D formado de um sub-domínio D1 , no qual o material dielétrico tenha permissividade dielétrica ε1 , e um sub-domínio D2 , caracterizado por um meio de permissividade ε2 , a equação do potencial varia de um meio a outro. Meio 1 : 011 =

     ∂ ∂ ∂ ∂+
     ∂ ∂ ∂ ∂ y V yx V x εε Meio 2 : 022 =

     ∂ ∂ ∂ ∂+
     ∂ ∂ ∂ ∂ y V yx V x εε Na interface, estas equações não são válidas e substituem-se por duas outras: da continuidade do potencial e da indução elétrica, que leva em conta a descontinuidade do campo elétrico: Interface: V1 = V2 e n V n V ∂ ∂= ∂ ∂ 2 2 1 1 εε (n normal a interface) Os problemas elípticos são características da análise de fenômenos de regime permanente, fenômenos do tipo estático (sem variação temporal) ou variável no tempo segundo função conhecida (senoidal, por exemplo). Os problemas parabólicos e hiperbólicos são ligados aos estudos de regime transitório (chamados algumas vezes de dinâmicos) e sua resolução permite analisar a evolução de um fenômeno físico no decorrer do tempo (regime transitório elétrico ou térmico, resposta mecânica a uma perturbação). Noção de um problema “Bem Definido” As equações a derivadas parciais e as condições de contorno associadas constituem o chamado problema a derivadas parciais. Entretanto, a associação de uma equação e das condições de contorno não conduzem forçosamente a um problema matemático simulando um fenômeno ou um processo físico. Com finalidade de caracterizar mais precisamente este tipo de problema, Hadamard introduziu a noção de problema “bem definido”. Um problema é “bem definido” desde que satisfaça as três condições seguintes: - O problema possui uma solução; 79 Apêndice F Visão Histórica Esta parte é baseada na contribuição original de [34]. O desenvolvimento moderno do método de elementos finitos começa na década de quarenta no campo da engenharia estrutural com o trabalho de HRENNIKOFF [1] (1941) e MCHENRY [2] (1943), usando um modelo de elementos em linhas, de uma dimensão tipo barras e feixes, para a solução de problemas de esforços em sólidos contínuos. Em artigos publicados em 1943, mas não largamente reconhecidos por muitos anos, COURANT [3] propôs usar a solução para os esforços, de uma forma variacional. Ele. então, introduz funções de interpolação sobre pedaços de sub-regiões triangulares, e une todas estas em um método, para obter uma solução numérica aproximada. LEVY [4] (1947), desenvolve o método de flexibilidade ou forças, e, em outro trabalho [5] (1943), sugere um método de rigidez ou deslocamentos, prometendo ser uma alternativa na análise de estruturas redundantes estáticas de aviões. Contudo, suas equações apresentavam dificuldade de solução manual, e seu método só se tornaria popular com o advento dos computadores digitais de alta velocidade. ARGYRIS e KELSEY [6], [7] (1954) desenvolvem um método matricial de análise estrutural, usando o princípio das energias. Este ilustraria a importância das regras, que o princípio das energias poderia ter no Método dos Elementos Finitos. O primeiro tratamento de elementos em duas dimensões foi dado por TURNER, CLOUGH, MARTIN e TOPP [8], em 1956. Eles derivaram matrizes de rigidez de elementos de pontos, feixes e elementos triangulares e retangulares de duas dimensões em planos de esforços, e esboçaram procedimentos comuns conhecidos, como o método da rigidez direta para obter a matriz total da estrutura rígida. Logo, com desenvolvimento de computadores digitais de alta velocidade, no começo da década de cinqüenta, o trabalho de TURNER, CLOUGH, MARTIN e TOPP [8] apontou para o desenvolvimento das equações de rigidez de elementos finitos expressas na notação matricial. O termo “Elementos Finitos” foi introduzida por CLOUGH [9] (1960), quando ambos, elementos triangulares e retangulares, foram usados na análise do plano de esforços. SILVERTER [58], nesta mesma época, cria avanços no modelamento matemático, trabalhando em esquemas de discretização e ajuste de curvas. A matriz de rigidez para elementos curvos de planos retangulares foi desenvolvida por MELOSH [10], em 1961, seguida pelo desenvolvimento de conchas (três dimensões) curvas em elementos curvos de matriz de rigidez, para conchas assimétricas e pressões em embarcações, por GRAFTON e STROME [11] (1963). Extensões do Método dos Elementos Finitos, para os problemas de três dimensões, com o desenvolvimento da matriz de rigidez tetrahédrica, foram realizadas por MARTIN [12] (1961), por GALLAGHER, PADLOG e BIJLAARD [13] (1962) e por MELOSH [14] (1963). Elementos tridimensionais adicionais foram estudados por ARGYRIS [15] (1964). Os casos especiais de sólidos assimétricos foram considerados por CLOUGH, RASHID [16] e WILSON [17], em 1965. A maioria dos trabalhos de elementos finitos, no começo dos anos 60, lidava com pequenos esforços e deslocamentos, elasticidade de materiais e cargas estáticas. Contudo, grandes deflexões e análises térmicas foram consideradas por TURNER, DILL, MARTIN e MELOSH [18] (1960), e a não linearidade dos materiais, por GALLAGHER, 80 PADLOG, e BIJLAARD [13] (1962), quando problemas de curvaturas foram inicialmente tratados por GALLAGHER e PADLOG [19] (1963). ZIENKIEWICZ, WATSON e KING [20] (1968) estenderam o método para problemas de visco-elasticidade. Em 1965, ARCHER [21] considera a análise dinâmica no desenvolvimento da matriz de massas de consistência (distribuições das cargas), aplicada para análise de sistemas de massa distribuída, tal como barras e feixes na análise estrutural. MELOSH [14] considera, em 1963, que o método de elementos finitos poderia ser colocado em termos da formulação variacional, e então, este começou a ser usado para a resolução de aplicações não estruturais. Problemas de cálculo de campo, tal como a determinação de torção de hastes, escoamento de fluidos e condição do calor foram resolvidos por ZIENKIEWICZ e CHEUNG [22] (1965), MARTIN [23] (1968), e WILSON e NICKEL [24] (1966). O trabalho que marcou a aplicação do MEF na Engenharia de eletricidade é creditado a SILVERTER & CHARI [74], em 1969, que promoveram seu desenvolvimento. A partir de então, uma série de pesquisadores dedicou esforços, no sentido de aplicá-lo na resolução dos maiores problemas da Engenharia de eletricidade, que é o cálculo de campos eletromagnéticos presentes nos dispositivos e sistemas elétricos. As extensões dos métodos foram possíveis, pela adaptação da técnica dos resíduos ponderados, primeiramente, por SZABO & LEE [25] (1969), que derivaram as equações de elasticidade usadas em análise estrutural, e posteriormente, por ZIENKIEWICZ & PAREKH [26] (1970), que a utilizaram nos problemas de campos transitórios. Entende-se, então, que a formulação direta e a variacional era difícil ou impossível de se usar, e logo, o método dos resíduos ponderados tornava-se mais apropriado. Por exemplo, em 1977, LYNESS, OWEN & ZIENKIEWICZ [27] aplicaram os resíduo ponderados para determinação de campos magnéticos. BELYTSCHKO [28], [29] (1976) estudou problemas associados a grandes deslocamentos em dinâmica não linear, e melhorou as técnicas numéricas para a resolução de sistemas de equações. Aplicações em novos campos como a bioengenharia foram tratadas pelo Métodos dos Elementos Finitos, embora outros, como materiais não lineares, geometria não lineares e outras complexidades, permaneceram ainda sem solução. Resumidamente, o Método dos Elementos Finitos é conhecido, no que se refere aos seus princípios, a, aproximadamente, meio século, e só veio a ser utilizado, com o aparecimento dos modernos meios informáticos. De fato, apesar das formulações integrais serem conhecidas a longo tempo, graças aos trabalhos de GALERKINE, RITZ, COURANT e HILBERT, suas aplicações não puderam se generalizar, de imediato, pela dificuldade em se resolver, facilmente, os sistemas algébricos lineares e não lineares de grandes dimensões. De fato, o trabalho necessário para se resolver um sistema linear, com algumas dezenas de equações e com o correspondente número de incógnitas, levou a maioria dos engenheiros a restringir este tipo de cálculo a um pequeno número de especialistas, como CHOLESK, SEITEL, JACOBI e GAUSS, que desenvolveram métodos astuciosos, recentemente, ainda utilizados. Os engenheiros civis e mecânicos, enfrentando problemas complexos de cálculo estrutural, foram os primeiros a se aproveitarem do desenvolvimento da informática e das 81 linguagens de programação de alto nível, para expressar, em termos de equações algébricas, os modelos de desempenho das estruturas mecânicas. Nos anos sessenta, o Método dos Elementos Finitos ganha rigor matemático e implementação computacional, passando a ser adotado por profissionais de diversas áreas, na resolução de equações diferenciais parciais. A partir de 1970, sob o impulso de numerosos pesquisadores da Engenharia e Matemática (ZIENKIEWICZ, GALLAGHER, ODEN, LIONS, RAVIART, SILVESTER, CHARI, TOUZOT e outros), este método tornou-se realmente popular entre os engenheiros de todas as especialidades, e sua utilização implicou no aparecimento de vários softwares comerciais, utilizados na mecânica (NASTRAN, ASKA), em problemas de temperatura (TITUS) e eletromagnetismo (FLUX, MAGNET11, PE2D, MATLAB), para citar apenas as aplicações mais desenvolvidas. No Brasil, o Método dos Elementos Finitos chega, na década de 70, por meio de inúmeras dissertações de mestrado e teses de doutorado na Universidade Federal do Rio de Janeiro (UFRJ), Universidade Estadual de São Paulo (USP), Universidade Federal de Santa Catarina (UFSC) e Universidade Federal de Minas Gerais (UFMG), originando a criação de disciplinas de pós-graduação sobre o tema. Em âmbito nacional, o primeiro trabalho sobre a aplicação deste método na engenharia de eletricidade foi desenvolvido na Escola Politécnica da USP, por JANISZEWSKI [75], em 1978, época em que surgem alguns grupos de pesquisa em universidades brasileiras, destacando-se o GRUCAD da UFSC, o grupo de pesquisa da UFMG e a equipe de Simulação de Fenômenos Eletromagnéticos da Escola Politécnica da USP. Em seqüência, na década de noventa, surgem cursos de pós-graduação ministrado por empresas comerciais, para suprir a grande necessidade de atualização no mercado de trabalho da Engenharia. E finalmente, após longos anos sendo ministrado apenas em cursos de pós-graduação, em 1995, a USP cria a disciplina optativa de Método dos Elementos Finitos, para o curso de graduação em Engenharia Civil, tornando-se pioneira no ensino de graduação, e assumida por diversos outros cursos de graduação em todo o país. As universidades brasileiras, a exemplo do que ocorreu nas universidades dos países desenvolvidos, foram as responsáveis pelo lançamento, no mercado nacional, dos primeiros produtos, visando atender as necessidades da comunidade acadêmica, no sentido de suprir as pesquisas neste setor, bem como assessorar o setor industrial em suas necessidades de estudos e projetos. Ferramentas computacionais são criadas, no Brasil, pelo produto das pesquisas técnico-científicas dos cursos de mestrado e doutorado realizadas nas universidades que implementam o Método dos Elementos Finitos. Muitas dessas pesquisas são convertidas em softwares comerciais e integradas a programas de desenho de CAD, gerando um novo tipo de ferramenta: a CAE – Computer Aided Enginneering. Esta, por sua vez, parece representar o futuro do projeto e concepção da Engenharia moderna. Outros métodos numéricos, derivados dos elementos finitos, surgem em diversas partes do mundo, destacando-se o método de volumes finitos e o método dos elementos de fronteira ou contorno. Em geral, eles são aplicados a problemas específicos de Engenharia e ciências físicas. Além disso, diversas variações na implementação e formulação dos elementos finitos são comuns, principalmente, aquelas voltadas para otimizações de ajustes e características particulares de uma dada aplicação. 84 dydzzfhxgxu x y



 

  +−+= 1 0 )()1()( (G.7) onde y e z são usados para denotar variáveis temporárias. Contudo, isso não é o principal interesse aqui. Está-se interessado em desenvolver um esquema para obter soluções aproximadas de (S) que serão aplicáveis para muitas situações complexas na qual a solução exata não é possível. Alguns métodos de aproximação começam diretamente com a formulação forte do problema. O exemplo mais notável é o Método das Diferenças Finitas. O Método dos Elementos Finitos requer uma formulação diferente, a qual é tratada no próximo item. G.4. Formulação Fraca ou Variacional Para definir a forma fraca ou variacional, contrapartida de (S), necessita-se caracterizar duas classes de funções. A primeira é composta de funções admissíveis. Neste conjunto de soluções admissíveis é importante impor que a condição de contorno u(1)=g seja satisfeita. A outra condição de contorno não será requerida na definição. É importante também, para que certas expressões façam sentido, impor que as derivadas primeiras das funções admissíveis tenham quadrado integrável. Isto é, se u é uma solução trivial, então
∞< 1 0 2 , )( dxu x (G.8) Funções que satisfazem (G.8) são chamadas de funções H1; escreve-se u ∈ H1. Algumas vezes o domínio é explicitado, isso é u ∈ H1([0, 1]). Então o conjunto de soluções admissíveis, denotado por γ, consiste de todas as funções as quais tem derivada primeira com quadrado integrável e que tenham o valor de g para x=1. Isso é representado como se segue: γ = { u | u ∈ H1 , u(1)=g } (Funções admissíveis) (G.9) O fato de que γ é um conjunto de funções é indicado pela chave em (G.9). A notação para um membro típico de um conjunto, neste caso u, vem em primeiro dentro no lado esquerdo das chaves, seguido à linha vertical (|) e a propriedade satisfeita pelo membro do conjunto. A segunda classe de funções é chamada de funções de teste ou funções de peso. Este conjunto é muito similar ao conjunto da solução admissível exceto que este requer a imposição de condições de contorno homogênea de g. Isto é, este conjunto requer que funções de teste, w, satisfaçam w(1)=0. O conjunto é denotado por υ e definida por υ = {w | w ∈ H1, w(1)=0 } (Funções de teste) (G.10) Isso simplifica o assunto no qual se quer ter f : Ω → R como sendo suave. 85 Em termos da definição anterior, pode-se agora estabelecer uma forma fraca apropriada, (W), do problema de valor de contorno.

  += ∈∈

1 0 1 0 ,x, )0(w w todopara que tal uAchar antes. como h, e g, f, Dando )( hwwfdxdxu W x υγ (G.11) Formulações deste tipo são freqüentemente chamadas de trabalho virtual ou deslocamento virtual na mecânica. Os w’s são os deslocamentos virtuais. A equação (G.11) é chamada de equação variacional ou (especificamente na mecânica) a equação do trabalho virtual. A solução de (W) é chamada de solução fraca ou generalizada. A definição dada da formulação fraca não é a única possível, mas é a mais natural para os problemas que serão considerados. G.5. Equivalência entre as Formas Forte e Fraca Claramente, há algum relacionamento entre a formulação forte e fraca do problema, caso contrário não existiria razão para introduzir a forma fraca. Mostra-se que a solução fraca e forte são idênticas. Isso se estabelece assumindo que todas as funções são suaves. Isso permitirá prosseguir sem envolver condições técnicas que complique a matemática envolvida. Prova deste tipo exige algumas vezes provas formais. A intenção aqui não é de apresentar uma prova completamente rigorosa, mas tornar plausível se acreditar na proposição. Com esta filosofia em mente, prova-se o seguinte: Proposição: a) Seja u uma solução de (S). Então u é também uma solução de (W). b) Seja u uma solução de (W). Então u é também uma solução de (S). Um outro resultado, o qual não se preocupou verificar, mas que de fato é facilmente estabelecido, é que ambas (S) e (W) possuem solução única. Então de (a) e (b), a solução forte e fraca são uma e a mesma. Conseqüentemente, (W) é equivalente a (S). Prova Formal: a) Como u é assumido ser uma solução de (S), pode-se escrever
+−= 1 0 , )(0 dxfuw xx (G.12) para qualquer w ∈ υ. Integrando G.12 por partes resulta em

−−= 1 0 1 0 0 1 ,,,0 xxx wuwfdxdxuw (G.13) Rearranjando e fazendo uso do fato que –u,x(0)=h e w(1)=0 resulta em 86

+= 1 0 1 0 ,, )0( hwwfdxdxuw xx (G.14) Além do mais, tem-se que u é uma solução de (S), satisfaz u(1)=g e portanto está em γ. Finalmente, desde que u também satisfaz (G.14) para todo w ∈ υ, u satisfaz a definição de uma solução fraca dado por (W). b) Agora u é assumido ser uma solução fraca. Então u ∈ γ; conseqüentemente u(1)=g, e

+= 1 0 1 0 ,, )0( hwwfdxdxuw xx para todo w ∈ υ. Integrando por partes e usando do fato de w(1)=0 resulta em
+++= 1 0 ,, ])0()[0()(0 huwdxfuw xxx (G.15) Provar que u é uma solução de (S) é suficiente mostrar que (G.15) implica em (estas equações são algumas vezes chamadas de equação de Euler-Lagrange da formulação fraca): 0)0(u . e ; em 0u . x, xx, =+ Ω=+ hii fi Primeiro, prova-se (i). Define w em (G.15) por w = φ( u,xx + f) (G.16) onde φ é suave; φ(x)>0 para todo x ∈ Ω = ]0, 1[; e φ(0) = φ(1) = 0. Por exemplo, pode-se tomar φ(x) = x(1-x), o qual satisfaz todos os requisitos estipulados. Segue que w(1) = 0 e então w ∈ υ, assim (G.16) define uma função membro legítima de υ. Substituindo (G.16) em (G.15) resulta em 0)(0 1 0 0 2 , ++=
≥ dxfu xx

 φ (G.17) Desde que φ > 0 em Ω, segue-se de (G.17) que (i) deve ser satisfeito. Agora que se tem estabelecido (i), passa-se a usar (G.15) para provar (ii), assim 0 = w(0)[u,x(0)+h] (G.18) Vê-se que w ∈ υ não impõem qualquer restrição no seu valor em x = 0. Portanto, pode-se assumir que o w em (G.18) é tal que w(0) ≠ 0. Então (ii) está também satisfeita, a qual completa a prova da proposição. 89 uh(1) = vh(1) + gh(1) = 0 + g (G.32) Assim (G.30) constitue uma definição de γh; que é, γh são todas as funções da forma de (G.30). O ponto importante a observar é que, acima disto às funções gh, γh e υh são compostas de conjuntos idênticos de funções. Agora escreve-se a equação variacional, da forma (G.21), em termos de uh ∈ γh e wh ∈ υh: a(wh,uh) = (wh,fh) + wh(0)h (G.33) Essa equação é considerada como definindo uma solução aproximada (fraca), uh. Substituindo (G.30) em (G.33), e a bilinearidade de a(.,.) possibilita escrever: a(wh,uh) = (wh,f) + wh(0)h – a(wh,gh) (G.34) O lado direito consiste totalmente de termos associados com dados fornecidos (isso é, f, g, e h). A equação (G.34) é usada para definir vh, a parte desconhecida de uh. A forma (Bubnov-) Galerkin do problema, denotada por (G) é representada como se segue:

  −+= ∈ ∈+= ),()0(),(),( w todopara que tal , vonde ,uachar antes, como h, e g, f, Dado )( h hh hhhhhh h hhh gwahwfwvwa gv G υ υ Note que (G) é justamente uma versão de (W) em termos de uma coleção de funções com dimensões finitas em υh. O método de Bubnov-Galerkin é comumente referenciado como simplesmente método de Galerkin, terminologia que se adotará para frente. A equação (G.34) é algumas vezes referida como Equação de Galerkin. O método de aproximações do tipo considerado são exemplos do chamado Método dos Resíduos Ponderados. G.8. As Equações Matriciais – A Matriz de Rigidez K O método de Galerkin conduz a um sistema de equações algébricas lineares. Para ver isso, necessita-se dar mais estrutura para a definição de υh. Assim υh consiste de toda combinação linear de funções dadas denotadas por NA:Ω→R, onde A = 1,2,....,n. Isso significa que se wh ∈ υh, então existe constantes cA, A = 1,2,...,n, tal que nn n A AA h NcNcNcNcNcw ++++==
= ...332211 1 (G.35) As NA’s são referidas com funções de forma, base ou interpolação. Requer-se que cada NA satisfaça NA(1) = 0, A = 1, 2, ...., n (G.36) 90 Da qual segue por (G.35) que wh(1)=0, como é necessário. υh é dito ter dimensões n por razões obvias. Para definir membros de γh necessita-se especificar gh. Para esse fim, introduz-se uma outra função de forma Nn+1: Ω→R, a qual tem a seguinte propriedade Nn+1(1) =1 (G.37) (Note Nn+1∉ υh.) Então gh é dado por gh = g Nn+1 (G.38) e então gh(1) = g (G.39) Com essas definições, um típico uh ∈ γh pode ser escrito como 1 1 + = +=+=
n n A AA hhh gNNdgvu (G.40) onde os dA’s são constantes e das quais é aparente que u h(1)=g. Substituindo (G.35) e (G.40) na Equação de Galerkin tem-se,
    −  +
     =
     + === ==




1 111 11 ,)0(, , n n A AA n A AA n A AA n B BB n A AA gNNcahNcfNc NdNca (G.41) Usando a bilinearidade de a(.,.) e (.,.), (G.41) torna
= = n A AAGc 1 0 (G.42) onde gNNahNfNdNNaG nAAA n B BBAA ),()0(),(),( 1 1 + = +−−=
(G.43) Agora a equação de Galerkin é garantida para todo wh ∈ υh. De (G.35), isso significa para todo cA’s, A = 1,2,...,n. Uma vez que os cA’s são arbitrários em (G.42), necessariamente segue que cada GA, A = 1,2,...,n, deve ser identicamente zero isto é de (G.43) gNNahNfNdNNa nAAA n B BBA ),()0(),(),( 1 1 + = +−=
(G.44) Note que todos os termos são conhecidos em (G.44) exceto os dB’s. Então (G.44) constitue um sistema de n equações e n incógnitas. Este pode ser escrito numa forma 91 mais concisa como se segue: Tem-se KAB = a(NA, NB) (G.45) FA=(NA, f) + NA(0)h – a(NA, Nn+1) g (G.46) Então (G.44) torna n..., 1,2, A , 1 ==
= A n B BAB FdK (G.47) Maior simplicidade é ganha adotando-se uma notação matricial. Tem-se







            == nnnn n n AB KKK KKK KKK KK ... ...... ...... ...... ... ... ][ 21 22221 11211 (G.48)







 







  == − n n A F F F F FF 1 2 1 . . . }{ (G.49) e







 







  == n B d d d dd . . . }{ 2 1 (G.50) Agora (G.47) pode ser escrita com Kd = F (G.51) As seguintes terminologias são freqüentemente aplicadas, especialmente quando o problema sobre consideração pertence a um sistema mecânico: K = Matriz de Rigidez F = Vetor de Forças d = Vetor Deslocamentos 94 Figura G.2 Funções Base para o espaço de elementos finitos lineares “Piecewise”. Um típico membro wh ∈ υh tem a forma
= n A AANc 1 e aparece como na figura G.3. Note que wh é contínua mas tem descontinuidades de declive (ou na derivada) em cada elemento de fronteira. Por essa razão, wh,x, derivada generalizada de w h, irá ser constante em pedaços (Piecewise), experimentando descontinuidades entre os elementos de contorno ou fronteira. (Tal função é chamada algumas vezes de uma função de passo generalizado.) Restringindo-se a cada elemento do domínio, wh é um polinômio linear em x. Em relação as condições de contorno essenciais homogêneas, wh(1)=0. Claramente, wh é identicamente zero se e somente se cada cA=0, A=1,2,...,n. Figura G.3 Um típico membro wh ∈ υh. Membros típicos de γh são obtidos pela adição de gh = g Nn+1 a membros típicos de υh. Isso assegura que uh(1)=g. As funções de elemento finito linear “Piecewise” são mais simples e mais largamente usadas por funções de elementos finitos em problemas unidimensionais. Exemplo: Considere a formulação fraca do problema no modelo unidimensional:

+= 1 0 1 0 ,, )0( hwwfdxdxuw xx (G.60) 95 onde w ∈ υ e u ∈ γ são assumidos ser suave nos interiores dos elementos (isso é, ]xA, xA+1[, A=1,2,... n), mas podem sofrer descontinuidades de declive entre os elementos de contorno. (Funções desta classe contém o espaço de elementos finitos lineares “Piecewise” descritos anteriormente.) De (G.60) e assumindo continuidade das funções, mostra-se que: )]()([)(])0()[0()(0 ,, 21 ,, 1 −+ == + −++++=

 + AxAx n A A n A x x xxx xuxuxwhuwdxfuw A A (G.61) Pode ser concluído que as condições de Euler-Lagrange de (G.61) são i. u,xx(x)+f(x)=0, onde x∈]xA, xA+1[ e A = 1,2, ...,n, ii. –u,x(0 +)=h; e iii. u,x(xA -)=u,x(xA +), onde A= 2,3,..., n. Observe que (i) é a equação diferencial restrita a elementos interiores, e (iii) é uma condição de continuidade entre os elementos de contorno. Isso pode ser contrastado com o caso no qual a solução é assumida suave. Neste caso a condição de continuidade é identicamente satisfeita e o somatório das integrais sobre os elementos interiores pode ser trocada por uma integral sobre todo domínio. Na formulação de elementos finitos de Galerkin, uma solução aproximada de (i)-(iii) é obtida. G.10. Propriedades da Matriz K As funções de forma NA, A= 1,2,..., n+1 são zero fora da vizinhança do “nó”. Como resultado, muitos termos de K são zero. Isso pode ser visto como se segue. Seja B > A+1. Então (Figura G.4) 0 1 0 0 ,, ==
dxNNK xBxAAB

 (G.62) NA NB A A+1 B Figura G.4 Se B > A+1, a porção diferente de zero de NB e NA não faz sobreposição. A simetria de K implica, adicionalmente, que (G.62) é garantida para A > B+1. E é dito que K é semidiagonal (isso é, seus valores diferentes de zero estão em uma banda sobre a diagonal principal). A figura G.5 mostra essa propriedade. Matrizes Semidiagonais têm vantagens significantes uma vez que os elementos zeros fora da banda não são armazenados e nem operados pelos computadores. A matriz rigidez obtida na análise de elementos finitos, em geral, a estreita banda, conduz ela mesma a uma formulação e solução mais econômica. 96








             = − −−−−− −−−−−− nnnn nnnnnn nnnnnn kk kkk kkk kkk kkk kk K 1, ,11,12,1 1,22,23,2 343332 232221 1211 0...0 0... 00.. ....... ..00 ...0 0...0 Figura G.5 Estrutura de Banda de K. Definição: Uma matriz A n x n é dita ser definida positiva se i. cTAc ≥ 0 para todo vetor c de ordem n; e ii. cTAc = 0 implica em c=0. Observa-se: 1. Uma matriz definida positiva e simétrica possui uma única inversa. 2. Os autovalores de uma matriz definida positiva são reais e positivos. Teorema: A matriz K n x n definida em (G.45) é definida positiva. Prova: i. Seja cA, A= 1, 2, ..., n, e os componentes de c (isso é, c={cA}), um vetor arbitrário. Usa-se estes cA’s para construir um membro de υh,
== n A AA h Ncw 1 , onde os NA’s são as funções bases para υh. Então cTKc = )K de (Definição ),( AB 1,1,

== = n BA BBAA n BA BABA cNNaccKc = a(.,.)) de dade(Bilineari , 11
    

== n B BB n A AA NcNca = ) wde (definição ),( hhh wwa = (2.19))(Por 0 )( 1 0 0 2 , ≥
≥ dxwhx
 ii. Assume cTKc=0. Por parte da prova de (i),
= 1 0 2 , 0)( dxw h x e conseqüentemente wh deve ser constante. Uma vez que wh ∈ υh, wh(1)=0. Combinando estes fatos, conclui-se que wh(x)=0 para todo x ∈ [0, 1], o qual é possível somente se cada cA = 0, A = 1,2,...,n. Então c=0. Note que à parte (ii) depende da definição de K e das condições de contorno essenciais zeros construída dentro da definição de υh. 99 Note também que (G.66) pode ser escrito em termos de (G.67):
= = 2 1 .)()( a e aa e xNx ξξ (G.68) Essa tem a mesma forma que as funções de interpolação. Por referência, nota-se o seguinte resultado: 2 )1( 2, a a aN −== ξ ξ (G.69) 2, e e hx =ξ (G.70) onde he = x2 e – x1 e e e ee x h x 2 )( 1,, == − ξξ (G.71) A descrição local e global dos e-simos elementos são esboçados na Figura G.6. Figura G.6 Descrição local e global do ézimo elemento. G.12. A Matriz de Rigidez e o Vetor Independente de um Elemento Genérico. Para desenvolver mais o ponto de vista dos elementos, assumi-se que o modelo em questão consiste de nel elementos, numerados como na figura G.7. Claro que para este caso nel = n. E seja e a variável índice para os elementos; então 1 ≤ e ≤ nel. 100 Número dos Elementos (e) 1 2 3 nel x1 x2 x3 x4 xn xn+1 Coordenadas dos “nós” Figura G.7 Os Elementos e seus “Nós”. Agora relembrando da definição (global) das matrizes de rigidez e do vetor força [ ] { }
1 , ×× == n A nn AB FFKK
 (G.72) onde
== 1 0 ,,),( dxNNNNaK xBxABAAB (G.73)

++ −+=−+= 1 0 ,1, 1 0 111 ),(),( dxgNNhfdxNgNNahfNF xnxAAAnAAAA δδ (G.74) Em (G.73) assumiu-se NA(x1)=δA1 , como para o espaço de elementos finitos lineares “Piecewise”. As integrais sobre [0,1] podem ser escritas como a soma de integrais sobre os elementos do domínio. Então [ ]eABe n e e KKKK el ==
= , 1 (G.75) { }eAe n e e FFF el ==
= F , 1 (G.76) onde
Ω == e dxNNNNaK xBxA e BA e AB ,,),( (G.77)

Ω + Ω + −+=−+= ee dxgNNhfdxNgNNahfNF xnxAAA e nAAe e A e A ,1,1111 ),(),( δδδ (G.78) e Ωe = [x1e, x2e], o domínio do ézimo elemento. Uma observação importante há fazer é que K e F podem ser construídos pela soma das contribuições das matrizes e vetores dos elementos, respectivamente. Pela definição de NA’s, tem-se que KAB e=0, se A ≠ e ou e+1 ou B ≠ e ou e+1 (G.79) e FA e=0, se A ≠ e ou e+1 (G.80) 101 A situação para um elemento típico, e, é mostrado na figura G.8. Na prática, não se deve é claro, adicionar os zeros mas meramente adicionar os termos diferentes de zeros na locação apropriada. Para esse propósito é muito útil definir a matriz de rigidez para o ézimo elemento ke e o vetor forca elemento fe como se segue: [ ]
{ }
1222 , ×× == ea ee ab e ffkk (G.81)
Ω == e dxNNNNak xbxa e ba e ab ,,),( (G.82)

  =− −= = +=
Ω el e a el a e a negk ne eh Nfdxf e 2 1 1,...,3,20 1δ (G.83)
1 1)(e Linha (e) Linha 1 Coluna ××







 







  = →+← →←          ↓↓ + = n e nn e X X F XX XX ee Coluna K











 Figura G.8 X’s indica termos diferentes de zero; todos os outros termos são zeros. Aqui ke e fe são definidos com respeito à ordenação local, onde Ke e Fe são definidos com respeito à ordenação global. Para determinar onde os componentes ke e fe irão ficar em K e F, respectivamente, requer informações adicionais. Isso é discutido no item seguinte. G.13. A Montagem da Matriz Global de Rigidez e do Vetor Independente. Em um programa computacional de elementos finitos, existe a tarefa de uma sub- rotina de elementos finitos para produzir ke e fe, e=1,2,...,nel, dos dados fornecidos e para prover a sub-rotina de montagem necessita-se de informação adicional tal que os termos em ke e fe possam ser adicionados em locações apropriadas em K e F, respectivamente. Estas informações para montagem são armazenadas em uma matriz chamada LM, a matriz de localizações. Constroe-se uma matriz LM para o problema acima considerado. As dimensões de LM são nel, o número de nós dos elementos, pelo número de elementos; no caso presente, os números são 2 e nel, respectivamente. Dando um número particular de graus de liberdade e um número de elementos (diga-se a e e, respectivamente), o valor retornado pela matriz LM é correspondente ao número global das equações, A, isso é 104 Diversos métodos numéricos são disponíveis na resolução de um sistema linear. Destaque pode ser dados as duas classes abaixo, devido à facilidade de implementação computacional e a sua simplicidade matemática: Métodos Diretos (Baseados no Escalonamento de Matrizes): - Método de Jordan; - Método de Gauss; - Método da Pivotação Parcial e - Método da Pivotação Completa. Métodos Iterativos: - Método de Jacobi; - Método de Gauss-Seidel e - Método SOR (Sucessive Over Relaxation). Nota-se que para os métodos iterativos tem a restrição de seu critério de convergência :
≠ = = n ji i ijii aa 1 |||| . Em certos casos extremos pode-se até usar análise de sistemas mal condicionados e refinamento de sistemas lineares. Além disso, alguns métodos só podem ser aplicados a sistemas especiais como por exemplo o método de Cholesky no qual a matriz A deve ser simétrica. 105 Bibliografia [1] HRENNIK OFF. A. Solution of problems in elasticity by frame work method. Journal of Appied Mechanics. Vol. 8, Nº 4, pp. 169-175, Dec, 1941. [2] MCHENRY, D. A lattice analogy for the s olution of plane stress Problems. Journal of Institution of Civil Engineers. Vol. 21, pp. 59-82, Dec, 1943. [3] COURANT, R., Variational methods for the solution of problems of Equilibrium and Vibrations. Bulletin of the American Mathematical Society. Vol. 49, pp. 1-23, 1943. [4] LEVY, S., Computation of influence coeficients for aircraft structures with discontinuities and sweepback. Journal of Aeronautical Sciences. Vol. 14, Nº 10, pp. 547-560, Oct, 1947. [5] LEVY, S. Structural analysis and influence coefficients for delta wings. Journal of Aeronautical Sciences, Vol 20, Nº 7, pp. 449-454, July, 1953. [6] ARGYRIS, J.H., Energy Theorems and structural analysis. Aircraft Engineering, Oct., Nov., Dec, 1954; and Feb, Mar, Apr, May, 1955. [7] ARGYRIS, J. H. & KELSEY, S. Energy theorems and structural Analysis, Butterworths. London, 1960 (collection of papers published in Aircraft Engineering in 1954 and 1955). [8] TURNER, M. J. & CLOUGH, R. W. & MARTIN, H. C. & TOPP, L. J. Stiffness and deflection analysis of complex structures. Journal of Aeronautical Sciences. Vol. 23, No. 9. Pp.805-824, Sept, 1956. [9] CLOUGH, R. W., The finite Element Method in Plane Stress Analysis. Proceedings. American Society of civil Engineers, Journal . 2nd Conference on Eletronic Computation, Pittsburgh, PA, pp. 345-378, Sept, 1960. [10] MELOSH, R. J., A stiffness matrix for the analysis of thin plates in bending. Journal of the Aerospace Sciences. Vol. 28, Nº 1, pp. 34- 42, Jan. 1961. [11] GRAFTON, P. E., & STROME, D. R. Analysis of axisymmetric shells by the direct stiffness methods. Journal of American Institute of Aeronautics and Astronautics. Vol.1, Nº 10, pp. 2342-2347, 1963. [12] MARTIN, H. C. Plane elasticity problems and the direct stiffness Methods. The trend in Engineering. Vol. 13, pp. 5-19, Jan, 1961. [13] GALLAGHER, R. H. & PADLOG, J. & BIJLAARD, P. P. Stress Analysis of heated complex shapes. Journal of the American Rocket Society. Vol. 32, pp. 700-707, May, 1962. [14] MELOSH, R. J. Structural analysis of solids. Journal of Structural Division. Proceedings of the American Society of Civil Enginners, pp. 205-223, Aug. 1963. [15] ARGYRIS, J. H., Recent advances in matrix methods of structural analysis. Progress in Aeronautical Science, Vol. 4, Pergamon Press, New York, 1964. [16] CLOUGH, R. W. & RASHID, Y. Finite element of axisymmetric Solids. Journal of the Enginnering Mechanics Division. Proceedings of the American Society of Civil Engineers, Vol. 91, pp. 71-85, Feb, 1965. 106 [17] WILSON, E. L. Structural analysis of axisymmetric solids. Journal of the American Institute of Aeronautics and Astronautics, Vol. 3, Nº 12, pp. 2269-2274, Dec. 1965. [18] TURNER, M. J. & DILL, E. H. & MARTIN, H. C. & MELOSH, R. J. Large deflections of structures subjected to heating and external Loads. Journal of Aeronautical Sciences. Vol. 27, Nº 2, pp. 97-107, Feb 1960. [19] GALLAGHER, R. H., AND PADLOG, J., Discrete element approach to structural stability analysis. Journal of the American Institute of Aeronautics and Astronautics. Vol. 1, Nº 6, pp. 1437-1439, 1963. [20] ZIENKIEWICZ, O. C. & WATSON, M. & KING I. P. A Numerical methods of visco-elastic stress analysis. International Journal of Mechanical Sciences, Vol. 10, pp. 807-827, 1968. [21] ARCHER, J. S. Consistent matrix formulations for structural analysis using finite-element techniques. Journal of the American Institute of Aeronautics and Astronautics. Vol. 3, Nº 10, pp. 1910-1918, 1965. [22] ZIENKIEWICZ, O. C. & CHEUNG, Y. K. Finite elements in the solution of field problems. The Enginner. pp. 507-510, Sept. 24, 1965. [23] MARTIN, H.C. Finite element analysis of fluids flows. Proceedings of the Second Conference on Matrix Methods in Structural Mechanics, Wright-Patterson Air Force Base, Ohio, pp. 517-535, Oct, 1968. (AFFDL-TR-68-150, Dec. 1968; AD-703-685, N.T.I.S.) [24] WILSON, E. L. & NICKEL, R. E. Application of the finite element methods to heat conduction analysis. Nuclear Engineering and Design. Vol. 4, pp 276-286, 1966. [25] SZABO, B. A. & LEE, G. C. Derivation of stiffness matrices for problems in plane elasticity by Galerkin’s Methods. International Journal of Numerical Methods in Enginneering. Vol. 1, pp. 301-310, 1969. [26] ZIENKIEWICZ, O. C., & PAREKH. C. J., Transient field problems: two-dimensional and three-dimensional analysis by isoparametric finite elements. International Journal of Numerical Methods in Enginnering. Vol. 2, Nº 1, pp. 61-71, 1970. [27] LYNESS, J. F. & OWEN, D. R. J. & ZIENKIEWICZ, O. C. Three- dimensional magnetic field determination using a scalar potential. A finite element solution. Transactions on Magnetics. Institute of Eletrical and Electronics Engineers, pp. 1649-1656, 1977. [28] BELYTSCHKO, T. A Survey of numerical methods and computer programs for dynamic structural analysis. Nuclear Engineering and Design. Vol. 37, Nº 1, pp. 23-34, 1976. [29] BELYTSCHKO, T. Efficient Large-Scale Nonlinear Transient Analysis by Finite Elements. International Journal of Numerical Methods in Enginnering. Vol. 10, Nº 3, pp. 579-596, 1976. [30] HUISKIES, R. & CHAO, E. Y. S. A survey of finite elements analysis in orthopedic biomechanics: The First Decade. Journal of Biomechanics. Vol. 16, Nº 6, pp. 385-409, 1983. [31] JOURNAL OF BIOMECHANICAL ENGINNERING. Transactions of american society of mechanical engineers. ( published quarterly ) (1st issue published 1977).