O Método dos Elementos Finitos e os métodos “meshfree” na solução de Equações Diferenciais Parciais

1. O Método dos Elementos Finitos e os métodos “meshfree” na solução de Equações Diferenciais Parciais José Márcio Machado UNESP - Campus de São José do Rio Preto

2. 1) Ilustrando as idéias básicas do método com um exemplo simples: Consideremos a equação de Laplace em domínios bidimensionais, que aparece comumente em vários tipos de problemas. Pede-se calcular o potencial u que é solução da e.d.p. (1.1) em um domínio com uma fronteira C, dividida em subcontornos C1 e C2 onde u satisfaz respectivamente a condições de Neumann e de Dirichlet, e tais que O princípio de mínima energia para o potencial requer que a energia do campo armazenada no domínio assuma um valor mínimo para a solução de (1.1). Em outras palavras, a solução procurada minimiza a funcional (1.2) Portanto, uma forma de resolver (1.1) é procurar uma expressão aproximada para W(u), supondo que u possa ser aproximado por uma expansão com funções simples adequadamente escolhidas, e coeficientes indeterminados. A minimização da energia determinará então esses coeficientes, e portanto a aproximação para u. Esta “abordagem indireta” através de uma formulação integral equivalente e expansões de aproximação é, em essência, a idéia bá- sica seguida por todos os Métodos de Elementos Finitos e suas variantes: substituir o problema de encontrar diretamente a solução de uma dada equação diferencial parcial, pelo problema equivalente de encontrar uma solução minimizante para uma certa integral.

3. Ao construir essa expansão para o potencial u, inicialmente decompomos o domínio do problema em uma união de subdomínios topológicamente regulares (os elementos finitos), disjuntos dois a dois, não superpostos e não intersectantes (ou seja, construímos uma cobertura do domínio), como ilustrado pela figura 1 abaixo. Seja então, um elemento individual dessa decomposição (fig. 2): Fig. 2 No caso mais simples, uma aproximação U para a solução u de (1.1) é dada no elemento da Fig. 2 simplesmente por (1.3) com os coeficientes a, b e c satisfazendo Fig. 1 OBS.: Note na figura que os elementos po- dem possuir dimensões e orientações distin- tas, se adaptando às particularidades da geo- metria do domínio original. (1.4)

4. Combinando (1.3) e (1.4), é imediato concluir que para um único elemento, (1.5) Com um pouco de Álgebra, a expansão (1.5) pode ser escrita na forma (1.6) onde, por exemplo, sendo A a área do elemento triangular. As demais funções alfa são fácilmente obtidas por permutação cíclica de índices. Essas funções, chamadas na terminologia do método de fun- ções de base ou funções de forma, são interpolatórias nos vértices do triângulo, satisfazendo a propriedade do delta de Kronecker: (1.7) Substituindo agora a expansão (1.6) na equação (1.2), a energia associada a um únicoele- mento é aproximada por

5. (1.8) Se definirmos os elementos de matriz (1.9) a equação (1.8) será escrita então como uma forma matricial quadrática: (1.10) onde U representa o vetor coluna com os valores do potencial nos vértices (pontos nodais) do elemento triangular. Para determinar a aproximação que inclua todos os elementos de uma dada discretização (malha de elementos finitos), é suficiente impor a continuidade do potencial em pontos nodais que são comuns a dois ou mais elementos. Vamos considerar na figura ao lado, as numera- ções disjunta (a) e conjunta (b), que também são chamadas na terminologia corrente do método, respectivamente, de numerações local (a) e global (b). Fig. 3

6. Para os dois elementos separados, o vetor das incógnitas se escreverá como (1.11) e a energia total assocciada ao par de elementos será (1.12) onde (1.13) A matriz S em (1.13) é conhecida como matriz de Dirichlet ou de rigidez. Se considerarmos agora a montagem conjunta com numeração global (elementos unidos), é necessário que os potenciais variem contínuamente através das interfaces entre elementos. Em outras palavras, ao juntarmos dois elementos, os valores de potencial em vértices (pontos nodais) correspon- dentes devem ser iguais. Na figura 3.a, isto quer dizer que os potenciais em 1 e 6 são iguais, assim como os potenciais em 2 e 4.

