Método dos Elementos Finitos

1 Problemas unidimensionais 1.1 Interpolação e projeção 1.3 Condições de contorno

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

	$\displaystyle-u^{\prime\prime}=f,\quad x\in I=[0,l],$		(1.63)
	$\displaystyle u(0)=u(l)=0,$		(1.64)

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 fraca⁴⁴4Uma representação integral da equação diferencial. do problema em um espaço de funções apropriado. No caso do problema (1.63)-(1.64), tomamos o espaço

V_{0}=\{v\in H^{1}(I):\leavevmode\nobreak\ v(0)=v(l)=0\}.

(1.65)

Ou seja, $v\in V_{0}$ é tal que $v\in H^{1}(I)$ , i.e. $\|v\|_{L^{2}(I)}<\infty$ e $\|v^{\prime}\|_{L^{2}(I)}<\infty$ , bem como $v$ satisfaz as condições de contorno do problema.

A formulação fraca é obtida multiplicando-se a equação diferencial (1.63) por uma função teste $v\in V_{0}$ (arbitrária) e integrando-se por partes, i.e.

	$\displaystyle\int_{I}fv\,dx=-\int_{I}u^{\prime\prime}v\,dx$		(1.66)
	$\displaystyle\text{}\quad=\int_{I}u^{\prime}v^{\prime}\,dx-u^{\prime}(l)v(l)+u% ^{\prime}(0)v(0)$		(1.67)

Donde, das condições de contorno das funções de teste $v\in V_{0}$ , temos

\int_{I}u^{\prime}v^{\prime}\,dx=\int_{I}fv\,dx.

(1.69)

Desta forma, o problema fraco associado a (1.63)-(1.64) lê-se: encontrar $u\in V_{0}$ tal que

a(u,v)=L(v),\quad\forall v\in V_{0},

(1.70)

onde

	$\displaystyle a(u,v)=\int_{I}u^{\prime}v^{\prime}\,dx$		(1.71)
	$\displaystyle L(v)=\int_{I}fv\,dx,$		(1.72)

são chamadas de forma bilinear e forma linear, respectivamente. A formulação fraca construida aqui é chamada de formulação de Galerkin⁵⁵5Boris Galerkin, 1871 - 1945, engenheiro e matemático soviético. Fonte: Wikipédia., 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.70) em um espaço de dimensão finita. Aqui, vamos usar o espaço $V_{h,0}$ das funções afins por partes em $I$ que satisfazem as condições de contorno, i.e.

V_{h,0}=\{v\in V_{h}:\leavevmode\nobreak\ v(0)=v(l)=0\}.

(1.73)

Então, substituindo o espaço $V_{0}$ pelo subespaço $V_{h,0}\subset V_{0}$ em (1.70), obtemos o seguinte problema de elementos finitos: encontrar $u_{h}\in V_{h,0}$ tal que

a(u_{h},v)=L(v),\quad\forall v\in V_{h,0}.

(1.74)

Observemos que o problema (1.74) é equivalente a: encontrar $u_{h}\in V_{h,0}$ tal que

a(u_{h},\varphi_{i})=L(\varphi_{i}),

(1.75)

para todo $i=0,1,\dotsc,n$ , onde $\varphi_{i}$ é a $i$ -ésima função base de $V_{h,0}$ . Então, como $u_{h}\in V_{h,0}$ , temos

u_{h}=\sum_{j=0}^{n}\xi_{j}\varphi_{j},

(1.76)

onde $\xi_{j}$ , $j=0,1,\dotsc,n$ , são as $n+1$ incógnitas a determinar. I.e., ao computarmos $\xi_{j}$ , $j=0,1,\dotsc,n$ , temos obtido a solução $u_{h}$ do problema de elementos finitos (1.74).

Agora, da forma bilinear (1.71), temos

	$\displaystyle a(u_{h},\varphi_{i})=a\left(\sum_{j=0}^{n}\xi_{j}\varphi_{j},% \varphi_{i}\right)$		(1.77)
	$\displaystyle\text{}\quad=\sum_{j=0}^{n}\xi_{j}a(\varphi_{j},\varphi_{i}).$		(1.78)

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

A\boldsymbol{\xi}=\boldsymbol{b},

(1.79)

onde $A=[a_{i,j}]_{i,j=0}^{n}$ é a matriz de rigidez com

a_{i,j}=a(\varphi_{j},\varphi_{i})=\int_{I}\varphi_{j}^{\prime}\varphi_{i}^{% \prime}\,dx,

