| | | |

Método dos elementos finitos

Ajude a manter o site livre, gratuito e sem propagandas. Colabore!

1.2 Problema modelo

Vamos fazer a introdução do método de elementos finitos pela aplicação no seguinte problema de valor de contorno (problema forte): encontrar u tal que

−u′′=f,x∈I=[0,l], (1.66)
u⁢(0)=u⁢(l)=0, (1.67)

onde f é uma função dada (fonte).

1.2.1 Formulação fraca

A derivação de um método de elementos finitos inicia-se da formulação fraca444Uma representação integral da equação diferencial. do problema em um espaço de funções apropriado. No caso do problema (1.66)-(1.67), tomamos o espaço

V0={v∈H1⁢(I):v⁢(0)=v⁢(l)=0}. (1.68)

Ou seja, v∈V0 é tal que v∈H1⁢(I), i.e. ‖v‖L2⁢(I)<∞ e ‖v′‖L2⁢(I)<∞, bem como v satisfaz as condições de contorno do problema.

A formulação fraca é obtida multiplicando-se a equação diferencial (1.66) por uma função teste v∈V0 (arbitrária) e integrando-se por partes, i.e.

∫If⁢v⁢𝑑x=−∫Iu′′⁢v⁢𝑑x (1.69)
=∫Iu′v′dx−u′(l)v(l)+u′(0)v(0) (1.70)

Donde, das condições de contorno das funções de teste v∈V0, temos

∫Iu′⁢v′⁢𝑑x=∫If⁢v⁢𝑑x. (1.72)

Desta forma, o problema fraco associado a (1.66)-(1.67) lê-se: encontrar u∈V0 tal que

a⁢(u,v)=L⁢(v),∀v∈V0, (1.73)

onde

a⁢(u,v)=∫Iu′⁢v′⁢𝑑x (1.74)
L⁢(v)=∫If⁢v⁢𝑑x, (1.75)

são chamadas de forma bilinear e forma linear, respectivamente. A formulação fraca construida aqui é chamada de formulação de Galerkin555Boris Galerkin, 1871 - 1945, engenheiro e matemático soviético. Fonte: Wikipédia: Boris Galerkin., em que o espaço da solução é igual ao espaço das funções teste.

1.2.2 Formulação de elementos finitos

Uma formulação de elementos finitos é um aproximação do problema fraco (1.73) em um espaço de dimensão finita. Aqui, vamos usar o espaço Vh⁢,0 das funções afins por partes em I que satisfazem as condições de contorno, i.e.

Vh⁢,0={v∈Vh:v⁢(0)=v⁢(l)=0}. (1.76)

Então, substituindo o espaço V0 pelo subespaço Vh⁢,0⊂V0 em (1.73), obtemos o seguinte problema de elementos finitos: encontrar uh∈Vh⁢,0 tal que

a⁢(uh,v)=L⁢(v),∀v∈Vh⁢,0. (1.77)

Observemos que o problema (1.77) é equivalente a: encontrar uh∈Vh⁢,0 tal que

a⁢(uh,φi)=L⁢(φi), (1.78)

para todo i=0,1,…,n, onde φi é a i-ésima função base de Vh⁢,0. Então, como uh∈Vh⁢,0, temos

uh=∑j=0nξj⁢φj, (1.79)

onde ξj, j=0,1,…,n, são as n+1 incógnitas a determinar. I.e., ao computarmos ξj, j=0,1,…,n, temos obtido a solução uh do problema de elementos finitos (1.77).

Agora, da forma bilinear (1.74), temos

a⁢(uh,φi)=a⁢(∑j=0nξj⁢φj,φi) (1.80)
=∑j=0nξja(φj,φi). (1.81)

Daí, o problema (1.77) é equivalente a resolvermos o seguinte sistema de equações lineares

A⁢𝝃=𝒃, (1.82)

onde A=[ai,j]i,j=0n é a matriz de rigidez com

ai,j=a⁢(φj,φi)=∫Iφj′⁢φi′⁢𝑑x, (1.83)

𝝃=(ξ0,ξ1,…,ξn) é o vetor das incógnitas e 𝒃=(b0,b1,…,bn) é o vetor de carga com

bi=L⁢(φi)=∫If⁢φi⁢𝑑x. (1.84)
Exemplo 1.2.1.

Consideramos o problema (1.66)-(1.67) com f≡10 e l=1, i.e.

−u′′=10,x∈I=[0,1], (1.85)
u⁢(0)=u⁢(1)=0. (1.86)