7. Relacionando agora as numerações local e global por meio de uma matriz de conexões, (1.14) Desta forma, substituindo (1.14) em (1.12), teremos (1.15) (1.16) A eq. (1.16) define a matriz global para o problema, fornecendo juntamente com U uma aproximação para a energia do domínio constituído pelos dois elementos triangulares conectados.

8. O passo que falta é a minimização da energia W: (1.17) onde os subacritos l e p representam, respectivamente, os potenciais livres para variar (as in- cógnitas) e os valores prescritos (condições de contorno do tipo Dirichlet). Uma vez que a di- ferenciação óbviamente só é possível com respeito aos potenciais livres, temos do que resulta por fim, (1.18) Note que nesta formulação, as condições de Dirichlet são satisfeitas exatamente (são os valores prescritos de potencial). Pode-se mostrar que as condições de Neumann homogêneas (derivada normal nula) são implícitamente satisfeitas para este caso, não de forma precisa, mas com um valor aproximado que depende essencialmente da natureza da aproximação po- linomial que fizemos para U. Neste exemplo ela é linear nas coordenadas, mas poderia ser quadrática, cúbica, etc. Uma demonstração deste fato será feita na seção 2. Os pontos essenciais a mencionar agora são os seguintes: 1) Toda a técnica apresentada acima pode ser estendida sem dificuldade para qualquer número de elementos, aumentando apenas a dimensão dos vetores e matrizes. A formação da matriz global também pode ser feita diretamente, sem matrizes de conexão;

9. 2) A matriz global em (1.16) é esparsa e simétrica. O padrão de distribuição dos elementos nulos depende da numeração global adotada, que pode ser inteiramente arbitrária; 3) O grau do polinômio em (1.3) pode ser aumentado para se obter uma aproximação tão boa quanto se queira. Na aplicação do MEF, o que se busca é um compromisso: soluções precisas com aproximações lineares requerem um grande número de elementos, mas aproximações quadráticas ou cúbicas aumentam a dimensão das matrizes de contribuição local (1.9). Aplicações correntes raramente empregam aproximações cúbicas ou de grau maior. No exemplo anterior (P. P. Silvester & R. L. Ferrari, Finite Elements for Electrical Engineers, Cambridge University Press) estão presentes portanto os seguintes elementos característicos do MEF: A) Uma funcional, oriunda de princípios variacionais ou da integral de um resíduo ponderado (métodos de Galerkin, Bubnov-Galerkin, etc.) cuja solução minimizante é a solução da e.d.p. correspondente, B) Uma discretização do domínio segundo uma certa malha de elementos finitos; C) Um conjunto de funções de base, interpolantes nos elementos da discretização, satisfazendo as propriedades(1.7). O ítem A) será discutido mais detalhadamente agora. É importante adiantar que, mesmo para problemas em que nào é possível encontrar uma formulação variacional correspondente, o método de Galèrkin e suas variantes permite a construção da integral procurada.

10. 2) Princípios Variacionais e o Método de Galerkin : Usando o Cálculo das Variações, é possível mostrar a validade do seguinte resultado: Sejam funções dadas em uma região G bidimen- sional, e funções definidas do segmento de arco s ao longo da fronteira C de G. A função u(x,y) que torna a expressão integral estacionária sob a condição associada em onde C1 é um segmento de C, será necessáriamente uma solução do problema de valores de contorno em G, sob as condições de contorno de Dirichlet em C1 e geral de Cauchy em C2, sendo

11. nx e ny os cosenos diretores da normal exterior relativamente ao contorno C, e sendo C2 o resto do contorno, de tal forma que Claramente temos uma vantagem aqui: a condição do tipo Cauchy não está presente na formulação integral de modo explícito. Apenas a condição de Dirhclet precisa ser levada em consideração, se formulamos o problema através de um princípio de extremo. (H. R. Schwarz, Finite Element Methods, Academic Press) (equação de Laplace) (equação de Poisson) (equação de Helmholtz) Existem contudo, muitos problemas para os quais não é possível encontrar funcionais derivadas pelo Cálculo das Variações. Nestes casos, uma técnica bem mais abrangente para se encontrar a integral procurada é o chamado método de Galerkin, uma variante do Método dos Resíduos Ponderados. O método de Galerkin irá produzir matrizes globais formalmente idênticas às encontradas, por exemplo, pelo teorema anterior em problemas que admitem uma formulação variacional.