(1.80)

$\boldsymbol{\xi}=(\xi_{0},\xi_{1},\dotsc,\xi_{n})$ é o vetor das incógnitas e $\boldsymbol{b}=(b_{0},b_{1},\dotsc,b_{n})$ é o vetor de carga com

b_{i}=L(\varphi_{i})=\int_{I}f\varphi_{i}\,dx.

(1.81)

Exemplo 1.2.1.

Consideramos o problema (1.63)-(1.64) com $f\equiv 10$ e $l=1$ , i.e.

	$\displaystyle-u^{\prime\prime}=10,\quad x\in I=[0,1],$		(1.82)
	$\displaystyle u(0)=u(1)=0.$		(1.83)

Neste caso, a solução analítica $u(x)=-5x^{2}+5x$ pode ser facilmente obtida por integração. Vamos computar uma aproximação de elementos finitos no espaço das funções afins por partes $V_{h,0}=\{v\in V_{h};\leavevmode\nobreak\ 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\in V_{0}=\{v\in H^{1}(I);\leavevmode\nobreak\ v(0)=v(1)=0\}$ tal que

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

(1.84)

onde

a(u,v)=\int_{I}u^{\prime}v^{\prime}\,dx,\quad L(v)=\int_{I}1\cdot v\,dx.

(1.85)

Então, a formulação de elementos finitos associada, lê-se: encontrar $u_{h}\in V_{h,0}$ tal que

a(u_{h},v_{h})=L(v_{h}),\quad\forall v_{h}\in V_{h,0}.

(1.86)

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 $u_{h}$ (computada pelo Código 4). Verifique!

Refer to caption — Figura 1.5: Soluções aproximada $u_{h}$ versus analítica $u=-5x^{2}+5x$ do problema (1.82)-(1.83).

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

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))

15# condição de contorno

16uD = fem.Function(V)

17uD.interpolate(lambda x: np.full(x.shape[1], 0.))

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)

27# problema mef

28u = ufl.TrialFunction(V)

29v = ufl.TestFunction(V)

31f = fem.Constant(domain, default_scalar_type(10.))

33a = ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx

34L = f * v * ufl.dx

36problem = LinearProblem(a, L, bcs=[bc],

37 petsc_options_prefix="ex_mef1d_modelo")

38uh = problem.solve()

1.2.3 Estimativa a priori

Em revisão

Existem dois tipos de estimativas do erro $e:=u-u_{h}$ . Estimativas a priori, são aquelas em que o erro é dado em relação da solução $u$ , enquanto que nas estimativas a posteriori o erro é expresso em relação a solução de elementos finitos $u_{h}$ .

Teorema 1.2.1.

(Ortogonalidade de Galerkin.) A solução de elementos finitos $u_{h}$ de (1.74) satisfaz a seguinte propriedade de ortogonalidade

a(u-u_{h},v):=\int_{I}(u-u_{h})^{\prime}v^{\prime}\,dx=0,\quad v\in V_{h,0},

(1.87)

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

Demonstração.

De (1.74), (1.70) e lembrando que $V_{h,0}\subset V_{0}$ , temos

a(u,v)=L(v)=a(u_{h},v)\Rightarrow a(u-u_{h},v)=0,

(1.88)

para todo $v\in V_{h,0}$ . ∎

Teorema 1.2.2.

(A melhor aproximação.) A solução de elementos finitos $u_{h}$ dada por (1.74) satisfaz a seguinte propriedade de melhor aproximação

\|(u-u_{h})^{\prime}\|_{L^{2}(I)}\leq\|(u-v)^{\prime}\|_{L^{2}(I)},\quad v\in V% _{h,0},

(1.89)

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

Demonstração.

Escrevendo $u-u_{h}=u-v+v-u_{h}$ para qualquer $v\in V_{h,0}$ e usando a ortogonalidade de Galerkin (Teorema 1.2.1), temos

$\displaystyle\\|(u-u_{h})^{\prime}\\|_{L^{2}(I)}^{2}$	$\displaystyle=\int_{I}(u-u_{h})^{\prime}(u-u_{h})^{\prime}\,dx$	(1.90)
	$\displaystyle=\int_{I}(u-u_{h})^{\prime}(u-v+v-u_{h})^{\prime}\,dx$	(1.91)
	$\displaystyle=\int_{I}(u-u_{h})^{\prime}(u-v)^{\prime}\,dx+\int_{I}(u-u_{h})^{% \prime}(v-u_{h})^{\prime}\,dx$	(1.92)
	$\displaystyle=\int_{I}(u-u_{h})^{\prime}(u-v)^{\prime}\,dx$	(1.93)
	$\displaystyle\leq\\|(u-u_{h})^{\prime}\\|_{L^{2}(I)}\\|(u-v)^{\prime}\\|_{L^{2}(I)}.$	(1.94)