Neste caso, a solução analítica u⁢(x)=−5⁢x2+5⁢x pode ser facilmente obtida por integração. Vamos computar uma aproximação de elementos finitos no espaço das funções afins por partes Vh⁢,0={v∈Vh;v⁢(0)=v⁢(1)=0} construído numa malha uniforme de 5 células no intervalo I=[0,1]. Para tanto, consideramos o problema fraco associado: encontrar u∈V0={v∈H1⁢(I);v⁢(0)=v⁢(1)=0} tal que

a⁢(u,v)=L⁢(v), (1.87)

onde

a⁢(u,v)=∫Iu′⁢v′⁢𝑑x,L⁢(v)=∫I1⋅v⁢𝑑x. (1.88)

Então, a formulação de elementos finitos associada, lê-se: encontrar uh∈Vh⁢,0 tal que

a⁢(uh,vh)=L⁢(vh),∀vh∈Vh⁢,0. (1.89)

A Figura 1.5 apresenta o esboço dos gráficos da solução analítica u e da sua aproximação de elementos finitos uh (computada pelo Código 4). Verifique!

Refer to caption
Figura 1.5: Soluções aproximada uh versus analítica u=−5⁢x2+5⁢x do problema (1.85)-(1.86).
Código 4: ex_mef1d_modelo.py
1from mpi4py import MPI
2import numpy as np
3import ufl
4from dolfinx import mesh
5from dolfinx import default_scalar_type
6from dolfinx.fem.petsc import LinearProblem
7
8# malha
9domain = mesh.create_unit_interval(MPI.COMM_WORLD,
10 nx = 5)
11# espaço
12from dolfinx import fem
13V = fem.functionspace(domain, ('P', 1))
14
15# condição de contorno
16uD = fem.Function(V)
17uD.interpolate(lambda x: np.full(x.shape[1], 0.))
18
19tdim = domain.topology.dim
20fdim = tdim - 1
21domain.topology.create_connectivity(fdim, tdim)
22boundary_facets = mesh.exterior_facet_indices(domain.topology)
23boundary_dofs = fem.locate_dofs_topological(V, fdim,
24 boundary_facets)
25bc = fem.dirichletbc(uD, boundary_dofs)
26
27# problema mef
28u = ufl.TrialFunction(V)
29v = ufl.TestFunction(V)
30
31f = fem.Constant(domain, default_scalar_type(10.))
32
33a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
34L = f * v * ufl.dx
35
36problem = LinearProblem(a, L, bcs=[bc],
37 petsc_options_prefix="ex_mef1d_modelo")
38uh = problem.solve()

1.2.3 Estimativa a priori

As estimativas de erro são classificadas como a priori, quando o erro é dado em relação a solução u. E, quando expresso em relação à solução de elementos finitos uh, é classificada como a posteriori. A seguir, apresentamos algumas propriedades de uh que nos permitem obter uma estimativa a priori do erro ‖(u−uh)′‖L2⁢(I).

Teorema 1.2.1.(Ortogonalidade de Galerkin)

A solução de elementos finitos uh de (1.77) satisfaz a seguinte propriedade de ortogonalidade

a⁢(u−uh,v)=0,v∈Vh⁢,0, (1.90)

onde u é a solução de (1.73).

Demonstração.

De (1.77), (1.73) e lembrando que Vh⁢,0⊂V0, temos

a⁢(u,v)=L⁢(v)=a⁢(uh,v) (1.91)
⇒a(u−uh,v)=0, (1.92)

para todo v∈Vh⁢,0. ∎

Teorema 1.2.2.(A melhor aproximação)

A solução de elementos finitos uh dada por (1.77) satisfaz a seguinte propriedade de melhor aproximação

‖(u−uh)′‖L2⁢(I)≤‖(u−v)′‖L2⁢(I),v∈Vh⁢,0, (1.93)

onde u é a solução de (1.73).

Demonstração.

Escrevendo u−uh=u−v+v−uh para qualquer v∈Vh⁢,0 e usando a ortogonalidade de Galerkin (Teorema 1.2.1), temos

‖(u−uh)′‖L2⁢(I)2=∫I(u−uh)′⁢(u−uh)′⁢𝑑x (1.94)
=∫I(u−uh)′(u−v+v−uh)′dx (1.95)
=∫I(u−uh)′(u−v)′dx+∫I(u−uh)′(v−uh)′dx (1.96)
=∫I(u−uh)′(u−v)′dx (1.97)
≤∥(u−uh)′∥L2⁢(I)∥(u−v)′∥L2⁢(I). (1.98)

∎

Teorema 1.2.3.(Estimativa a priori)

O erro em se aproximar a solução u de (1.73) pela solução de elementos finitos uh dada por (1.77) satisfaz a seguinte estimativa a priori

‖(u−uh)′‖L2⁢(I)2≤C⁢∑i=1nhi2⁢‖u′′‖L2⁢(Ii)2. (1.99)
Demonstração.

