Comparação entre o Método de Análise Isogeométrica e o Método dos Elementos Finitos

RESUMO

O Método de Análise Isogeométrica (AIG) é uma combinação entre o Método dos Elementos Finitos (MEF) e NURBS (Non-Uniform Rational B-Splines) onde as funções NURBS, além de representarem o domínio do problema, definem a base do espaço no qual aproximamos a solução de uma Equação Diferencial Parcial. Neste trabalho, a formulação básica do AIG é apresentada e discutiremos brevemente a formulação fraca do problema a partir da base de funções NURBS. Além disso, mostramos comparações numéricas entre o método dos elementos finitos e o método isogeométrico em termos de convergência. Por último, destacamos em um exemplo, uma das vantagens obtidas pelo AIG frente ao MEF devido à modelagem do domínio pelas funções NURBS. Um dos objetivos deste trabalho é servir como referência em português para uma introdução ao método AIG.

Palavras-chave:
Método Isogeométrico; Método dos Elementos Finitos; NURBS; Equações Diferenciais Parciais

ABSTRACT

The Isogeometric Analysis (IGA) is a combination between Finite Elements Method (FEM) and NURBS (Non-Uniform Rational B-Splines). The NURBS functionsmodel the problem domain and they define the basis of the space in which the Partial Differential Equation solution is approximated. In this article, we present the basic formulation of IGA and we discuss about the weak form of the problem. Moreover, we show numerical comparisons between the convergence of the finite elements method and the isogeometric analysis. Finally, we illustrate by an example, the advantages of modeling the domain using the NURBS functions, comparing both methods.

Keywords:
Isogeometric Analysis; Finite Elements Method; NURBS; Partial Differential Equations

1 INTRODUÇÃO

Desde que foi idealizado, por volta dos anos 50, o Método dos Elementos Finitos (MEF) tem passado por diversos avanços. Dentre os avanços mais recentes está o Método de Análise Isogeométrico (AIG), apresentado pela primeira vez por Hughes, Cottrell e Bazilevs em ¹⁴, cuja ideia principal é unir os conceitos CAD (Computer Aided Design) - MEF (Método dos Elementos Finitos) de maneira a modelar o domínio utilizando NURBS e então utilizar as mesmas NURBS como funções base para aproximação da solução do problema envolvendo Equações Diferenciais Parciais (EDP).

Em 2006, Bazilevs et al.⁴ apresentam a teoria e o método do h-refinimento no espaço utilizado pelo método isogeométrico, mostrando a convergência alcançada através desse refinamento. Em 2014, Veiga et al.²⁰ publicam na Acta Numérica uma revisão dos textos publicados sobre o AIG além da análise de convergência do método, inclusive utilizando T-splines, extensão das splines que permitem refinamento local (para detalhes ver ¹⁹).

Do ponto de vista dos problemas da engenharia, o método isogeométrico tem se destacado em muitas áreas onde podemos citar sua utilização para problemas de análise de placas e cascas (e.g. ver ¹⁵), mecânica não linear dos sólidos (e.g. ver ⁸), problemas de escoamento (e.g. ver ²), otimização de forma estrutural (e.g. ver ²⁴), interação fluido-estrutura (e.g. ver ³), problemas de eletromagnetismo (e.g. ver ⁵), entre outros.

A respeito dos trabalhos que comparam AIG com MEF, todos citam diferenças, vantagens edesvantagens, no entanto não apresentam exemplos que ilustram essas diferenças. Nosso principal objetivo aqui é ilustrar com tabelas e gráficos, geradas a partir de exemplo, as diferenças entre o AIG e MEF e trazer uma referência em português com detalhes sobre a implementação do AIG (ver ¹⁰).

Inicialmente apresentaremos de forma breve os conceitos principais das NURBS e da solução de EDPs elípticas utilizando o método isogeométrico, e então será apresentada uma comparação entre MEF e AIG. Para maiores detalhes veja ¹⁰.

2 B-SPLINES, NURBS, FUNÇÃO GEOMÉTRICA E REFINAMENTO

A principal diferença entre o método isogeométrico e o método dos elementos finitos está nas funções bases utilizadas pois, como dito anteriormente, o método isogeométrico utiliza asNURBS como funções base. As NURBS são coeficientes ponderados entre B-splines, que são funções polinomiais definidas por partes de forma recursiva da seguinte forma: Seja Ξ =(ξ_1, ξ_2,..., ξ_n ) um vetor de nós não decrescente. Definimos a i-ésima B-spline de grau p através das equações (2.1) e (2.2).

(2.1)

(2.2)

Para garantir que as funções base B-splines tenham a propriedade da partição da unidade é necessário a repetição dos nós inicial e final p + 1 vezes. Além disso, a repetição de um nó gera a redução do grau da função base correspondente, como ilustra a Figura 1.