12. O método de Galerkin: Como já vimos, para um único elemento finito escrevemos a aproxima- ção onde supomos agora que a aproximação pode ser de ordem mais elevada, com p maior que 3. (2.1) Ora, a forma global da aproximação U em todo o domínio discretizado é uma combinação das a- proximações U(e)(x,y) sobre todos os elementos, ou seja, é a união das aproximações (2.1) sobre todos os elementos. Numerando todos os pontos nodais do domínio discretizado de 1 a n, a apro- ximação global U será dada então por (2.2) são as funções de base globais. Em (2.2), é a combinação (união) das funções de base individuais que tem o valor 1 no ponto nodal Pk, que é associado com a variável nodal Uk. Claramente a função de base global índice k tem valores não nulos apenas em elementos que contém o ponto nodal Pk, um subdomínio muito restrito, seu suporte local. A expressão (2.2) é uma aproximação de Ritz, com funções que possuem apenas suporte local. Contribuições integrais globais serão dadas naturalmente pela soma das contribuições individuais, por exemplo, Para exemplificar o método, consideremos como exemplo o problema de condução de calor em regime estacionário, dado pelas equações

13. da chamada forma mista, e em termos de operadores vetoriais, como (2.3) (2.4) Usando (2.2) para aproximar cada um dos componentes do vetor u em (2.3), podemos escrever então que (2.5) que é a integral da ponderação do resíduo tanto da equação (2.3) quanto de suas condições de contorno (2.4), resíduos estes obtidos ao substituirmos o valor exato da solução pela aproximação (2.2). Quando escolhemos as funções de ponderação wj como sendo as funções de base globais da aproximação (2.2), temos o método de Galerkin, ou de Bubnov-Galerkin.

14. Como podemos ver na equação (1.18), problemas de campo estacionário se reduzem em geral a uma equação matricial do tipo sendo a matriz S bem esparsa. Exemplo: discretização de um domínio, e matriz de Dirichlet global resultante: Fig. 2.2 Fig. 2.1 Matriz de Dirichlet global para funções de base quadráticas nos elementos triangulares. Numeração global arbitrária. OBS.: As regiões em branco na matriz são ocupadas por zeros.

15. 3) O Método EFGM (Element Free Galerkin Method) : • Na técnica EFGM, aproximantes MLS são usadas para construir uma aproximação uh( x ) para a solução exata y( x ) de uma dada equação diferencial parcial. Para tanto, um conjunto discreto de N pontos nodais é tomado sobre o domínio do problema, e uma função w de suporte compacto é associada a cada ponto nodal ou nó. Para cada nó existe um subdomínio, chamado domínio de influência, no qual o ponto nodal respectivo contribui para a solução. Este subdomínio consiste simplesmente da região em torno do ponto em que a função w (função pêso) é diferente de zero. Ao contrário do MEF, em que os elementos devem ser disjuntos dois a dois, domínios de influência na técnica EFGM podem se superpor. Para simplificar, consideremos um problema no plano e domínios de influência circulares ou retangulares, que correspondem respectivamente às duas funções pêso mais comumente encontradas, dadas por exponenciais gaussianas ou produtos tensoriais:

16. O domínio geométrico do problema é o contorno omega, e os círculos e retângulos superpostos são os domínios de influência. A solução procurada é avaliada nos pontos indicados, que continuaremos chamando de nodais mas que agora não fazem mais parte de uma malha de elementos conectados, mas são centrais a regiões que podem ser superpostas.