Tomando v=π⁢u no teorema da melhor aproximação (Teorema 1.2.2), obtemos

‖(u−uh)′‖L2⁢(I)≤‖(u−π⁢u)′‖L2⁢(I). (1.100)

Daí, da estimativa do erro de interpolação (Proposição 1.1.2), temos

‖(u−uh)′‖L2⁢(I)2≤C⁢∑i=1nhi2⁢‖u′′‖L2⁢(Ii)2. (1.101)

∎

Observação 1.2.1.(Estimativa do erro na norma L⁢2)

Lembrando que e⁢(0)=e⁢(1)=0 para o problema modelo (1.66)-(1.67), a desigualdade de Poincaré666Henri Poincaré, 1854 - 1912, matemático francês. Fonte: Wikipédia: Henri Poincaré. garante que

‖u−uh‖L2⁢(I)≤C⁢‖(u−uh)′‖L2⁢(I). (1.102)

Assim, a estimativa a priori do erro ‖u−uh‖L2⁢(I) pode ser obtida a partir da estimativa (1.99), ou seja,

‖u−uh‖L2⁢(I)≤C⁢h⁢‖u′′‖L2⁢(I), (1.103)

assumindo uma malha uniforme, i.e. hi=h para todo i=1,2,…,n. Verifique!

Como observaremos no próximo Exemplo 1.2.2, esta não é uma estimativa a priori ótima. Na verdade, para elementos finitos P1, pode-se mostrar que [3, Section 5.4]

‖u−uh‖L2⁢(I)≤C⁢h2⁢‖u′′‖L2⁢(I), (1.104)

sob hipóteses adicionais sobre a regularidade da solução u.

Exemplo 1.2.2.

A Tabela 1.3 apresenta uma análise numérica da convergência de malha na solução de elementos finitos do problema (1.85)-(1.86).

Tabela 1.3: Análise da convergência de malha na solução de elementos finitos do problema (1.85)-(1.86).
n ‖u−uh‖L2⁢[0,1] taxa ‖(u−uh)′‖L2⁢[0,1] taxa
5 3.65⁢e−2 — 5.77⁢e−1 —
10 9.13⁢e−3 2.00 2.89⁢e−1 1.00
20 2.28⁢e−3 2.00 1.44⁢e−1 1.00
40 5.71⁢e−4 2.00 7.22⁢e−2 1.00

O cálculo dos erros ‖u−uh‖L2⁢(I) e ‖(u−uh)′‖H1⁢(I) pode ser feito pelo seguinte código, considerando uh computado pelo Código 4. Verifique!

Código 5: ex_mef1d_modelo_erro.py
1# solução analítica
2V2 = fem.functionspace(domain, ('P', 2))
3ua = fem.Function(V2)
4ua.interpolate(lambda x: 5. * x[0] * (1 - x[0]))
5
6# erro L2
7error_L2_local = fem.assemble_scalar(fem.form((uh - ua)**2 * ufl.dx))
8error_L2 = np.sqrt(domain.comm.allreduce(error_L2_local, op=MPI.SUM))
9print(f"Erro L2: {error_L2:.2e}")
10
11# erro H1
12error_H1_local = fem.assemble_scalar(fem.form(ufl.dot(ufl.grad(uh - ua), ufl.grad(uh - ua)) * ufl.dx))
13error_H1 = np.sqrt(domain.comm.allreduce(error_H1_local, op=MPI.SUM))
14print(f"Erro H1: {error_H1:.2e}")

1.2.4 Estimativa a posteriori

Vamos obter uma estimativa a posteriori para o erro e=u−uh da solução de elementos finitos uh do problema modelo (1.66)-(1.67).

Teorema 1.2.4.(Estimativa a posteriori)

A solução de elementos finitos uh satisfaz

‖(u−uh)′‖L2⁢(I)2≤C⁢∑i=1nηi2⁢(uh), (1.105)

onde ηi⁢(uh) é chamado de elemento residual e é dado por

ηi⁢(uh)=hi⁢‖f+uh′′‖L2⁢(Ii). (1.106)
Demonstração.

Tomando e=u−uh e usando a ortogonalidade de Galerkin (Teorema 1.2.1) temos

‖e′‖L2⁢(I)2=∫Ie′⁢(e−π⁢e)′⁢𝑑x (1.107)
=∑i=1n∫Iie′(e−πe)′dx. (1.108)

Então, aplicando integração por partes

‖e′‖L2⁢(I)2=∑i=1n∫Ii(−e′′)⁢(e−π⁢e)⁢𝑑x+[e′⁢(e−π⁢e)]xi−1xi. (1.109)

Daí, observando que e−π⁢e=0 nos extremos dos intervalos Ii e que −e′′=−(u−uh)′′=−u′′+uh′′=f+uh′′, temos