∎

Teorema 1.2.3.

(Estimativa a priori.) O erro em se aproximar a solução $u$ de (1.70) pela solução de elementos finitos $u_{h}$ dada por (1.74) satisfaz a seguinte estimativa a priori

\|(u-u_{h})^{\prime}\|_{L^{2}(I)}^{2}\leq C\sum_{i=1}^{n}h_{i}^{2}\|u^{\prime% \prime}\|_{L^{2}(I_{i})}^{2}.

(1.95)

Demonstração.

Tomando $v=\pi u$ no teorema da melhor aproximação (Teorema 1.2.2), obtemos

\|(u-u_{h})^{\prime}\|_{L^{2}(I)}\leq\|(u-\pi u)^{\prime}\|_{L^{2}(I)}.

(1.96)

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

\|(u-u_{h})^{\prime}\|_{L^{2}(I)}^{2}\leq C\sum_{i=1}^{n}h_{i}^{2}\|u^{\prime% \prime}\|_{L^{2}(I_{i})}^{2}.

(1.97)

∎

Exemplo 1.2.2.

A Figura 1.6 apresenta o esboço da evolução do erro $\|(u-u_{h})^{\prime}\|_{L^{2}(I)}$ da solução de elementos finitos do problema (1.82)-(1.83) para malhas uniformes com $n=2,4,8,\ldots,128$ células.

Com o FEniCS, a computação do problema de elementos finitos pode ser feita com o seguinte código:

from __future__ import print_function, division
from fenics import *
import numpy as np
import matplotlib.pyplot as plt

def boundary(x,on_boundary):
    return on_boundary

def solver(n):
    # malha
    mesh = IntervalMesh(n,0,1)

    # espaco
    V = FunctionSpace(mesh, ’P’, 1)

    bc = DirichletBC(V,Constant(0.0),boundary)

    #MEF problem
    u = TrialFunction(V)
    v = TestFunction(V)
    f = Constant(1.0)
    a = u.dx(0)*v.dx(0)*dx
    L = f*v*dx

    #computa a sol
    u = Function(V)
    solve(a == L, u, bc)

    return u, mesh

#sol analitica
ua = Expression(’-x[0]*x[0]/2+x[0]/2’,
                degree=2)


lerrors=[]
for n in [2,4,8,16,32,64,128]:
    u, mesh = solver(n)
    e = errornorm(u,ua,norm_type=’H10’,mesh=mesh)
    lerrors.append(e)

plt.plot([2,4,8,16,32,64,128],lerrors)
plt.xscale(’log’,basex=2)
#plt.yscale(’log’,base=2)
plt.xlabel(r"$n$")
plt.ylabel(r"$|\!|(u-u_h)’|\!|_{L^2(I)}$")
plt.xlim((2,128))
plt.xticks([2,4,8,16,32,64,128],[2,4,8,16,32,64,128])
plt.grid(’on’)
plt.show()

1.2.4 Estimativa a posteriori

Em revisão

Aqui, vamos obter uma estimativa a posteriori para o erro $e=u-u_{h}$ da solução de elementos finitos $u_{h}$ do problema (1.63)-(1.64).

Teorema 1.2.4.

A solução de elementos finitos $u_{h}$ satisfaz

\|(u-u_{h})^{\prime}\|_{L^{2}(I)}^{2}\leq C\sum_{i=1}^{n}\eta_{i}^{2}(u_{h}),

(1.98)

onde $\eta_{i}(u_{h})$ é chamado de elemento residual e é dado por

\eta_{i}(u_{h})=h_{i}\|f-u_{h}^{\prime\prime}\|_{L^{2}(I_{i})}.

(1.99)

Demonstração.

Tomando $e=u-u_{h}$ e usando a ortogonalidade de Galerkin (Teorema 1.2.1) temos

\|e^{\prime}\|_{L^{2}(I)}^{2}=\int_{I}e^{\prime}(e-\pi e)^{\prime}\,dx=\sum_{i% =1}^{n}\int_{I_{i}}e^{\prime}(e-\pi e)^{\prime}\,dx.

(1.100)

Então, aplicando integração por partes