17. Uma aproximação local uhl(x) é construída como um polinomio de ordem m com coeficientes não constantes, como segue: (3.1) A escolha mais simples é uma base linear pT, em uma dimensão: Para calcular os coeficientes a ( x ) e resolver a aproximação acima, uma norma discreta ponderada L2 (denotada por J) é construída de tal forma que (3.2) onde a soma é extendida a todos os nós para os quais a função pêso w é diferente de zero, e ulé um parâmetro nodal que controla a função u(x) no ponto xl. Note que para a técnica EFGM, os parâmetros nodais ulsão diferentes dos valores uh(xl) da aproximação.

18. Calculando as derivadas parciais de J com respeito às componentes de a(x), resulta do critério de minimização que (3.3) onde as matrizes A e B e o vetor u, para uma base polinomial linear, são dados por

19. Por substituição de (3.3) em (3.1), resulta onde definimos (3.4) A abordagem EFGM consiste portanto em utilizar uma expansão como (3.4), e seus elementos constituintes, para a análise de uma dada EDP. Como um primeiro exemplo, consideremos a equação escalar de Helmholtz homogênea em um domínio bidimensional que descreve a secção reta de uma guia de onda retangular. Na análise modal da equação acima, uma expansão como (3.4) é usada na funcional correspondente, que é modificada quando condições de Dirichlet são impostas. Por exemplo, as soluções para modos TM na guia devem satisfazer condições de Dirichlet homogêneas no contorno do domínio. As soluções procuradas serão expressas por aproximantes EFGM, de modo que

20. Por exemplo, as soluções para modos TM na guia devem satisfazer condições de Dirichlet homogêneas no contorno do domínio. As soluções procuradas serão expressas por aproximantes EFGM, de modo que A funcional modificada incluirá multiplicadores de Lagrange para imposição das condições de contorno: (3.5) As condições de contorno homogêneas requerem que seja zero na fronteira, de modo que o problema de auto-estados correspondente para a funcional F é dado por ((5)

21. (3.6) onde Nas integrais acima, s é o comprimento de arco ao longo do contorno e as funções N(s) são interpolantes de Lagrange. Dependendo da geometria do problema, uma soma discreta também pode ser usada na imposição das CC: Neste caso, a expressão para G é simplesmente

22. Neste caso, a expressão para G é simplesmente A equação (3.6) pode ser então resolvida para determinar o vetor u. Para o caso de modos TE, condições de Dirichlet não são dadas explícitamente, o uso de multiplicadores de Lagrange é desnecessário e temos então que

23. Conhecidos os valores de u, usamos novamente a expansão (3.7) para determinar os valores do campo escalar. Observe que (3.7) fornece os auto-estados físicos do problema, as frequências correspondentes serão Note que o uso dos multiplicadores de Lagrange se deve ao fato de que as funções de forma vistas até aqui não são interpolantes, como no MEF, não satisfazendo a condição do delta de Kronecker.

24. A inconveniência disto é que os multiplicadores que aparecem em (5), embora essenciais para assegurar a satisfação das condições de contorno, não se relacionam diretamente com a solução física do problema. Isto representa um dispêndio extra de memória e processamento. Mesmo assim, resultados razoáveis podem ser conseguidos com esta abordagem. Como exemplo de desempenho, consideremos uma guia de onda retangular, com 10x20 mm, nós regularmente espaçados e integração por quadratura gaussiana 4x4 em cada célula. As soluções foram obtidas empregando-se funções de base bidimensionais lineares, tais como Como funções pêso, empregaram-se produtos tensoriais de splines cúbicas. Os gráficos seguintes ilustram as isolinhas de campo obtidas com o auxílio da técnica EFGM, e os êrros relativos no cálculos de alguns modos TE e TM (comparação com as soluções analíticas), feitos pela aplicação do MEF e da técnica EFGM, respectivamente. As figuras ilustrando as isolinhas demonstram a plena satisfação das condições de contorno.

25. Isolinhas para dois modos TE ..

26. Fig.1 Erros relativos nas frequencias, obtidos para uma guia retangular por EFGM e FEM em primeira e segunda ordens: (a) modo TE01, (b) modo TE23, (c) modo TE25, e (d) modo TE48 .

27. Fig.2. Erros relativos nas frequencias, obtidos para uma guia retangular por EFGM e FEM em primeira ordem, como função dos graus de liberdade: (a) modo TM11, (b) modo TM23, (c) modo TM25, e (d) modo TM55.

28. Outro exemplo: a equação de Schrödinger unidimensional Consideremos um domínio unidimensional dado por e no qual um poço ou barreira de potencial quântico é descrito por uma função de potencial arbitrária V(x). Na chamada aproximação de massa efetiva, a função de onda y para um elétron no domínio é dada pela equação de Schrödinger unidimensional (3.8) onde y é a função de onda do elétron, meff é a sua massa efetiva, E é o auto-valor da energia e Como antes, a funcional modificada por multiplicadores de Lagrange se escreve como

29. (3.9) Novamente, o problema de auto-valores se escreve como já que a função de onda (estados confinados) deve se anular nos extremos do domínio.

30. Casos de teste: 1) Poço de potencial com GaAs/AlxGa1-xAs: me= (0.067 + 0.083x)mo, x = 0.3