Figura 1
Funções base B-splines de grau 3 geradas pelo vetor de nós (0, 0, 0, 0, 0.1, 0.2, 0.3, 0.3, 0.4, 0.5, 0.6, 0.7, 0.7, 0.7, 0.8, 0.9, 1, 1, 1, 1). Note que a cada vez que um nó se repete a função tem seu grau reduzido.

Definindo agora W = (w ₁, w ₂, ..., ^_wm ) um vetor de pesos (um peso para cada função base B-spline), podemos definir as funções base NURBS de grau p em uma dimensão através da equação (2.3):

(2.3)

As funções base NURBS bidimensionais são definidas por

(2.4)

onde _^Ni,p e _^Mj,q são funções base B-splines não necessariamente de mesmo grau, definidas a partir dos vetores de nós (ξ_1,..., ξ_n-p-1 ) e (η ₁,...,_{^ηm-q- 1}), respectivamente, p e q são os graus das funções base, e _^wij ’s são os pesos respectivos a cada função base.

Em ¹⁸ é apresentada uma dedução completa e em português das derivadas das B-splines, que são dadas pela equação abaixo

(2.5)

Já a derivada das NURBS em uma dimensão pode ser obtida através da derivada da função (2.3), de onde obtemos

(2.6)

e para duas dimensões, obtemos a partir das derivadas da equação (2.4), apresentadas nas equações (2.7) e (2.8).

(2.7)

(2.8)

Baseado na definição de derivada e observando a Figura 1 percebemos que nem sempre as derivadas das funções base estão bem definidas, porém, da forma que está definida a equação (2.5) temos que a derivada em um “bico” assume o valor da derivada à direita (isso vem da equação (2.1)).

Definidas as funções base NURBS, podemos definir a função geométrica, que é a função que define a geometria do domínio em função de um domínio de referência, chamado domínio paramétrico, ver Figura 2. Através das funções base NURBS e de pontos de controle, que são pontos no espaço ℝ² no caso do domínio bidimensional, a função geométrica é definida na equação (2.9), onde _^Pij é um ponto de controle dado da forma (x, y).

(2.9)

Figura 2
Exemplo de função geométrica (do domínio paramétrico Ω₀ no domínio geométrico Ω).

A utilização da função geométrica é a maior diferença entre o método isogeométrico e o método dos elementos finitos, pois ela dispensa a utilização de algoritmos para geração de malhas, de maneira que as NURBS utilizadas para definir o domínio são as mesmas utilizadas como funções base para solução do problema. Assim como as malhas geradas pela função geométrica estão diretamente ligadas às NURBS, o refinamento também está, e no caso do AIG temos três tipos de refinamentos: h-refinamento (refinamento do tamanho dos elementos), p-refinamento (aumento do grau de continuidade das bases) e k-refinamento (combinação dos dois métodos anteriores). Será apresentado aqui apenas o h-refinamento e mais detalhes sobre os outros tipos podem ser encontrados em ⁷. O h-refinamento se dá por uma sucessão de inserções de nós no vetor de nós, de maneira que estes nós são adicionados exatamente no meio entre os já existentes, e geram alteração nos pontos de controle e nos pesos de tal forma que a geometria do domínio é preservada.

