Baixe Elementos Finitos - Apostilas - Análise de Sistemas Parte1 e outras Notas de estudo em PDF para Informática, somente na Docsity! 1 Capítulo III 3.2. Método das Diferenças Finitas Este item 3.2 foi extraído da referência [55] Carnahan 1969, sofrendo modificações. Em todo o desenvolvimento da análise numérica, utiliza-se as diferenças finitas, com um roteiro a ser seguido numa seqüência lógica. Partindo da definição de operadores numéricos de diferenças finitas (∆, ∇, µδ e δ), introduz-se o conceito de interpolação através das fórmulas de Gregory-Newton e Stirling, que utilizam estes operadores. Em seguida, introduz-se a derivação numérica e a integração numérica (Quadratura) por meio da derivação e integração da fórmula de Gregory-Newton, chegando às fórmulas de Newton-Cotes. Seguindo essa técnica, praticamente todos os tópicos da análise numérica podem ser introduzidos por meio das diferenças finitas, e depois, eles são desenvolvidos para além das diferenças finitas. No estudo de análise numérica de equações diferenciais não é diferente. Sugere-se uma introdução por meio dos métodos das diferenças finitas (do inglês: Finite Difference Methods ou FDM), e posteriormente, o desenvolvimento do assunto para além das diferenças finitas, como por exemplo, a introdução do Método dos Elementos Finitos. Com isso, segue-se um desenvolvimento didático muito próximo do desenvolvimento histórico, uma vez que, conforme foi citado no item 1.2 da Introdução, o método das diferenças finitas surgiu antes do Método dos Elementos Finitos. O método das diferenças finitas pode ser utilizado para resolver problemas de valor de contorno ou valor inicial, envolvendo equações diferenciais ordinárias ou parciais. Assim, este método pode ser usado para solucionar as equações de modelos a parâmetros concentrados ou distribuídos. A técnica consiste em substituir cada derivada ou diferencial das equações diferenciais por aproximação de diferenças finitas ou acréscimo finitos das variáveis, como mostra as equação 3.1 abaixo: ∆∆ ∆≈ ∂∂ ∂ ∆ ∆≈ ∂ ∂ ∆ ∆≈ ∂ ∂ ∆ ∆≈ ∂ ∂ ∆ ∆≈ ∆ ∆≈ ∆ ∆≈ ∆≈∆≈ yx u yx u x u x u x u x u x u x u x y dx yd x y dx yd x y dx dy ydyxdx . ,,, ,, , 22 3 3 3 3 2 2 2 2 3 3 3 3 2 2 2 2 (3.1) O Método de Elementos Finitos é bem mais recente que o anterior, sendo mais genérico, e podendo ser aplicado a complexas estruturas geométricas e a ambientes com várias mudanças de meio. Ele possui uma formulação matemática mais trabalhada, sendo, portanto, um conjunto de técnicas e métodos que se baseia na discretização do problema em elementos pequenos e na aproximação de cada elemento por um conjunto de polinômios. Considere, primeiramente, o problema formado por equações diferenciais ordinárias (EDOs). Existem dois tipos. Um deles, é o problema de valor inicial, que 2 assume a forma geral abaixo: F[ t, y(t), y´(t) ] =0, t>0 t=t0, y=yo, y´=y´0 (3.2) Onde, o t é a variável independente, usualmente o tempo; y é um vetor de variáveis dependentes; y´ é a sua derivada em relação a t; F é um vetor de funções de t, y, e y´; e, finalmente, yo e y´o são vetores que representam as condições iniciais do problema. Note- se, que o domínio da variável t é semi-infinito, e que a solução deste problema deverá ser obtida marchando-se no tempo, a partir da condição inicial. Caso exista, pelo menos, uma função dentro do vetor F que não dependa de nenhum elemento do vetor y´, a equação representa um sistema de equações algébrico-diferenciais (sistema de EAD). [68] O outro tipo de problema é o do valor de contorno, que assume a seguinte forma geral para sistemas de segunda ordem: F[ x, y(x), y´(x), y´´(x) ] = 0, xo<x<xf x=xo, go(x,y,y´)=0 (3.3) x=xf, gf(x,y,y´)=0 Onde, o x é a variável independente, usualmente, uma coordenada espacial; y é o vetor de variáveis dependentes; y´ e y´´ são as suas derivadas: primeira e segunda, respectivamente, em relação a x; F é um vetor de funções; e go e gf são vetores de funções que representam as condições de contorno nos limites do domínio do sistema de equações. O objetivo do Método das Diferenças Finitas é transformar um problema composto por equações diferenciais em um problema formado por equações algébricas. O primeiro passo, nesta direção, é a chamada discretização do domínio da variável independente. A discretização consiste em dividir o domínio de cálculo em um certo número de subdomínios. Para um domínio semi-infinito, existem infinitos subdomínios. Quando o domínio é finito, o número de subdomínios também o é, e digamos que seja J. Em qualquer caso, estipulam-se os pontos que delimitam os subdomínios, que, no caso de um domínio finito, são iguais a (J+1), em número. Note-se, que os subdomínios podem ter o mesmo tamanho, gerando uma malha uniforme, ou então, formando uma malha não-uniforme. Embora as discretizações baseadas no primeiro tipo de malha sejam mais simples, existem vantagens numéricas, em muitos casos, no uso de malhas não-uniformes. O segundo passo é gerar aproximações para as derivadas das variáveis dependentes que aparecem nas diferenciais, nos pontos discretos, xj ou tj, isto é, obter y´j e y´´j, utilizando apenas os valores de y nestes pontos discretos: yj. Finalmente, aplicam- se as equações diferenciais ordinárias aos pontos discretos xj, substituindo as aproximações obtidas para y´j e y´´j. Isto gera sistemas de equações algébricas na forma: f( yj ) =0, (3.4) Onde, o f é um vetor de equações algébricas que depende dos valores desconhecidos yj, sendo que esta dependência varia conforme o tipo de problema, de contorno ou inicial. Este sistema de equações, quer seja ele linear ou não linear, pode ter 5 ou (2+h) yj+1 - 4(h 2+1) yj + (2-h) yj-1 =0, j=1,2,3,4,.....,J-1 (3.15) e yo=0 e yJ=1 Logo, caímos em um sistema linear de equações algébricas. As condições de contorno do problema dado acima são chamadas de primeiro tipo, isto é, definem o valor da variável no contorno, sendo facilmente incorporadas ao sistema algébrico das equações discretizadas. Diversos outros tipos de condições de contorno são possíveis, sendo que, a sua utilização no sistema de equações algébricas discretizadas, torna-se um pouco mais elaborada. Em geral, a condição de contorno pode ser não-linear, na qual go e gt são funções arbitrárias de x, y e y´. Entretanto, são três, os tipos existentes de condições de contorno lineares. Diz-se, que a condição de contorno é do primeiro tipo, quando o valor da variável dependente é dado no contorno, sendo facilmente utilizada nas equações discretizadas. O problema acima serve de exemplo, e a forma geral é dada por: x=xc, y=yc (3.16) Quando a condição de contorno é de segundo tipo, o valor da derivada da variável dependente é dado no contorno, isto é: x=xc, y´=y´c (3.17) Esta condição de contorno tem que ser discretizada, para ser combinada com o sistema algébrico discretizado, fazendo: h yy y jjj 2 11 −+ −=′ (3.18) A condição de contorno é dita de terceiro tipo, quando tem a seguinte forma geral: x=xc, ay´+by=c (3.19) Nela, a, b e c são constantes conhecidas. O seu tratamento é similar ao dado às condições de contorno de segundo tipo. Considere, agora, o seguinte problema de valor inicial que envolve apenas uma diferencial ordinária na sua forma normal: y´=f(t,y) com t=0, y(0)=yo (3.20) Sendo o intervalo genérico entre tj e tj+1, considere diferentes formas de aproximar a derivada primeira, utilizando, como informação conhecida, apenas o ponto j, isto é y(tj)=yj . Utilizando a aproximação de diferenças finitas para frente, para y´j, obtém-se : 6 Yj+1=yj+hf(t,yj) + O(h 2), h=tj+1-tj (3.21) Fórmula que permite calcular yj+1 a partir de yj , com erro da ordem de h 2. Esta equação é explícita no valor desconhecido de yj+1, sendo, pois, o método denominado de explícito. Especificamente, representa o método explícito de Euler. Não caímos em um sistema linear, tendo a solução direta. Caso, por outro lado, resolvermos utilizar a aproximação de diferenças finitas para trás de y´j+1 , dada com j+1 no lugar de j, podemos aplicar no ponto tj+1 e escrever: yj+1=yj + h f(tj+1, yj+1) + O(h 2), h=tj+1-tj (3.22) Fórmula que calcula yj+1 a partir de yj, com erro da ordem de h 2. Note-se que, no caso geral, a equação é não linear no valor desconhecido de yj+1, sendo, pois, necessário, utilizar-se um método adequado à solução de problemas não lineares, para se obter o valor de yj+1.Assim, como yj+1 não pode ser explicitado a partir da equação, o método é denominado de implícito. Mais especificamente, este é chamado método implícito de Euler. Além dos dois métodos acima, podemos, ainda, obter um terceiro, a partir da aproximação por diferença central de y´j+1/2, no intervalo considerado. Aplicando ao meio do intervalo, temos: yj+1=yj + h f(tj+h/2, yj+1/2) + O(h 3), h=tj+1-tj (3.23) Fórmula que ainda não pode ser usada para obter yj+1, porque o valor de f no ponto considerado, não é conhecido. Entretanto, expandindo fj e fj+1 em série de Taylor, em torno do ponto tj+1/2=tj+h/2, tem-se que: yj+1=yj+h/2 . [ f(tj,yj) + f(tj+1,yj+1) ] + O(h 3), h= tj+1-tj (3.24) Fórmula que permite calcular yj+1, ainda que de forma implícita. Este método é ainda implícito, sendo denominado de método trapezoidal (ou de Crank-Nicholson). Considere a solução numérica do problema de valor inicial, abaixo: y´=-y2, t>0 com t=0, y=1 (3.25) Ela utiliza os métodos de Euler, até o ponto t=1, com passos uniformes de integração de 0,2. A solução analítica deste problema é: y(t)=1/(1+t) (3.26) Aplicando os métodos, temos: 7 Tabela 3.1 Comparação dos métodos de diferenças finitas aplicados à equação 3.25. [68]. T Euler Explícito Euler Implícito Trapezoidal Solução Analítica 0,0 1,0000 1,0000 1,0000 1,0000 0,2 0,8000 0,8541 0,8310 0,8333 0,4 0,6720 0,7434 0,7113 0,7143 0,6 0,5817 0,6572 0,6220 0,6250 0,8 0,5140 0,5880 0,5528 0,5556 1,0 0,4612 0,5316 0,4975 0,5000 No caso de problema de valor inicial com equações de ordem maior do que 1, uma redução de ordem deve ser feita para uma ordem menor, com a inserção de uma variável no sistema. Exemplo: y´´=3y´+y+5x com y(0)=1 e y´(0)=2 (3.27) Introduz-se a variável z com z=y´, derivando z´=y´´. Assim, temos agora, dois P.V.I. (Problema de Valor Inicial) de primeira ordem: = =′ = ++=′ 1)0( 2)0( 53 y zy z xyzz (3.28) Resolvemos, simultaneamente, os dois P.V.I., construindo uma tabela de x, y e z. [52] Os métodos de Runge-Kutta são de ponto simples, explícitos, mas com diversos estágios, de modo a se obter uma maior ordem de aproximação. A idéia básica deste tipo de método é definir que a variação da variável dependente, no passo em questão, é dada por uma média ponderada de variações desta variável, calculadas com avaliações diferentes da função derivada, isto é: >++∆=∆ ∆=∆ ∆=− = + 1 com , ),(. ),(.1 1 1 ibyatfty ytfty yCyy ijiji jj n i iijj (3.29) Nela, Ci, ai e bi são coeficientes a serem determinados. Usando a equação acima, e ou igualando-a à série de Taylor, determinamos estes coeficientes para a diversa ordem de Runge-Kutta. Mostra-se, que Runge-Kutta é, na verdade, uma variação do Método das Diferenças Finitas. Por exemplo: Para segunda ordem: C1+C2=1, aC2=1/2 e bC2=1/2. Se escolhermos C2=1/2, chegamos a Euler Modificado, onde C1=C2=0,5 e a=b=1. 10 2 ,1,,1 2 ,1,,1 2 , 2 )( 22 x uuu h uuu x u jijijijijijiji ∆ +− = +− = ∂ ∂ −+−+ (3.37) A última etapa para a solução de uma equação diferencial parcial por diferenças finitas, é substituir as aproximações das derivadas na equação e nas suas condições de contorno, gerando um sistema algébrico, cuja solução fornece a solução aproximada do problema original. Dado o exemplo abaixo, considere uma barra fina metálica de 1 metro de comprimento, e estude a condução de calor, segundo a equação: 2 2 x u t u ∂ ∂= ∂ ∂ (3.38) u(x,t) é a temperatura da barra na posição x e instante t. Condição Inicial: u(x,0)=0oC Condições de Contorno: u(0,t)=u(1,t)=100oC Passos: ∆x=0,2metros e ∆t=0,01segundos Substituindo, na equação diferencial, as aproximações de diferenças finitas das equações 3.34 e 3.37, e trocando o índice i por j para a derivada em relação a t, tem-se: 2 ,1,,1,1, )( 2 x uuu t uu jijijijiji ∆ +− = ∆ − −++ (3.39) Isolando o termo ui,j+1 , que corresponde ao tempo seguinte, e substituindo os valores de ∆x e ∆t, tem-se : 4 2 ,1,,1 1, jijiji ji uuu u −++ ++ = (3.40) A equação acima fornece a temperatura do tempo j+1 em função do tempo anterior j. Quando se tem uma equação deste tipo, onde o novo valor é calculado em função de valores conhecidos, dá-se o nome de método explícito, pois pode-se explicitar um valor desconhecido em função de outros já conhecidos. Em seguida, pode-se chegar na tabela abaixo: Tabela 3.2: Solução numérica do problema de condução de calor. [55] I 0 1 2 3 4 5 J t \ x 0,0 m 0,2 m 0,4 m 0,6 m 0,8 m 1,0 m 0 0,00 0 0 0 0 0 0 1 0,01 100 0 0 0 0 100 2 0,02 100 25 0 0 25 100 3 0,03 100 37,5 6,25 6,25 37,5 100 4 0,04 100 45,3 14,0 14,0 45,3 100 5 0,05 100 51,1 21,8 21,8 51,1 100 6 0,06 100 56,0 29,1 29,1 56,0 100 11 Veja-se, agora, um outro exemplo: considere uma placa quadrada fina metálica de 1 metro de comprimento, e encontre a temperatura de equilíbrio (após um tempo muito grande), segundo a equação de Laplace: 0 2 2 2 2 = ∂ ∂+ ∂ ∂ y u x u (3.41) u(x,y) é a temperatura da barra na posição x e y como na figura abaixo: y 1 x 0 1 Figura 3.2 Placa Quadrada do Exemplo Condições de Contorno: u(0,y)=u(x,0)=100oC e u(1,y)=u(x,1)=0oC Passos: ∆x=∆y=0,25 metros Substituindo, na equação diferencial, as aproximações de diferenças finitas da equação 3.37, e trocando o índice i por j para a derivada em relação a y, tem-se: 0 )( 2 )( 2 2 1,,1, 2 ,1,,1 = ∆ +− + ∆ +− −+−+ y uuu x uuu jijijijijiji (3.42) Lembrando que ∆x = ∆y, corta os denominadores e isolando o termo ui,j , tem-se : 4 1,1,,1,1 , −+−+ +++= jijijijiji uuuu u (3.43) A equação acima fornece a temperatura na posição i, j , em função da média aritmética das temperatura de cima, de baixo, da direita e da esquerda. Quando se tem uma equação deste tipo, onde o novo valor é calculado em função de valores desconhecidos, dá-se o nome de método implícito. Em seguida, pode-se chegar a um sistema linear, onde a equação 3.43 é expandida para cada ponto do interior da placa, seguindo a numeração dada na figura abaixo: y 100 0 0 0 0 100 u7 u8 u9 0 100 u4 u5 u6 0 100 u1 u2 u3 0 100 100 100 100 100 0 0,25 0,5 0,75 1 x Figura 3.3 : Malha da discretização da placa. Expandindo a equação 3.43 para cada ponto do interior, de u1 a u9, chega-se às equações: 12 =+−− =−+−− =−+− =−+−− =−−+−− =−−+− =−+− =−−+− =−− +++ = +++ = +++ = +++ = +++ = +++ = +++ = +++ = +++ = 04 04 1004 04 04 1004 1004 1004 2004 4 00 4 0 4 0100 4 0 4 4 100 4 1000 4 100 4 100100 986 9875 874 9653 86542 7541 632 5321 421 68 9 579 8 48 7 395 6 2846 5 175 4 62 3 513 2 42 1 uuu uuuu uuu uuuu uuuuu uuuu uuu uuuu uuu uu u uuu u uu u uuu u uuuu u uuu u uu u uuu u uu u (3.44) Resolvendo o sistema linear da equação 3.44, chega-se à solução abaixo: 000,50428,71714,85 571,28000,50428,71 285,14571,28000,50 321 654 987 === === === uuu uuu uuu (3.45) Observa-se a simetria em relação à diagonal principal da placa, e constata-se que os valores da temperatura, nos vértices da placa, não entraram no cálculo. [55] 15 Figura 3.5 - Exemplos de elementos empregados em (a) uma, (b) duas, e (c) três dimensões. 3.3.1.2. Equações dos Elementos (Processamento) A seguir, desenvolvem-se equações, a fim de aproximar a solução de cada elemento. Isso envolve dois sub-passos. Primeiro, escolhe-se uma função apropriada com coeficientes desconhecidos, que serão usados para aproximar a solução. Por último, avaliam-se os coeficientes, em que as funções se aproximam da solução, de forma considerada ótima. [54] Escolha das Funções de Aproximação Considera-se que, por serem de fácil manipulação matemática, os polinômios são freqüentemente empregados para este propósito. Para o caso unidimensional, a Elemento Linear (a) Unidimensional Nó Linha Nodal Elemento Elemento Quadrílateral Triangular (b) Bidimensional Elemento Hexaédrico Plano Nodal Elemento Tetraédrico (c) Tri-dimensional 16 alternativa mais simples é um polinômio de primeira ordem, ou uma linha reta: u(x) = a0 + a1 x (3.46) Nesta fórmula, u(x) é a variável dependente; a0 e a1 são constantes; e x é a variável independente. Essa função deve passar através dos valores u(x) nos pontos finais dos elementos em x1 e x2. Portanto: u1 = a0 + a1 x1 u2 = a0 + a1 x2 Onde u1 = u (x1) e u2 = u (x2). Estas equações podem ser resolvidas, usando a regra de Cramer, onde: a0 = (u1 x2 u2 x1) / (x2 x1) a1 = (u2 u1) / (x2 x1) Este resultado pode, então, ser substituído na Eq. (3.46), a qual, depois de se arrumar os termos, pode ser escrita como: u = N1 u1 + N2 u2 (3.47) onde N1 = (x2 x) / (x2 x1) (3.48) e N2 = (x x1) / (x2 x1) (3.49) A equação (3.47) é chamada função de aproximação ou de forma, e N1 e N2 são chamados de funções de interpolação. Inspecionando melhor, percebe-se que a Eq. (3.47) é, de fato, o polinômio interpolador de primeira ordem de Lagrange. Ela fornece um significado, para predizer valores intermediários (que é interpolar) entre valores dados u1 e u2 nos nós. 17 Figura 3.6 (b) Uma função de aproximação ou forma para (a) um elemento linear. As correspondentes funções de interpolação são mostradas em (c) e (d). Note-se, que a soma das funções de interpolação (N1+N2) são iguais a 1. Em adição, lidar com equações lineares facilita operações como a diferenciação e integração. Tais manipulações, mais à frente, serão importantes em outros itens. A derivação da Eq. (3.47) é: 2 2 1 1 u dx dN u dx dN dx du += (3.50) De acordo com as Eq. (3.48) e (3.49), as derivadas de N1 e N2 podem ser calculadas como: 12 2 12 1 1 1 xxdx dN xxdx dN − = − −= (3.51) E, portanto, a derivada de u é: )( 1 21 12 uu xxdx du +− − = (3.52) Em outras palavras, essa é a diferença dividida, representando a inclinação da reta conectada nos nós. A integral pode ser expressa como: Nó 1 Nó 2 (a) u1 u u2 x1 x2 (b) 1 N1 x1 x2 (c) N2 1 x1 x2 (d) 20 u(x,y) é a temperatura da barra na posição x e y como na figura abaixo: y 1 x 0 1 Figura 3.7 Placa Quadrada do Exemplo. Condições de Contorno: u(x,0)=0oC, u(x,1)=100oC, u(0,y)=75oC e u(1,y)=50oC 3.3.3. Problema Bidimensional Este item foi extraído da referência [54]. Embora o número de operações matemáticas cresça bastante, a extensão da aproximação de elementos finitos para duas dimensões é, conceitualmente, similar à aplicação unidimensional discutida. Assim, ela segue os mesmos passos, como mostrado no subitem 3.2.1. 3.3.3.1. Discretização Uma variedade de simples elementos, como triângulos e quadriláteros, são usualmente empregados para fracionar os elementos finitos em duas dimensões. Na discussão presente, limitar-se-á a elementos triangulares do tipo descrito na Fig. 3.8. y 3 2 1 x Figura 3.8 Um elemento triangular 3.3.3.2. Equações dos Elementos Como no caso unidimensional, o próximo passo é desenvolver uma equação para aproximar a solução ao elemento. Em um elemento triangular, a aproximação mais simples é um polinômio linear [Compare com a Eq. (3.46)]. u(x, y) = a0 + a1,1x + a1,2y (3.58) Onde u é a variável dependente, os as representam os coeficientes, e x e y são variáveis independentes. Essa função deve passar pelos valores de u(x, y) nos nós do triângulo (x1, y1), (x2, y2), e (x3, y3). Portanto: u1(x, y) = a0 + a1,1x1 + a1,2 y1 u2(x, y) = a0 + a1,1x2 + a1,2 y2 21 u3(x, y) = a0 + a1,1x3 + a1,2 y3 Ou na forma matricial: = 3 2 1 2,1 1,1 0 33 22 11 1 1 1 u u u a a a yx yx yx Este pode ser resolvida a: )]()()([ 2 1 1221331132233210 yxyxuyxyxuyxyxuA a e −+−+−= (3.59) )]()()([ 2 1 2131323211,1 yyuyyuyyuA a e −+−+−= (3.60) )]()()([ 2 1 1233122312,1 xxuxxuxxuA a e −+−+−= (3.61) Onde Ae é a área do elemento triangular: )]()()[( 2 1 122131132332 yxyxyxyxyxyxAe −+−+−= As equações (3.59) até (3.61) podem ser substituídas na Eq. (3.58). Depois de agrupar os termos, o resultado pode ser expresso como: u = N1u1 + N2u2 + N3u3 (3.62) Onde: ])()()[( 2 1 233223321 yxxxyyyxyxA N e −+−+−= ])()()[( 2 1 311331132 yxxxyyyxyxA N e −+−+−= ])()()[( 2 1 122112213 yxxxyyyxyxA N e −+−+−= A equação (3.62) fornece uma maneira de predizer valores intermediários para o elemento, com base nos valores dos nós. A Figura 3.9 mostra a função de aproximação com as correspondentes funções de interpolação. Percebe-se, que a soma das funções de interpolação é sempre igual a um. 22 Figura 3.9 (a) Uma linear função de aproximação para um elemento triangular. As correspondentes funções de interpolação são mostradas em (b) até (d). Fonte [54] Também no caso unidimensional, vários métodos são disponíveis para desenvolver as equações dos elementos, baseadas na subjacente PDE e nas funções de aproximação. As equações resultantes são, consideravelmente, mais complicadas que as do caso unidimensional. Contudo, devido às funções de aproximação serem usualmente polinômios de mais baixa ordem, como na Eq. (3.58), os termos da matriz final do elemento consistirá a polinômios de baixa ordem e constantes. 25 Figura 3.10 - Domínio a dois elementos triangulares. A escolha de cada elemento das funções de aproximação definirá o seu tipo, e caracterizará, pela seqüência, a natureza das funções de aproximação. Há uma ligação matemática rigorosa entre a escolha da natureza das funções de aproximação (lineares, quadráticas, cúbicas) e a forma dos elementos, sempre definidas por trechos sobre a discretização. Na aplicação deste método, é preciso escolher um conjunto de funções de projeção (ou funções de ponderação) Φ1, Φ2,..., ΦNN antes de escrever as equações de projeção L(u*) sobre cada uma destas funções. Na técnica de Galerkin tem-se N1=Φ1 , N2=Φ2 e assim por diante.[37] 0 1 1 =Ω − Φ = dfNuL NN j jj 0 1 2 =Ω − Φ = dfNuL NN j jj (3.64) ... 0 1 =Ω − Φ = dfNuL NN j jjNN Obtém-se, de novo, um sistema de NN equações algébricas a resolver para se encontrar as NN incógnitas u1, u2,..., uNN. Observações Importantes: A aplicação do método dos elementos finitos conduz à substituição de uma equação ou sistema de equações a derivadas parciais por um sistema de equações algébricos, contento coeficientes da função de aproximação, que são, na realidade, os valores das funções de aproximação nos nós do domínio discretizado. Como no Método dos Resíduos Ponderados as funções de ponderação são idênticas às funções de aproximação, o sistema de equações obtido é idêntico àquele obtido pela formulação variacional. Se o operador L é linear, então a funcional é quadrática, implicando na linearidade do sistema de equações algébricas obtido. Da mesma forma, no Método dos Resíduos Ponderados, a linearidade de L implica na linearidade das equações (3.64), pois, em cada uma delas, os coeficientes ui podem ser colocadas em evidência na integral. 3 1 4 2 26 ( )( ) = Ω = ΩΦ−ΩΦ=Ω −Φ NN i iijj NN j jji fddNLudfNuL 11 3.3.3.4. Condições de Contorno e Montagem A incorporação das condições de contorno e a montagem do sistema matricial também se tornam mais complicados, quando a técnica dos elementos finitos é aplicada em problemas de duas ou três dimensões. Contudo, como na derivação da matriz elemento, as dificuldades relativas ao mecanismo dos processos são maiores do que a complexidade conceitual. Por exemplo, o estabelecimento da topologia do sistema, o qual era trivial para o caso de uma dimensão, torna-se um problema de grande importância em duas ou três dimensões. Em particular, a escolha do esquema de numeração irá ditar a estrutura do sistema matricial resultante e, portanto, a eficiência com a qual ele pode ser resolvido. A Figura 3.11 mostra um esquema desenvolvido para uma placa plana em equilíbrio térmico, solucionada por meio do método dos elementos finitos. Figura 3.11 - Um esquema de numeração dos nós e elementos para uma aproximação por elementos finitos, para uma placa plana em equilíbrio térmico. [54] Logo para o exemplo tem-se: 2 2 2 2 y T x T R ∂ ∂+ ∂ ∂= Integrando pelo Método dos Resíduos ponderados e associado à técnica de Galerkin, tem-se: 0 2 2 2 2 =Ω ∂ ∂+ ∂ ∂ Ω d y T x T Ni Y X 100 oC O oC 50 oC 75 oC 27 onde Ni são as funções de teste no caso as funções de Lagrange N1, N2 e N3 da função de aproximação T = T1 . N1 + T2 . N2 + T3 . N3. Após uma integração por partes (Teorema de Green), obtém-se: dS n T Nd y N y T x N x T S i ii ∂ ∂−=∆ ∂ ∂ ∂ ∂+ ∂ ∂ ∂ ∂ ∆ para cada elemento triangular. O termo do lado direito representa o fluxo sobre o contorno do elemento triangular (Percorrendo no sentido anti-horário). Assim para um elemento genérico 1,2 e 3 pode-se escrever: [ ]( ) [ ]( ) [ ]( )312312 F 31 F 23 F 12 ladonoFluxoladonoFluxoladonoFluxodSn T N S i ++=∂ ∂− Pode-se aplicar as equações acima em cada elemento da malha: Para o Elemento (1) tem-se: Nó 7 (0,25;0,25) Nó 1 (0;0) Nó 2 (0,25;0) Figura 3.12 Elemento 1 da Malha com as coordenadas (x,y) de cada nó. Para i=1: dS n T Ndydx y N y T x N x T S x ∂ ∂−= ∂ ∂ ∂ ∂+ ∂ ∂ ∂ ∂ 1 25,0 0 0 11 Substituindo a integral do lado direito pelo fluxo em cada lado do triangulo e lembrando que N1 é zero para o lado oposto ao vértice 1 (F27=0) tem-se: 7112 25,0 0 0 11 FFdydx y N y T x N x Tx += ∂ ∂ ∂ ∂+ ∂ ∂ ∂ ∂ Resolvendo o lado esquerdo tem-se: 711221 5,05,0 FFTT +=− (3.65) Para i=2 tem-se de maneira similar: 30 3.3.3.5. Solução e Apresentação do Resultado (Pos-Processamento) Embora o mecanismo seja complicado, o sistema matricial é meramente um conjunto de n equações simultâneas, que pode ser solucionada, para se encontrar os valores das variáveis dependente nos n nós. A Figura 3.13 mostra uma solução possível, no caso de uma placa em equilíbrio térmico. Figura 3.13 - A distribuição da temperatura de uma placa em equilíbrio térmico, sendo calculada por meio do Método dos Elementos Finitos. [54] 3.3.3.6. Comparação com a Solução pelo MDF Mesmo neste exemplo com meio homogêneo e geometria simples, nota-se uma vantagem no MEF, pois não se é obrigado a usar malhas retangulares de mesmo tamanho como no MDF. Em outros problemas com meios não homogêneos, geometrias complexas e parâmetros não lineares as vantagens do MEF aumentam bastante. O que torna o MEF um método bem mais amplo e com uma aplicabilidade prática muito maior que o MDF. Assim o aumento na complexidade matemática dos elementos finitos é compensada por sua portabilidade e eficiência na solução de equações diferenciais parciais encontradas na engenharia. Por esta razão o MEF é o método mais comumente implementado pelas ferramentas de CAE e para muitos autores é a base da Engenharia Assistida por Computador. 31 Apêndice A Ajuste Ótimo da Função de Aproximação em uma Dimensão A.1. Solução de Elementos Finitos para Molas em Séries Figura A.1 (a) Uma série de molas interconectadas. Uma extremidade é fixada na parede, enquanto a outra é submetida a uma força constante F. (b) Representação em elementos finitos. Cada mola é representada por um elemento. Portanto, o sistema consiste de quatro elementos e cinco nós. [54] Descrição do problema: a figura A.1 mostra uma série de molas Interconectadas. Uma extremidade é fixada a uma parede, enquanto a outra é sujeita a uma força constante F. Usando, passo a passo, os procedimentos do Método dos Elementos Finitos, determina-se o deslocamento das molas. Solução − Discretização: o modo de particionar esse sistema é, obviamente, tratar cada mola como um elemento. Assim, o sistema consiste de quatro elementos e cinco nós (Fig. A.1b). Equações dos Elementos: como este sistema é muito simples, suas equações dos elementos podem ser escritas diretamente, sem o recurso da aproximação matemática. Este é um exemplo de aproximação direta para os elementos derivados. 32 Figura A.2 - Um diagrama de corpo livre de um sistema de molas. [54] A Figura A.2 mostra um elemento individual. O relacionamento entre a força F e o deslocamento x pode ser representado, matematicamente, pela lei de Hooke: F = k x Onde, k representa a constante da mola, a qual pode ser interpretada como a força requerida para causar uma unidade de deslocamento. Se uma força F1 é aplicada no nó 1, o seguinte balanço de força (reação) deve segurar: F = k (x1 x2) Onde, x1 é deslocamento do nó um da sua posição de equilíbrio; e x2 o deslocamento do nó dois da sua posição de equilíbrio. Assim, x2 x1 representa o quanto a mola é alongada ou comprimida e relativa ao equilíbrio (Fig. A.2). Essa equação pode também ser escrita como: F1 = k x1 k x2 Para um sistema estacionário, um balanço de forças também necessita que F1=F2 e, portanto: F2 = -k x1 + k x2 Estas duas simultâneas equações especificam o comportamento do elemento em resposta às forças aplicadas. Podem ser escritas numa forma matricial, como: = − − 2 1 2 1 F F x x kk kk Ou então: [ k ] { x } = { F } , Onde a matriz [ k ] é a matriz propriedade do elemento. Neste caso, é também referenciada com a matriz rigidez do elemento. Note-se, que esta última equação tem sido 35 Figura A.3 (a) O diagrama do sistema original. (b) O sistema depois da aplicação da força constante. Os deslocamentos são indicados no espaço entre os dois sistemas. [54] A.2. Exemplo de uma Haste Sendo Aquecida A figura A.4 mostra um sistema modelado pela equação de Poisson, na forma unidimensional: )( 2 2 xf dx Td −= (A.1) Onde, f(x) é uma função que define a fonte de calor ao longo da haste, e onde, as pontas da haste, são mantidas numa temperatura fixa T(0, t) = T1 e T(L, t) = T2. Essa não é uma equação diferencial parcial, mas uma equação diferencial ordinária com condições de contorno. Esse modelo simples é usado, porque ele permitirá introduzir a aproximação por elementos finitos, sem algumas das complicações, como por exemplo, um PDE de duas dimensões. Figura A.4 (a) Uma longa e fina haste sujeita a fixas condições de contorno e a uma contínua fonte de calor ao longo do seu eixo. (b) A representação de elementos finitos consistindo de quatro elementos de igual comprimento e cinco nós.[54] Exemplo: solução analítica para uma haste sendo aquecida. Descrição do Problema − Resolver a Eq. (A.1) para uma haste de 10 cm, com as seguintes condições de contorno: T(0, t) = 40 oC e T(10, t) = 200 oC e uma fonte de calor 36 uniforme f(x) = 10. Solução − A equação a ser resolvida é: 10 2 2 −= dx Td Assume-se uma solução na forma: T = a x2 + b x + c Ela é diferenciada duas vezes, resultando T´´ = 2 a. Substituindo este na equação diferencial, temos: a = -5. As condições de contorno são usadas para avaliar os coeficientes restantes. Da primeira condição, para x=0: 40 = -5 (0)2 + b(0) + c ou c = 40. Similarmente, para a segunda condição: 200 = -5(10)2 + b(10) + 40 A qual é solucionada dando b = 66. Portanto, a solução final é: T = -5 x2 + 66 x + 40 O resultado é mostrado na figura A.5 Figura A.5 - Distribuição de temperatura ao longo de uma haste aquecida por uma fonte uniforme de calor e mantida fixa a temperatura nos extremos da haste. Temperatura da Haste 0 40 80 120 160 200 240 280 x 5 x (cm) T ( G ra u s C el su s) 37 A.3. Discretização Uma configuração simples modeladora do sistema é uma série de elementos de igual comprimento (Fig. A.4b). Assim, o sistema é tratado com quatro elementos de igual comprimento e cinco nós. A.4. Equações dos Elementos Um elemento individual é mostrado na Fig. A.6a. A distribuição de temperatura para o elemento é representada pela função de aproximação: 2211 ~ TNTNT += (A.2) Onde N1 e N2 são funções lineares de interpolação especificadas pelas Eq. (3.48) e (3.49), respectivamente. Assim, como descrito na Fig. A.6b, a função de aproximação é uma interpolação linear entre as duas temperaturas dos nós. Nó 1 Nó 2 (a) T ~ T2 T1 x1 x2 (b) Figura A.6 (a) Um elemento individual (b). A função aproximação usada para caracterizar a distribuição de temperatura ao longo do elemento. A.5. Aproximação Direta Como descrito no subitem 3.2.1, há uma variedade de aproximações para desenvolver as equações elementos. Nesse subitem, emprega-se duas delas. Primeiro, uma aproximação direta será usada em um caso simples, onde f(x)=0. Então, em função de sua aplicabilidade geral na Engenharia, devotar-se-á a maior parte da abordagem para o método dos resíduos ponderados. No caso onde f(x)=0, o método direto é empregado para gerar as equações elementos. O relacionamento entre o fluxo de calor e o gradiente de temperatura é assim representado pela lei de Fourier: dx dT kq ′−=