31. Resultados para U = 0.225 eV: Auto-valores de energia para diferentes modos ligados, em função da largura do poço:

32. Poço quântico de GaAs/AlxGa1-xAs, com campo elétrico • externo: me= (0.0665 + 0.0835x)mo, x = 0.32, d = 95 Angstrons, db= 98 Angstrons, U = 0.34 eV, com campos variando de 0.0 a 12x10**4 V/cm:

33. Estados fundamentais para elétrons e lacunas (N = 1) relatiamente aos potenciais no centro do poço quântico, como função do campo elétrico aplicado:

34. Análise de erros numéricos: e Erros dos auto-valores de energia para os modos N = 1, 2, 3, e 4 em função do número de pontos nodais.

35. Conclusão: Os resultados obtidos para estados ligados em estruturas com poços de potencial retangulares, e para estados quase-ligados em poços submetidos a campos elétricos externos apresentam boa concordância com dados obtidos na Literatura através de outros métodos. Isto justifica a continuação da pesquisa para simular dispositivos eletro-ópticos e componentes semicondutores pela técnica EFGM.. Possíveis melhorias incluem a utilização de aproximantes MLS interpolantes, que possibilitam a eliminação de multiplicadores de Lagrange nas funcionais das EDPs.

36. Aproximantes MLS interpolantes: Como já mencionamos, em geral aproximantes MLS do tipo (3.4) não interpolam dados, uma vez que a propriedade do delta Kronecker não é satisfeita para o conjunto de pontos nodais. Contudo, já foi demonstrado que esta propriedade pode ser satisfeita tomando-se funções w singulares na definição das funções aproximantes (3.4). A construção dessas funções pêso singulares pode ser feita da seguinte forma: I) Seleciona-se um conjunto original de funções pêso w tais que (3.10)

37. II) Constroi-se um novo conjunto de funções pêso singulares tomando-se o conjunto da etapa I), aplicando-se em seguida a definição (3.11) Claramente, isto resulta em podendo-se mostrar que então, Existe um processo recursivo para geração das funções de forma neste caso, conhecido na Literatura como abordagem de consistência (consistency approach).

38. Segundo os cálculos da abordagem de consistência, as funções de forma serão dadas agora por (3.12) Com as funções pêso modificadas, a expressão acima se escreve como

39. onde p, q e r são sempre diferentes de i. De acordo com as propriedades das funções pêso singulares, sempre que e portanto Contudo, funções pêso singulares podem resultar em um mau condicionamento do algoritmo. Esse problema é comumente superado na implementação computacional tomando-se aproximantes MLS com funções pêso de Shepard normalizadas. Neste caso, define-se na vizinhança de um nó uma nova classe de funções pêso, Wi, dada por

40. Que tipo de funções podemos tomar e que satisfazem a condição (3.10)? Há um grande número de escolhas possíveis para funções pêso mencionadas na Literatura, mas as mais comuns em geral são splines e funções de Schwarz. Se os domínios de influência são retangulares por exemplo, podemos tomar produtos tensoriais de funções de Schwarz em que, para cada coordenada, a variável r é a norma da diferença entre o valor da coordenada para o ponto de observação e para o j-ésimo ponto nodal: (3.13) e uma vez que fica claro que o produto tensorial de funções como (3.13) (no plano, uma função para x e outra para y) satisfaz a condição do delta de Kronecker.