Dado o vetor de nós (ξ₁, ξ₂, ..., ξ_n ) e os pontos de controle ponderados [_{^Pw 1}, _{^Pw 2}, ..., _{^Pwn-p- 1}], que são o produto dos pontos de controle por seus respectivos pesos, vamos inserir o nó ¹⁷), suponha que _k , ξ_{k+ 1}), então é necessário alterar os ^_Pwi , ∀_i ∈[k - p + 1, k]. Portanto, sejam [_{^Qw 1}, _{^Qw 2}, ..., ^_Qwn-p ] os novos pontos de controle ponderados após a inserção do nó . Conforme em (, temos que, ∀i ∈[k -p + 1, k] ∈ [ξ

Qwi = αi Pwi + (1 − αi)Pwi−1

onde,

Os pontos Qwi ∀i <k - p + 1 são iguais aos Pwi enquanto Qwj = Pwj -1, ∀j > k. O mesmo raciocínio dos pontos de controle é aplicado aos pesos para se encontrar os novos valores.

A Figura 3 ilustra o h-refinamento aplicado um vez. Note que apesar do número de elementos ter quadruplicado, o número de pontos de controle apenas duplicou, mostrando um comportamento típico do h-refinamento no AIG, que é o fato do número de elementos crescer muito mais rápido que o número de funções base. Para mais detalhes sobre NURBS ver (¹⁷),(¹⁸).

Figura 3
Exemplo de h-refinamento de um domínio bidimensional onde estão representadas as divisões dos elementos e os pontos de controle (círculos). A Figura superior mostra o domínio original e a inferior o mesmo domínio após passar por um refinamento.

3 O MÉTODO AIG APLICADO À EQUAÇÕES DIFERENCIAIS PARCIAIS

Para solução de Equações Diferenciais Parciais o método isogeométrico é semelhante ao método dos elementos finitos no que diz respeito à forma variacional de busca por soluções discretas, com diferença apenas na construção das funções base e na definição dos elementos. O AIG utiliza a função geométrica para o mapeamento, enquanto o MEF utiliza uma malha pré-definida.

3.1 Solução de Equações Diferenciais Parciais

Como exemplo, vamos mostrar a formulação variacional para encontrar a solução da equação de Poisson com condições de contorno de Dirichlet através do método isogeométrico. O problema é encontrar u: Ω → ℝ tal que

(3.1)

onde ∆u é o Laplaciano de u ou seja, ∆u = ∂ ² u/∂x ² + ∂ ² u/∂y ² e ∂Ω indica a fronteira de Ω.

A formulação variacional deste problema obtém-se multiplicando a primeira equação de (3.1) por uma função qualquer no clássico espaço H ¹ ₀(Ω) e integrando ambos os lados, obtendo

(3.2)

Pelo teorema da divergência temos

onde η(x, y) é um vetor normal exterior a Ω em (x, y). Então, como v(x, y) = 0 para (x, y) ∈ ∂Ω, temos que ∫_{∂ Ω} v(∇u.η)dS = 0, e assim o problema (3.2) tem a seguinte formulação variacional

(3.3)

Para discretizar o problema vamos definir um espaço de dimensão finita 𝒱^h tal que 𝒱^h está contido em H ¹ ₀(Ω) e 𝒱^h converge para H ¹ ₀(Ω) no sentido L ²(Ω) quando h converge para zero (veja ⁶ para maiores detalhes). A ideia é aproximar a solução exata que pertence a H ¹ ₀(Ω) por uma solução em 𝒱^h . Mais precisamente a solução encontrada em 𝒱^h é a projeção ortogonal na norma H ¹ da solução exata no espaço 𝒱^h . Definindo {ϕ ₁₁, ..., _^ϕnm } uma base para 𝒱^h e substituindo u e v por _{^{uh =}} ∑_{^n,mi,j =1 ^αijϕij} (x, y) e _^vh = ∑_{^n,mi,j =1 ^βijϕij} (x, y) em (3.3) temos,

(3.4)

Como (3.4) deve valer para qualquer função _^vh ∈ _^𝒱h temos que

(3.5)

para 0 ≤ k ≤ n e 0 ≤ l ≤ m.

As funções base NURBS são definidas para o domínio paramétrico e com isso faz-se necessário relacionar as funções base _^ϕij com as funções _^Rij definidas na equação (2.4). Essa relação éfeita através da função geométrica da seguinte forma:

(3.6)

Substituindo (3.6) em (3.5) temos

Fazendo agora uma mudança de variável do domínio geométrico (aquele que devemos resolver o problema) para o domínio paramétrico através da função geométrica F definida pela equação (2.9), temos

(3.7)

onde DF é a matriz jacobiana de F e J=det(DF(ξ, η)). A equação (3.7) vale para 0 ≤ k ≤ n e 0 ≤ l ≤ m.

Calcular essas integrais na maioria das vezes não é uma tarefa fácil. O método da Quadratura de Gauss tem se mostrado um método muito eficiente para a integração numérica e em particular é bastante eficiente na integração de NURBS ⁷.

A partir de (3.7) é possível montar um sistema linear onde as incógnitas são os _^αij ’s, e então aplicar as condições de contorno, da mesma maneira que se faz para o métodos dos elementos finitos, para assim obter a solução do sistema

(3.8)

3.2 Implementação Computacional

Vamos agora entender o processo de implementação do método isogeométrico para o problema apresentado na seção anterior, que se resume em resolver a equação, ou melhor, o sistema de equações apresentado em (3.7).

Primeiramente, por facilidade, vamos reescrever a equação (3.7) simplificando os índices.

(3.9)

para 0 ≤ j ≤ K e K = n ≤ m.

Podemos então escrever (3.9) de forma matricial.

(3.10)

onde

Vamos então calcular _^Aij para um elemento do domínio (lembrando que _^Aij é uma integral definida em todo domínio, porém devido ao suporte compacto das funções base, na maioria dos elementos seu valor é nulo). Para realizar a integração vamos utilizar o método da Quadratura de Gauss, conforme a equação (3.11).

(3.11)

onde _^xi é um ponto onde deve ser calculada f e _^pi é um peso pelo qual a mesma deve ser multiplicada (esse peso não tem nenhuma relação com os pesos utilizados nas funções NURBS).

Como nem todos os domínios de integração estão limitados entre [-1, 1], as vezes é necessário fazer uma mudança de variáveis, onde os novos valores de _^xi , agora chamados de b - a)/2 _^xi + (a + b)/2, onde a e b são os limites de integração inferior e superior, respectivamente. = (, assumem a forma

Para resolver uma integral dupla, que é o nosso caso, a ideia é a mesma aplicada a cada uma das variáveis. Como exemplo, podemos pensar na integração dentro de um dos elementos do domínio paramétrico através do método da Quadratura de Gauss usando 5 pontos de integração, que se resume na transformação de um domínio [-1, 1] × [-1, 1] para [_{^{ξi, ξi +}} 1] × [_{^{ηj, ηj +}} 1], conforme mostra a Figura 4. Com isso, a integração se tornou um somatório duplo e _^Aij agora tomou a forma apresentada na equação (3.12), onde e são os pontos da Quadratura de Gauss já transformados para o domínio de integração.

(3.12)

Figura 4
Transformação de domínio para integral dupla através do método da Quadratura de Gauss.

Para calcular DF(, ) basta seguir o seguinte procedimento:

1. Calcular as funções bases B-splines N(ξ) no ponto M(η) no ponto equações (2.1) e (2.2). e através das

2. Calcular as funções base NURBS bidimensionais de acordo com a equação (2.4).

3. Calcular as derivadas em ξ e η das funções base NURBS através da equação (2.7) e (2.8), respectivamente.

4. Calcular DF através das derivadas da equação (2.9).

Para calcular ∇_^Ri () basta seguir os passos de 1 a 3 do procedimento anterior, porém ao invés de calcular todas as funções base B-splines, é necessário calcular apenas as respectivas à função base Ri. Obviamente ∇Rj(, ) segue da mesma forma. Para calcular J basta calcular o determinante de DF, que já foi calculado.,

Calculados esses itens é possível calcular parte de Aij referente ao elemento em que estamos interessados. Resta calcular a parte respectiva aos outros elementos e somá-las. O cálculo dos bi’s segue de maneira análoga ao apresentado para Aij.

3.3 Erro e Convergência

Uma pergunta natural a se fazer é em que condições o AIG converge e qual é a taxa de convergência. Para o MEF essas questões já estão bem respondidas em (⁶). Em (⁷) temos uma ótima descrição da convergência do MEF para problemas elípticos: “No Método de Elementos Finitos clássico, a estimativa de erro para problemas elípticos lineares, expressa como a norma da diferença entre a solução exata, u, e a solução aproximada pelo método, uh, é da forma,

||u - uh||m ≤ Chβ ||u ||r

onde ||•||m e ||•||r são as normas correspondentes aos espaços de Sobolev Hm(Ω) e Hr(Ω), respectivamente, h está relacionado ao tamanho da malha, β =min(p + 1 - m, r - m), p é o grau da base polinomial, e C é uma constante que não depende de u ou h”.

Para o AIG temos um resultado semelhante, porém um pouco mais complexo tendo em vista que as funções base NURBS são mais complexas que as funções base do MEF. A descrição apresentada em (¹⁷) para a convergência do AIG é: “No Método de Análise Isogeométrica, a estimativa de erro para problemas elípticos lineares, é da forma,

onde k e l são índices inteiros tais que 0 ≤ k ≤ l ≤ p + 1, u ∈ H (Ω), he está relacionado ao tamanho de cada elemento, ∏k é o operador projeção ortogonal de Hk(Ω) no espaço gerado pelas bases da análise isogeométrica e C é uma constante que depende apenas de p e da forma do domínio (mas não do tamanho).”

Perceba que ∏ku é igual à uh definida em (3.8), uma vez que ambas são as projeções ortogonais relativas à norma Hk no espaço gerado pelas funções base ϕij. No entanto optamos por manter a notação apresentada em (⁷). Note também que a estimativa para o AIG é feita elemento a elemento e utiliza, da mesma forma que na resolução numérica de EDPs, a função geométrica para avaliar o domínio. As dificuldades na demonstração deste resultado se dão pelo fato das propriedades de aproximação das bases NURBS serem mais difíceis de determinar que das bases polinomiais convencionais (isto está relacionado com os pesos que são informações provenientes da geometria do domínio e não podem ser modificados). Além disso, as bases NURBS tem suporte em vários elementos e podem ter variação do grau de continuidade, o que dificulta estabelecer uma estimativa de convergência. A demonstração deste resultado, enunciado de forma mais precisa, pode ser encontrada em (⁴).

Vale destacar ainda que no AIG o somatório realizado elemento a elemento nos leva a crer que a convergência deste método é pior que a do MEF, porém Cottrell, Hughes e Bazilevs em (⁷) apresentam a seguinte conclusão: “A solução do Método de Análise Isogeométrica obtida usando NURBS de ordem p tem a mesma ordem de convergência que a esperada no Método dos Elementos Finitos clássico usando bases de funções clássicas com polinômios de ordem p.” Este resultado é extremamente importante, pois a cada h-refinamento, ou seja, dividindo todos os elementos ao meio, nós temos que o número de funções base do MEF aumenta muito mais do que do AIG (quanto maior o grau das funções base, maior essa diferença). Por exemplo, se tivermos apenas um elemento bidimensional de grau 2, ambos MEF e AIG possuem 9 funções base, porém ao passar por um h-refinamento, temos que o MEF passa a ter 25 funções base, enquanto o AIG apenas 16. Já para o grau 4, se tivermos apenas um elemento inicial, ambos tem 25 funções base e após apenas um refinamento, o MEF passar a ter 289, enquanto o AIG apenas 36. Esse fato implica em sistemas lineares com matrizes muitos menores para serem resolvidos no AIG, mas em contrapartida a implementação do AIG é mais cara computacionalmente para um mesmo número de funções base.

É importante destacar que a convergência tanto do MEF quanto do AIG dependem da regularidade da solução exata do problema, da conformidade do espaço gerado pelas funções base, e também da regularidade da solução exata. Uma boa referência sobre regularidade de soluções de problemas elípticos em domínios suaves é (⁹). O caso de domínios não suaves, principalmente polígonos e poliedros, é muito bem discutido em (¹²). A grosso modo a regularidade da solução de um problema elíptico depende da natureza da EDP, da regularidade da geometria do domínio, da regularidade do dados de fronteira e da regularidade do forçante, ou seja “o lado direito da EDP”.

Para ilustrar que o MEF e o AIG convergem na mesma taxa vamos comparar o erro da solução de uma EDO através dos dois métodos usando funções base de graus 1, 2 e 3. O erro percentual foi calculado na norma L2 de acordo com (3.13). Perceba que o exemplo a seguir tem objetivo apenas de ilustrar a convergência, pois o método AIG não tem efetividade em dimensão um, uma vez que a fronteira tem geometria óbvia.

(3.13)

Exemplo 1: Encontrar u: [0, 1] → ℝ tal que

A solução exata desse problema é u(x) = (10xsen(8πx) + x(x − 1))². Perceba que a solução exata é suave (C ^∞) e portanto as taxas de convergência do MEF e do AIG dependem da regularidade do espaço de funções teste.

A Figura 5 apresenta a comparação do erro entre o MEF e o AIG para funções base de graus 1, 2 e 3. Os dados utilizados para construir o gráfico se encontram na Tabela 1. Note que os primeiros dados da Tabela foram suprimidos no gráfico para uma melhor visualização.

Figura 5
Comparação do erro entre o MEF e o AIG para funções base de graus 1, 2 e 3 (log x log). Os triângulos tem hipotenusa de mesma inclinação que as retas e representam a taxa de convergência para cada grau.

Tabela 1
Comparação do erro entre o método dos elementos finitos e o método isogeométrico.

	Erro %
	Grau 1		Grau 2		Grau 3
h	MEF	AIG	MEF	AIG	MEF	AIG
1	7.474e-1	7.474e-1	1.191e+2	1.191e+2	1.191e+2	1.191e+2
2-1	1.128e+2	1.128e+2	1.590e+2	1.590e+2	1.590e+2	1.590e+2
2-2	4.186e+0	4.186e+0	3.917e+0	3.917e+0	3.917e+0	3.917e+0
2-3	9.921e-1	9.921e-1	2.863e-1	5.775e-1	2.850e-1	5.329e-1
2-4	1.462e-1	1.462e-1	1.125e-1	5.391e-1	6.156e-3	1.294e-1
2-5	1.360e-1	1.360e-1	1.272e-2	2.426e-2	1.218e-3	7.829e-3
2-6	3.296e-2	3.296e-2	1.488e-3	1.908e-3	7.199e-5	2.700e-4
2-7	8.195e-3	8.195e-3	1.800e-4	1.980e-4	4.531e-6	1.290e-5
2-8	2.052e-3	2.052e-3	2.280e-5	2.360e-5	2.854e-7	7.186e-7
2-9	5.160e-4	5.160e-4	2.842e-6	2.882e-6	1.783e-8	4.267e-8
2-10	1.260e-4	1.260e-4	3.547e-7	3.569e-7	1.115e-9	2.619e-9

Pode-se observar primeiramente que a partir de um h suficientemente pequeno os gráficos em escala (log x log) são retas. A inclinação destas determina a taxa de convergência. Para o grau 1 os dois métodos geram resultados exatamente iguais (de fato, para funções base de grau 1 os métodos são iguais) e têm taxa de convergência igual a 2, para grau 2 a partir de h < 0.01 as curvas se sobrepõem, mostrando uma taxa de convergência de mesma ordem, 3, e para grau 3 observa-se que para um mesmo h o erro do MEF é um pouco menor mas a taxa de convergência é a mesma do AIG e é igual a 4, pois as retas são paralelas, como previsto. Outra coisa importante a se destacar é que a medida que p aumenta a reta tem inclinação maior, ou seja convergêcia mais rápida.

3.4 Comparação entre o MEF e o AIG na Solução de um Problema Bidimensional

Tendo em vista que a vantagem proposta pelo método isogeométrico é o fato do mesmo mapear exatamente uma maior quantidade de domínios que as funções base clássicas do método dos elementos finitos (a saber o primeiro descreve com exatidão as cônicas e as quádricas), vamos aqui apresentar um problema em um domínio onde o AIG descreve exatamente a geometria, enquanto para o MEF o domínio precisa ser aproximado polinomialmente. Nosso exemplo foi pensado de forma que a solução exata tenha uma variação muito rápida perto da fronteira, para ilustrar a diferença entre o AIG e o MEF em termos de aplicação.

Exemplo 2: Encontrar u: Ω → ℝ tal que

onde o domínio é um disco de raio unitário, e f é igual à

Para cada a, a solução exata do problema é

Perceba que quanto mais próximo de 1 for o valor de a maior é a variação da solução próximo da fronteira, como podemos observar na Figura 6 para a = 1.20 e para a = 2.00. Além disso a solução é radial e suave para todo a > 1.

Para modelar o domínio foram utilizadas NURBS de grau 2, vetores de nós Ξ_{^{U =}} Ξ_V = (0, 0, 0, 1, 1, 1) e os pontos de controle e pesos apresentados na Tabela 2, que foram obtidos em ²³.

Figura 6
Solução exata para o problema. À esquerda está a solução para a = 1.20 e à direita a solução para a = 2.00.

Tabela 2
Pontos de Controle e pesos usando para modelar o domínio circular de raio unitário.

Pts de Controle	Pesos
	1
	1
(0,0)	1
	1
	1

Para o MEF, utilizando funções base de grau 1 e aproximando o domínio por polígonos de 10, 30 e 50 lados obtivemos os dados da Tabela 3 e geramos o gráfico da Figura 7. Observe que para o caso do polígono de 10 lados a aproximação é muito grosseira e isso gera erros muito grandes, porém mesmo para o caso onde o domínio foi aproximado por um polígono de 50 lados, o erro aumenta com o aumento do número de funções base, ao invés de diminuir. Esse fato é devido ao crescimento rápido da solução próximo da fronteira (para a = 1.20), o que faz com que o erro obtido próximo da fronteira seja grande, e pela natureza no problema o erro se propaga para todo o domínio, como mostra a Figura 8.

Tabela 3
Erro em função do número de funções base de grau 1 para o MEF em domínio aproximados por polígonos de 10, 30 e 50 lados (a = 1.20).

10 lados		30 lados		50 lados
fun. base	Erro(%)	fun. base	Erro(%)	fun. base	Erro(%)
14	123.11	160	27.34	194	15.20
43	111.53	607	37.98	723	16.26
149	113.29	2365	40.10	2789	16.70
553	116.68	9337	40.06	10953	16.46
2129	117.67
8353	117.89

Figura 7
Comparação do erro relativo obtido através de funções base de grau 1 do MEF para a = 1.20.

Figura 8
Erro absoluto obtido através de 8353 funções base de grau 1 do MEF para a = 1.20 em um domínio aproximado por polígono de 10 lados.

Para o método isogeométrico, o erro relativo obtido está apresentado na Tabela 4. Note que exceto pelo primeiro refinamento, o erro reduz sempre de um refinamento para o outro, como é esperado de acordo com os teoremas de convergência ^4),(7),(20.

Tabela 4
Análise erro relativo para o método isogeométrico (a = 1.20)

fun. base	Erro(%)
9	82.098
16	61.840
36	52.004
100	22.076
324	7.025
1156	1.645

Evidente que o erro obtido pelo MEF com funções base de grau 1 é muito maior do que o obtido pelo AIG, já que para o domínio circular é necessário o uso de funções base NURBS de grau 2. Vamos então verificar como se comporta o erro do MEF com funções base de grau 2.

Da mesma forma que para as funções de grau 1, temos que aproximar o domínio por um polígono, já que não é possível mapear exatamente com essas funções. Vale ressaltar que o que vamos chamar agora de lado do polígono é uma curva de grau dois formada pela interpolação de três pontos sobre a fronteira do domínio.

Na Tabela 5 estão apresentados os erros relativos para as funções base de grau 2 do MEF em domínios aproximados de 8, 12 e 16 partes (lados curvos dos elementos finitos de grau 2). A Figura 9 apresenta os dados dessa Tabela comparados com o erro relativo para o AIG apresentado anteriormente na Tabela 4.

Tabela 5
Erro relativo para o MEF com funções bases de grau dois em domínios aproximados de 8, 12 e 16 lados para a = 1.20.

8 lados		12 lados		16 lados
fun. base	Erro(%)	fun. base	Erro(%)	fun. base	Erro(%)
25	47.580	45	48.212	225	41.133
81	41.821	153	36.272	865	12.355
289	13.494	561	13.155	3393	2.874
1089	7.148	2145	5.612	13441	0.830
4225	7.929	8385	6.699	53505	0.384
16641	11.831	33133	13.871

Figura 9
Comparação do erro relativo obtido através de funções base de grau 2 do MEF e do AIG para a = 1.20.

É possível notar que para esses domínios aproximados o AIG tem o menor erro se compararmos o mesmo número de funções base (mesmo para um domínio de 40 lados, aqui não apresentado, o erro para 921 funções base de grau 2 é de 7,15% e para o próximo refinamento, 3601 funções base, é de 1,58%). Além disso, percebe-se que após certa quantidade de refinamentos o erro do MEF deixa de decrescer no caso dos domínios aproximados por polígonos de 8 e 12 lados, e para o domínio aproximado de 16 lados nota-se que a curva do erro teve uma ligeira redução na sua inclinação, indicando um possível comportamento semelhante aos outros casos.

Esse comportamento apresentado é devido ao fato da aproximação do domínio. Temos a garantia que o MEF converge para a solução exata, porém a solução exata no domínio aproximado não é a mesma solução do domínio exato, sendo assim, por mais refinamentos que se faça, o erro tende a ser a diferença entre a solução exata no domínio aproximado pela solução exata no domínio exato. Em muitos casos, essa diferença é muito pequena, mas em casos, assim como o exemplo apresentado, em que a solução cresce muito rápido próximo a fronteira e a fronteira não é bem aproximada, não se tem garantia a respeito da qualidade da solução obtida via MEF.

Como não conhecemos o comportamento das soluções dos problemas que desejamos resolver, é importante o uso de uma aproximação muito boa da fronteira, o que gera uma grande número de elementos no domínio e com isso um alto custo computacional. Sendo assim, utilizar o método isogeométrico em domínios mapeados exatamente por NURBS e não por polinômios é uma grande vantagem na convergência da solução. Um ponto importante a ser destacado é que o custo computacional do MEF é bem menor, para o mesmo número de funções base, do que o custo do AIG. Isso se deve essencialmente por conta da construção do sistema a ser resolvido, devido às integrais a serem calculadas. Sendo assim, não é possível a priori saber se é mais vantajoso em termos computacionais resolver o problema com uma aproximação muito fina com MEF ou usar o AIG para aqueles domínios em que podemos descrever de maneira exata o domínio via AIG.

Um fato importante a se destacar é que o AIG é de fato uma generalização do MEF, no seguinte sentido: Se o domínio pode ser descrito exatamente com polinômios então o AIG neste caso é o MEF. Outra coisa importante a se destacar é que não é fácil encontrar os pesos e pontos de controle para descrever um domínio dado, porém Kiendl et al. cita em ¹⁵ um plug-in para o programa Rhino (programa CAD baseado em NURBS utilizado por arquitetos e designers - http://www.rhino3d.com) capaz de transformá-lo em um pré-processador para o método isogeométrico.

O método AIG ganhou certa atenção da comunidade de métodos numéricos para solução EDPs, inclusive da comunidade de decomposição de domínios para desenvolvimento de pré-condicionadores para matrizes de sistemas oriundos da formulação do AIG, como exemplos temos ^13),(21),(22.

REFERÊNCIAS

• ¹
  R.A. Adams. “Sobolev Spaces”, Academic Press, New York (1975).
• ²
  Y. Bazilevs, V.M. Calo, J.A. Cottrell, T.J.R. Hughes, A. Reali & G. Scovazzi. Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows, Comp. Meth. in Appl. Mech. and Eng., 197 (2007), 173-201.
• ³
  Y. Bazilevs, V.M. Calo, T.J.R. Hughes & Y. Zhang. Isogeometric fluid-structure interaction: Theory, algorithms, and computations, Comp. Meth. in Appl. Mech. and Eng., 43 (2008), 3-37.
• ⁴
  Y. Bazilevs, L.B. da Veiga, J.A. Cottrell, T.J. Hughes & G. Sangalli. Isogeometric analysis: approximation, stability and error estimates for h-refined meshes, Math. Mod. and Meth. in Appl. Sci., 16 (2006), 1031-1090.
• ⁵
  A. Buffa, G. Sangalli & R. Vázquez. Isogeometric analysis in electromagnetics: B-splines approximation, Comp. Meth. in Appl. Mech. and Eng., 199 (2010), 1143-1152.
• ⁶
  P.G. Ciarlet. “The Finite Element Method for Elliptic Problems”, Siam, [S.l.], (2002).
• ⁷
  J.A. Cottrell, T.J. Hughes & Y. Bazilevs. “Isogeometric Analysis: Toward Integration of CAD and FEA”, John Wiley & Sons, [S.l.], (2009).
• ⁸
  T. Elguedj, Y. Bazilevs, V. Calo & T.J.R. Hughes. B, and F projection methods for nearly incompressible linear and non-linear elasticity and plasticity using higher-order NURBS elements, Comp. Meth. in Appl. Mech. and Eng., 197 (2008), 2732-2762.
• ⁹
  L.C. Evans. “Partial differential equations”, American Mathematical Society, Providence, Rhode Island (2010).
• ¹⁰
  H. Gomes Filho. “Método dos Elementos Finitos através da Análise Isogeométrica: Uma Introdução”, Dissertação de Mestrado, UFES, Vitória, ES (2016).
• ¹¹
  E. Gonçalves Júnior. “Preconditioners for Elliptic Control Problems”, Tese de Doutorado, IMPA, Rio de Janeiro, RJ (2009).
• ¹²
  P. Grisvard. “Elliptic Problems in Nonsmooth Domains”, SIAM (2011).
• ¹³
  C. Hofer & U. Langer. Dual-Primal Isogeometric Tearing and Interconnecting Solvers for large-scale systems of multipatch continuous Galerkin IgA equations, arXiv:1511.07183, (2015).
• ¹⁴
  T.J. Hughes, J.A. Cottrell & Y. Bazilevs. Isogeometric analysis: Cad, finite elements, NURBS, exact geometry and mesh refinement, Comp. Meth. in Appl. Mech. and Eng., 149 (2005), 4135-4195.
• ¹⁵
  J. Kiendl, K.-U. Bletzinger, J. Linhard & R. Wuchner. Isogeometric shell analysis with Kirchhoff-Love elements, Comp. Meth. in Appl. Mech. and Eng., 198 (2009), 3902-3914.
• ¹⁶
  V.P. Nguyen, C. Anitescu, S.P. Bordas & T. Rabszuk. Isogeometric analysis: an overview and computer implementation aspects, Math. and Comp. in Simulation, 117 (2015), 89-116.
• ¹⁷
  L. Pielg & W. Tiller. “The NURBS Book”, Springer Science & Business Media, [S.l.], (1996).
• ¹⁸
  F.F. Rocha. “NURBS e o Método Isogeométrico”, Dissertação de Mestrado, UFES, Vitória, ES (2016).
• ¹⁹
  M. Scott. “T-splines as a Design-Through-Analysis Technology”. Tese de Doutorado, The University of Texas at Austin (2011).
• ²⁰
  L.B. da Veiga, A. Buffa, G. Sangalli & R. Vázquez. Mathematical analysis of variational isogeometric methods, Acta Numerica, 23 (2014), 157-287.
• ²¹
  L.B. da Veiga, D. Cho, L.F. Pavarino & S. Scacchi. Overlapping additive Schwarz methods for iso geometric analysis, SIAM J. Num. Anal., 50 (2012), 1394-1416.
• ²²
  L.B. da Veiga, D. Cho, L.F. Pavarino & S. Scachi. BDDC preconditioners for isogeometric analysis, Math. Mod. and Meth. in Appl. Sci., 23, (2013), 1099-1142.
• ²³
  A.V. Vuong, C. Heinrich & B. Simeon. ISOGAT: A 2D tutorial MATLAB code for isogeometric analysis, Comp. Aid. Geom. Des., 27 (2010), 644-655.
• ²⁴
  W.A. Wall, M.A. Frenzel & C. Cyron. Isogeometric structural shape optimization, Comp. Meth. in Appl. Mech. and Eng., 197 (2008), 2976-2988.

Datas de Publicação

• Publicação nesta coleção
  Jan-Apr 2017

Histórico

• Recebido
  18 Out 2016
• Aceito
  01 Fev 2017