‖e′‖L2⁢(I)2=∑i=1n∫Ii(f+uh′′)⁢(e−π⁢e)⁢𝑑x. (1.110)

Agora, usando as desigualdades de Cauchy-Schwarz e a estimativa padrão de interpolação (1.28), obtemos

‖e′‖L2⁢(I)2≤∑i=1n‖f+uh′′‖L2⁢(Ii)⁢‖e−π⁢e‖L2⁢(Ii)⁢d⁢x (1.111)
≤C∑i=1nhi∥f+uh′′∥L2⁢(Ii)∥e′∥L2⁢(Ii) (1.112)
≤C(∑i=1nhi2∥f+uh′′∥L2⁢(Ii)2)1/2(∑i=1n∥e′∥L2⁢(Ii)2)1/2 (1.113)
=C(∑i=1nhi2∥f+uh′′∥L2⁢(Ii)2)1/2∥e′∥L2⁢(I), (1.114)

donde segue o resultado desejado. ∎

Observação 1.2.2.

No caso da solução de elementos finitos no espaço das funções lineares por partes, temos uh′′=0. Logo, o elemento residual se resume em ηi⁢(uh)=hi⁢‖f‖L2⁢(Ii).

1.2.5 Exercícios

E. 1.2.1.

Obtenha uma aproximação por elementos finitos P1 da solução de

−u′′=π2⁢sen⁡(π⁢x),∀x∈(0,1), (1.115)
u⁢(0)=u⁢(1)=0. (1.116)

Faça o esboço dos gráficos da solução analítica u=sen⁡(π⁢x) e da sua aproximação de elementos finitos uh. Então, faça uma análise da convergência de malha e verifique as taxas de convergência dos erros ‖u−uh‖L2⁢(I) e ‖(u−uh)′‖H1⁢(I).

E. 1.2.2.

Obtenha uma aproximação por elementos finitos P2 da solução do problema dado no E.1.2.1. Faça o esboço dos gráficos da solução analítica u=sen⁡(π⁢x) e da sua aproximação de elementos finitos uh. Então, faça uma análise da convergência de malha e estime as taxas de convergência dos erros ‖u−uh‖L2⁢(I) e ‖(u−uh)′‖H1⁢(I).

E. 1.2.3.

Obtenha uma aproximação por elementos finitos P⁢1 da solução de

−u′′+u=2⁢sen⁡x,∀x∈(−π,π), (1.117)
u⁢(−π)=u⁢(π)=0. (1.118)

Faça o esboço dos gráficos da solução analítica u=sen⁡x e da sua aproximação de elementos finitos uh. Então, faça uma análise da convergência de malha e estime as taxas de convergência dos erros ‖u−uh‖L2⁢(I) e ‖(u−uh)′‖H1⁢(I).

E. 1.2.4.

Considere o seguinte problema de difusão-advecção

−ϵ⁢u′′+u′=10,x∈(0,1), (1.119)
u⁢(0)=u⁢(1)=0, (1.120)

onde ϵ>0. A solução analítica deste problema é dada por

u⁢(x)=10⁢x−10⁢(ex/ϵ−1)e1/ϵ−1. (1.121)

Estude aproximações por elementos finitos para este problema com diferentes valores de ϵ. O que acontece com a solução de elementos finitos quando ϵ é muito pequeno? Faça uma análise da convergência de malha para ϵ≪1.

E. 1.2.5.

Estude aproximações por elementos finitos para o seguinte problema de difusão-advecção-reação com coeficientes variáveis

−[(1+x)⁢u′]′+x⁢u′+u=f⁢(x),x∈(0,1), (1.122)
u⁢(0)=u⁢(1)=0, (1.123)

onde f⁢(x)=π⁢(x−1)⁢cos⁡(π⁢x)+[1+π2⁢(1+x)]⁢sen⁡(π⁢x). A solução analítica deste problema é dada por u⁢(x)=sen⁡(π⁢x). Faça uma análise da convergência de malha e estime as taxas de convergência dos erros ‖u−uh‖L2⁢(I) e ‖(u−uh)′‖H1⁢(I).


Envie seu comentário

Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

Opcional. Preencha seu nome para que eu possa lhe contatar.
Opcional. Preencha seu e-mail para que eu possa lhe contatar.
As informações preenchidas são enviadas por e-mail para o desenvolvedor do site e tratadas de forma privada. Consulte a política de uso de dados para mais informações.

Licença Creative Commons
Este texto é disponibilizado nos termos da Licença Creative Commons Atribuição-CompartilhaIgual 4.0 Internacional. Ícones e elementos gráficos podem estar sujeitos a condições adicionais.

Pedro H A Konzen
| | | |