41. 3) Um teste da formulação interpolante: A Fig.1 mostra um problema de teste (C-type magnet) resolvido pelo MEF e pela técnica EFGM com funções de forma interpolantes. Os parâmetros são Jz é a componente z do vetor densidade de corrente, e é a relutividade magnética do vácuo. Os pontos do domínio usados para a análise dos resultados obtidos com os dois métodos são mostrados na Fig. 2. Fig. 1 Fig. 2

42. As equações de Maxwell para campos magnetostáticos são reduzidas neste caso a uma única equação de Poisson para o potential vetor magnético, que devido à simetria do problema só tem a componente z: . (3.14) Uma vez que estamos usando aproximações EFGM interpolantes, multiplicadores de Lagrange são desnecessários para a imposição das condições de contorno. Obtemos agora o seguinte sistema de equações lineares: (3.15)

43. O domínio foi discretizado em 3773 pontos nodais, e as integrações numéricas executadas com o uso de quadratura Gaussiana de quarta ordem com 1856 células. Resultados de ensaios preliminares em problemas de eletrostática mostraram que uma discretização menos refinada mas com uma ordem de integração mais elevada é preferível, ao invés de um refinamento maior na discretização mas com ordens de integração mais baixas. A figura 3 abaixo mostra as isolinhas do potencial obtidas pela técnica EFGM para o problema em questão. Fig. 3

44. Fig. 4: Erros absolutos (tomando a solução FEM como referência) para os pontos de amostragem da Fig. 2. A aproximação FEM foi feita em 3779 pontos nodais e a EFGM, em 3773. Valores de erro absoluto em Wb/m2. Notar a razoável concordância entre os dois métodos. Como esperado, o erro é maior em regiões onde o fluxo de campo é mais intenso (núcleo de ferro, perto dos enrolamentos).

45. Fig. 5: Erros absolutos (tomando a solução FEM como referência) para os pontos de amostragem da Fig. 2. A aproximação FEM foi feita em 7993 pontos nodais e a EFGM, em 3773. Valores de erro absoluto em Wb/m2. Uma comparação dos erros relativos por sua vez, permite uma compreensão melhor do desempenho desses métodos:

46. Fig. 6: Erros relativos para os pontos da Fig. 2. A aproximação FEM foi feita em 3779 pontos nodais e a EFGM, em 3773. Um número de conectividade em torno de 6-11 mostrou-se uma boa escolha neste caso, mas a técnica de truncamento do domínio, devido a descontinuidade do material, pode requerer outras escolhas em outros casos. Note que os erros relativos são menores (melhores) para a solução FEM com melhor discretização (Fig. 7, adiante).

47. Fig. 7:Erros relativos para os pontos da Fig. 2. A aproximação FEM foi feita em 7993 pontos nodais e a EFGM, em 3773.

48. A Fig. mostra a descontinuidade da norma do campo magnético, H, calculada sobre a linha de avaliação da Fig. 1, que concorda razoávelmente bem tanto qualitativa quanto quantitativamente com o resultado obtido via FEM. Fig. 8: Descontinuidade na norma do campo segundo o método EFGM.

49. Outra aplicação: um problema de Magnetohidrodinâmica Fluxos magnetohidrodinâmicos em condição estacionária são normalmente descritos por um sistema de equações diferenciais parciais, composto pelas equações de Maxwell que descrevem o campo eletromagnético, e pelas equações de Navier-Stokes que descrevem o fluído condutor. Ou seja,

50. Afora outras hipóteses adicionais, variáveis adimensionais são introduzidas para simplificar as equações diferenciais: Após algumas manipulações algébricas, usando os índices 1 e 2 para denotar regiões de fluído e paredes do ducto respectivamente, e incluindo as variáveis adimensionais acima, o problema MHD será descrito por um conjunto de EDPs acopladas. Para a região de fluído: