Método dos Elementos Finitos

1 Problemas unidimensionais 1.2 Problema modelo 1.4 Aplicação: EDP não-linear

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

1.3 Condições de contorno

Vamos estudar sobre soluções de elementos finitos para o problema modelo

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}-u^{\prime% \prime}=f,\quad x\in I=[0,L],

(1.121)

com diferentes condições de contorno.

1.3.1 Condições de Dirichlet

Consideramos o seguinte problema com condições de contorno de Dirichlet⁷⁷7Johann Peter Gustav Lejeune Dirichlet, 1805 - 1859, matemático alemão. Fonte: Wikipédia: Johann Peter Gustav Lejeune Dirichlet.: encontrar $u$ tal que

	$\displaystyle-u^{\prime\prime}=f,\quad\forall x\in I=(0,L),$		(1.122)
	$\displaystyle u(0)=u_{0},\quad u(L)=u_{L},$		(1.123)

com $u_{0},u_{L}\in\mathbb{R}$ e $f$ dados.

Tomando uma função teste $v\in\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}V_{0}:=H^% {1}_{0}(I):=\{v\in H^{1}(I);\leavevmode\nobreak\ v(0)=v(L)=0\}$ e multiplicando-a em (1.122), obtemos

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

(1.124)

Aplicando a integração por partes, temos

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

(1.125)

Desta forma, definimos o seguinte problema fraco associado: encontrar $u\in\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}V:=\{v\in H% ^{1}(I);\leavevmode\nobreak\ v(0)=u_{0},\leavevmode\nobreak\ v(L)=v_{L}\}$ tal que

\displaystyle\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}% a(u,v)=L(v),\quad\forall v\in V_{0},

(1.126)

onde $a(u,v)$ é a forma bilinear

a(u,v)=\int_{I}u^{\prime}v^{\prime}\,dx

(1.127)

e $L(v)$ é a forma linear

L(v)=\int_{I}fv\,dx.

(1.128)

Exemplo 1.3.1.

Consideramos o problema

	$\displaystyle-u^{\prime\prime}=1,\quad x\in I=[0,1],$		(1.129)
	$\displaystyle u(0)=1/2,\quad u(1)=1.$		(1.130)

Sua solução analítica é $u(x)=-x^{2}/2+x+1/2$ .

Para obtermos uma aproximação de elementos finitos, consideramos o seguinte problema fraco: encontrar $u\in V:=\{v\in H^{1}(I);\leavevmode\nobreak\ v(0)=1/2,\leavevmode\nobreak\ v(1% )=1\}$ tal que

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

(1.131)

para todo $v\in V_{0}=\{v\in H^{1}(I);\leavevmode\nobreak\ v(0)=v(1)=0\}$ , onde

	$\displaystyle a(u,v)=\int_{I}u^{\prime}v^{\prime}\,dx,$		(1.132)
	$\displaystyle L(v)=\int_{I}fv\,dx.$		(1.133)

Então, o problema de elementos finitos no espaço das funções afins por partes lê-se: encontrar $u_{h}\in V_{h}=\{v\in P_{1}(I);\leavevmode\nobreak\ v(0)=1/2,\leavevmode% \nobreak\ v(1)=1\}$ tal que

a(u_{h},v_{h})=L(v_{h}),

(1.134)

para todo $v_{h}\in V_{h,0}=\{v\in H^{1}(I);\leavevmode\nobreak\ v(0)=v(1)=0\}$ .

Este problema pode ser implementado com um código semelhante ao Código 4, mas com os devidos ajustes para as condições de contorno de Dirichlet não homogêneas. Estes ajustes estão presentes no código a seguir. Faça os ajustes e verifique!

Código 6: ex_mef1d_dirichlet.py

⬇

1# condição de contorno

2import numpy as np

3uD = fem.Function(V)

4# uD(0) = 0.5, uD(1) = 1

5def dirichlet_bc(x):

6 y = np.full(x.shape[1], 1.)

7 y[np.isclose(x[0], 0.)] = 0.5

8 return y

9uD.interpolate(dirichlet_bc)

1.3.2 Condições de Neumann

Consideramos o seguinte problema com condições de contorno de Neumann⁸⁸8Carl Gottfried Neumann, 1832 - 1925, matemático alemão. Fonte: Wikipédia: Carl Neumann. homogênea em $x=L$ : encontrar $u$ tal que

	$\displaystyle-u^{\prime\prime}=f,\quad\forall x\in I=[0,L],$		(1.135)
	$\displaystyle u(0)=u_{0},\quad u^{\prime}(L)=0,$		(1.136)

com $u_{0}$ e $f$ dados. Trata-se de um problema com condição de contorno de Dirichlet⁹⁹9Johann Peter Gustav Lejeune Dirichlet, 1805 - 1859, matemático alemão. Fonte: Wikipédia: Johann Peter Gustav Lejeune Dirichlet. à esquerda e condição de contorno de Neumann¹⁰¹⁰10Carl Gottfried Neumann, 1832 - 1925, matemático alemão. Fonte: Wikipédia: Carl Neumann. homogênea à direita.

Tomando uma função teste $v\in\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}V:=\{v\in H% ^{1}(I);\leavevmode\nobreak\ v(0)=0\}$ e multiplicando-a em (1.135), obtemos

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

(1.137)

Aplicando a integração por partes, temos

	$\displaystyle\int_{I}u^{\prime}v^{\prime}\,dx-\underbrace{u^{\prime}(L)v(L)}_{% u^{\prime}(L)=0}+\underbrace{u^{\prime}(0)v(0)}_{v(0)=0}$		(1.138)
	$\displaystyle\text{}\quad=\int_{I}fv\,dx.$		(1.139)

Desta forma, definimos o seguinte problema fraco associado: encontrar $u\in\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\tilde{V}% :=\{v\in H^{1}(I);\leavevmode\nobreak\ v(0)=u_{0}\}$ tal que

\displaystyle a(u,v)=L(v),\quad\forall\color[rgb]{0,0,1}\definecolor[named]{% pgfstrokecolor}{rgb}{0,0,1}v\in V,

(1.140)

onde $a(u,v)$ é a forma bilinear

a(u,v)=\int_{I}u^{\prime}v^{\prime}\,dx

(1.141)

e $L(v)$ é a forma linear

L(v)=\int_{I}fv\,dx.

(1.142)

Exemplo 1.3.2.

Consideramos o problema

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

Sua solução analítica é $u(x)=-x^{2}/2+x$ .

Podemos construir uma aproximação por elementos finitos do seguinte problema fraco associado: encontrar $u\in V=\{v\in H^{1}(I);\leavevmode\nobreak\ v(0)=0\}$ tal que

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

(1.145)

para todo $v\in V$ , com as formas bilinear $a(\cdot,\cdot)$ e linear $L(\cdot)$ dadas em (1.141) e (1.142).

Então, considerando elementos afins por partes, temos o seguinte problema de elementos finitos: encontrar $u_{h}\in V_{h}=\{v_{h}\in P_{1}(I);\leavevmode\nobreak\ v_{h}(0)=0\}$ tal que

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

(1.146)

A implementação deste problema de elementos finitos é semelhante a do Código 4, mas com os devidos ajustes para as condições de contorno. Estes ajustes estão presentes no código a seguir. Faça os ajustes e verifique!

Código 7: ex_mef1d_neumann.py

⬇

1import numpy as np

2uD = fem.Function(V)

3uD.interpolate(lambda x: np.zeros(x.shape[1]))

4# uD(0) = 0 (Dirichlet); x=1: Neumann homogênea (natural)

5boundary_dofs = fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0.))

6bc = fem.dirichletbc(uD, boundary_dofs)

Agora, consideramos o seguinte problema com condições de Neumann não-homogênea em $x=L$ : encontrar $u$ tal que

	$\displaystyle-u^{\prime\prime}=f,\quad\forall x\in I=(0,L),$		(1.147)
	$\displaystyle u(0)=u_{0},\quad u^{\prime}(L)=\alpha,$		(1.148)

com $u_{0}$ , $\alpha$ e $f$ dados.

Tomando uma função teste $v\in V:=\{v\in H^{1}(I);\leavevmode\nobreak\ v(0)=0\}$ e multiplicando-a em (1.147), obtemos

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

(1.149)

Aplicando a integração por partes, temos

\int_{I}u^{\prime}v^{\prime}\,dx-\alpha v(L)=\int_{I}fc\,dx.

(1.150)

Desta forma, definimos o seguinte problema fraco associado: encontrar $u\in\tilde{V}:=\{v\in H^{1}(I);\leavevmode\nobreak\ v(0)=u_{0}\}$ tal que

\displaystyle a(u,v)=L(v),\quad\forall v\in V,

(1.151)

onde $a(u,v)$ é a forma bilinear

a(u,v)=\int_{I}u^{\prime}v^{\prime}\,dx

(1.152)

e $L(v)$ é a forma linear

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}L(v)=\int_{I}% fv\,dx+\alpha v(L).

(1.153)

Exemplo 1.3.3.

Consideramos o problema

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

Sua solução analítica é $u(x)=-x^{2}/2+2x$ .

Agora, consideramos o seguinte problema fraco associado: encontrar $u\in V=\{v\in H^{1}(I);\leavevmode\nobreak\ v(0)=0\}$ tal que

a(u,v)=L(v),\quad\forall v\in V,

(1.156)

com

a(u,v)=\int_{I}u^{\prime}v^{\prime}\,dx

(1.157)

L(v)=\int_{I}fv\,dx+1\cdot v(1).

(1.158)

O problema de elementos finitos no espaço das funções lineares por partes lê-se: encontrar $u_{h}\in V_{h}=\{v_{h}\in P_{1}(I);\leavevmode\nobreak\ v_{h}(0)=0\}$ tal que

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

(1.159)

A implementação deste problema pode ser feita com um código semelhante ao Código 4, mas com os devidos ajustes para as condições de contorno e para a forma linear $L(\cdot)$ , que agora inclui um termo de contorno. Estes ajustes estão presentes no código a seguir. Faça os ajustes e verifique!

Código 8: ex_mef1d_neumann_nh.py

⬇

1# condição de contorno

2import numpy as np

3uD = fem.Function(V)

4uD.interpolate(lambda x: np.zeros(x.shape[1]))

5# uD(0) = 0 (Dirichlet); x=1: Neumann homogênea (natural)

6boundary_dofs = fem.locate_dofs_geometrical(V, lambda x: np.isclose(x[0], 0.))

7bc = fem.dirichletbc(uD, boundary_dofs)

9# problema mef

10import ufl

11from dolfinx import default_scalar_type

12from dolfinx.fem.petsc import LinearProblem

13u = ufl.TrialFunction(V)

14v = ufl.TestFunction(V)

16f = fem.Constant(domain, default_scalar_type(1.))

17g = fem.Constant(domain, default_scalar_type(1.))

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

20L = f * v * ufl.dx

21L += g * v * ufl.ds

1.3.3 Condições de Robin

Consideramos o seguinte problema com condições de contorno de Robin¹¹¹¹11Victor Gustave Robin, 1855 - 1897, matemático francês. Fonte: Wikipedia: Victor Gustave Robin.: encontrar $u$ tal que

	$\displaystyle-u^{\prime\prime}=f,\quad\forall x\in I=[0,L],$		(1.160)
	$\displaystyle u^{\prime}(0)=r_{0}(u(0)-s_{0}),\leavevmode\nobreak\ -u^{\prime}% (L)=r_{L}(u(L)-s_{L}),$		(1.161)

com $r_{0}$ , $r_{L}$ , $s_{0}$ , $s_{L}$ e $f$ dados.

Tomando uma função teste $v\in V=H^{1}(I)$ e multiplicando-a em (1.160), obtemos

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

(1.162)

Aplicando a integração por partes, temos

	$\displaystyle\int_{I}u^{\prime}v^{\prime}\,dx-\underbrace{u^{\prime}(L)v(L)}_{% -u^{\prime}(L)=r_{L}(u(L)-s_{L})}+\underbrace{u^{\prime}(0)v(0)}_{u^{\prime}(0% )=r_{0}(u(0)-s_{0})}$		(1.163)
	$\displaystyle\text{}\quad=\int_{I}fc\,dx.$		(1.164)

ou, mais adequadamente,

	$\displaystyle\int_{I}u^{\prime}v^{\prime}\,dx+r_{L}u(L)v(L)+r_{0}u(0)v(0)$		(1.165)
	$\displaystyle\text{}\quad=\int_{I}fv\,dx+r_{L}s_{L}v(L)+r_{0}s_{0}v(0).$		(1.166)

Desta forma, definimos o seguinte problema fraco associado: encontrar $u\in H^{1}(I)$ tal que

\displaystyle a(u,v)=L(v),\quad\forall v\in V,

(1.167)

onde $a(u,v)$ é a forma bilinear

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}a(u,v)=\int_{% I}u^{\prime}v^{\prime}\,dx+r_{L}u(L)v(L)+r_{0}u(0)v(0)

(1.168)

e $L(v)$ é a forma linear

\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}L(v)=\int_{I}% fv\,dx+r_{L}s_{L}v(L)+r_{0}s_{0}v(0).

(1.169)

Exemplo 1.3.4.

Consideramos o problema

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

Sua solução analítica é $u(x)=-x^{2}/2+5x/6+5/6$ .

Aqui, tomamos o seguinte problema fraco: encontrar $u\in V=H^{1}(I)$ tal que

a(u,v)=L(v),\quad\forall v\in V,

(1.172)

onde

a(u,v)=\int_{I}u^{\prime}v^{\prime}\,dx+u(1)v(1)+u(0)v(0)

(1.173)

L(v)=\int_{I}fv\,dx+1\cdot v(1).

(1.174)

Então, uma aproximação por elementos finitos lineares por partes pode ser obtida resolvendo o seguinte problema: encontrar $u_{h}\in V_{h}=P_{1}(I)$ tal que

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

(1.175)

A implementação deste problema de elementos finitos é semelhante a do Código 4, mas com os devidos ajustes para as condições de contorno e para as formas bilinear e linear. Para tanto, precisamos marcar as bordas do domínio para podermos integrar os termos de contorno. Estes ajustes estão presentes no código a seguir. Faça os ajustes e verifique!

Código 9: ex_mef1d_robin.py

⬇

1# boundary colors

2import numpy as np

3from dolfinx.mesh import locate_entities

4from dolfinx.mesh import meshtags

5boundaries = [(0, lambda x: np.isclose(x[0], 0.)),

6 (1, lambda x: np.isclose(x[0], 1.))]

7facet_indices, facet_markers = [], []

8fdim = domain.topology.dim - 1

9for (marker, locator) in boundaries:

10 facets = locate_entities(domain, fdim, locator)

11 facet_indices.append(facets)

12 facet_markers.append(np.full_like(facets, marker))

13facet_indices = np.hstack(facet_indices).astype(np.int32)

14facet_markers = np.hstack(facet_markers).astype(np.int32)

15sorted_facets = np.argsort(facet_indices)

16facet_tag = meshtags(domain, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])

18# problema mef

19import ufl

20from dolfinx import default_scalar_type

21from dolfinx.fem.petsc import LinearProblem

22u = ufl.TrialFunction(V)

23v = ufl.TestFunction(V)

25f = fem.Constant(domain, default_scalar_type(1.))

26g = fem.Constant(domain, default_scalar_type(1.))

28ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag)

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

31a += u * v * ds(1) + u * v * ds(0)

32L = f * v * ufl.dx

33L += g * v * ds(1)