\|e^{\prime}\|_{L^{2}(I)}^{2}=\sum_{i=1}^{n}\int_{I_{i}}(-e^{\prime\prime})(e-% \pi e)\,dx+[e^{\prime}(e-\pi e)]_{x_{i-1}}^{x_{i}}.

(1.101)

Daí, observando que $e-\pi e=0$ nos extremos dos intervalos $I_{i}$ e que $-e^{\prime\prime}=-(u-u_{h})^{\prime\prime}=-u^{\prime\prime}+u_{h}^{\prime% \prime}=f+u_{h}^{\prime\prime}$ , temos

\|e^{\prime}\|_{L^{2}(I)}^{2}=\sum_{i=1}^{n}\int_{I_{i}}(f+u_{h}^{\prime\prime% })(e-\pi e)\,dx.

(1.102)

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

$\displaystyle\\|e^{\prime}\\|_{L^{2}(I)}^{2}$	$\displaystyle\leq\sum_{i=1}^{n}\\|f+u_{h}\\|_{L^{2}(I_{i})}\\|e-\pi e\\|_{L^{2}(I_% {i})}\,dx$	(1.103)
	$\displaystyle\leq C\sum_{i=1}^{n}h_{i}\\|f+u_{h}\\|_{L^{2}(I_{i})}\\|e^{\prime}\\|% _{L^{2}(I_{i})}$	(1.104)
	$\displaystyle\leq C\left(\sum_{i=1}^{n}h_{i}^{2}\\|f+u_{h}\\|_{L^{2}(I_{i})}^{2}% \right)^{1/2}\left(\sum_{i=1}^{n}\\|e^{\prime}\\|_{L^{2}(I_{i})}^{2}\right)^{1/2}$	(1.105)
	$\displaystyle=C\left(\sum_{i=1}^{n}h_{i}^{2}\\|f+u_{h}\\|_{L^{2}(I_{i})}^{2}% \right)^{1/2}\\|e^{\prime}\\|_{L^{2}(I)},$	(1.106)

donde segue o resultado desejado. ∎

Observação 1.2.1.

No caso da solução de elementos finitos no espaço das funções lineares por partes, temos $u_{h}^{\prime\prime}=0$ . Logo, o elemento residual se resume em $\eta_{i}(u_{h})=h_{i}\|f\|_{L^{2}(I_{i})}$ .

1.2.5 Exercícios

Em revisão

E. 1.2.1.

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

	$\displaystyle-u^{\prime\prime}+u=2\operatorname{sen}x,\quad\forall x\in(-\pi,% \pi),$		(1.107)
	$\displaystyle u(-\pi)=u(\pi)=0.$		(1.108)

Código FENiCS.

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!

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.

$\displaystyle\\|e^{\prime}\\|_{L^{2}(I)}^{2}$	$\displaystyle\leq\sum_{i=1}^{n}\\|f+u_{h}\\|_{L^{2}(I_{i})}\\|e-\pi e\\|_{L^{2}(I_% {i})}\,dx$	(1.103)
	$\displaystyle\leq C\sum_{i=1}^{n}h_{i}\\|f+u_{h}\\|_{L^{2}(I_{i})}\\|e^{\prime}\\|% _{L^{2}(I_{i})}$	(1.104)
	$\displaystyle\leq C\left(\sum_{i=1}^{n}h_{i}^{2}\\|f+u_{h}\\|_{L^{2}(I_{i})}^{2}% \right)^{1/2}\left(\sum_{i=1}^{n}\\|e^{\prime}\\|_{L^{2}(I_{i})}^{2}\right)^{1/2}$	(1.105)
	$\displaystyle=C\left(\sum_{i=1}^{n}h_{i}^{2}\\|f+u_{h}\\|_{L^{2}(I_{i})}^{2}% \right)^{1/2}\\|e^{\prime}\\|_{L^{2}(I)},$	(1.106)

Política de uso de dados

Política de uso de dados

Método dos Elementos Finitos

1.2 Problema modelo

1.2.1 Formulação fraca

1.2.2 Formulação de elementos finitos

Exemplo 1.2.1.

1.2.3 Estimativa a priori

Teorema 1.2.1.

Demonstração.

Teorema 1.2.2.

Demonstração.

Teorema 1.2.3.

Demonstração.

Exemplo 1.2.2.

1.2.4 Estimativa a posteriori

Teorema 1.2.4.

Demonstração.

Observação 1.2.1.

1.2.5 Exercícios

E. 1.2.1.

Envie seu comentário