Método dos Elementos Finitos

1 Problemas Unidimensionais 1.2 Problema Modelo 1.4 Malhas Auto-Adaptativas

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

1.3 Condições de Contorno

Em revisão

Nesta seção, vamos discutir sobre soluções de elementos finitos para a equações diferencial

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

(1.105)

com diferentes condições de contorno.

1.3.1 Condições de Dirichlet

Em revisão

Consideramos o seguinte problema com condições de contorno de Dirichlet⁵⁵5Johann 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.106)
	$\displaystyle u(0)=u_{0},\quad u(L)=u_{L},$		(1.107)

com $u_{0}$ , $u_{L}$ 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.106), obtemos

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

(1.108)

Aplicando a integração por partes, temos

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

(1.109)

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

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

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

(1.111)

e $L(v)$ é a forma linear

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

(1.112)

Exemplo 1.3.1.

Consideramos o problema

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

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

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.116)
	$\displaystyle L(v)=\int_{I}fv\,dx.$		(1.117)

Então, o problema de elementos finitos no espaço das funções lineares 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.118)

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

Código 4: ex_mef1d_dirichlet.py

⬇

1from mpi4py import MPI

3# malha

4from dolfinx import mesh

5domain = mesh.create_unit_interval(MPI.COMM_WORLD,

6 nx = 5)

7# espaço

8from dolfinx import fem

9V = fem.functionspace(domain, ('P', 1))

11# condição de contorno

12import numpy as np

13uD = fem.Function(V)

14def dirichlet_bc(x):

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

16 y[x[0,:] > 0.5] = 1.

17 return y

18uD.interpolate(dirichlet_bc)

20tdim = domain.topology.dim

21fdim = tdim - 1

22domain.topology.create_connectivity(fdim, tdim)

23boundary_facets = mesh.exterior_facet_indices(domain.topology)

24boundary_dofs = fem.locate_dofs_topological(V, fdim,

25 boundary_facets)

26bc = fem.dirichletbc(uD, boundary_dofs)

28# problema mef

29import ufl

30from dolfinx import default_scalar_type

31from dolfinx.fem.petsc import LinearProblem

32u = ufl.TrialFunction(V)

33v = ufl.TestFunction(V)

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

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

38L = f * v * ufl.dx

40problem = LinearProblem(a, L, bcs=[bc])

41uh = problem.solve()

43# armazena para visualização (paraview)

44from dolfinx import io

45from pathlib import Path

46results_folder = Path("results")

47results_folder.mkdir(exist_ok=True, parents=True)

48filename = results_folder / "u"

49with io.XDMFFile(domain.comm, filename.with_suffix(".xdmf"), "w") as xdmf:

50 xdmf.write_mesh(domain)

51 xdmf.write_function(uh)

1.3.2 Condições de Neumann

Em revisão

Consideramos o seguinte problema com condições de contorno de Neumann⁶⁶6Carl 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.119)
	$\displaystyle u(0)=u_{0},\quad u^{\prime}(L)=0,$		(1.120)

com $u_{0}$ e $f$ dados. Trata-se de um problema com condição de contorno de Dirichlet à esquerda e condição de contorno de Neumann⁷⁷7Carl Gottfried Neumann, 1832 - 1925, matemático alemão. Fonte: Wikipédia: Carl Neumann. homogênea à direita.

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

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

(1.121)

Aplicando a integração por partes, temos

\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}=\int_{I}fv\,dx.

(1.122)

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

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

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

(1.124)

e $L(v)$ é a forma linear

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

(1.125)

Exemplo 1.3.2.

Consideramos o problema

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

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

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

Então, considerando elementos lineares 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.129)

Código 5: ex_mef1d_neumann.py

⬇

1from mpi4py import MPI

3# malha

4from dolfinx import mesh

5domain = mesh.create_unit_interval(MPI.COMM_WORLD,

6 nx = 5)

7# espaço

8from dolfinx import fem

9V = fem.functionspace(domain, ('P', 1))

11# c.c. dirichlet

12import numpy as np

13from dolfinx.fem import dirichletbc, locate_dofs_geometrical

14uD = fem.Function(V)

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

17def boundary_D(x):

18 return np.isclose(x[0], 0.)

20dofs_D = locate_dofs_geometrical(V, boundary_D)

21bc = dirichletbc(uD, dofs_D)

23# problema mef

24import ufl

25from dolfinx import default_scalar_type

26from dolfinx.fem.petsc import LinearProblem

27u = ufl.TrialFunction(V)

28v = ufl.TestFunction(V)

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

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

33L = f * v * ufl.dx

35problem = LinearProblem(a, L, bcs=[bc])

36uh = problem.solve()

38# armazena para visualização (paraview)

39from dolfinx import io

40from pathlib import Path

41results_folder = Path("results")

42results_folder.mkdir(exist_ok=True, parents=True)

43filename = results_folder / "u"

44with io.XDMFFile(domain.comm, filename.with_suffix(".xdmf"), "w") as xdmf:

45 xdmf.write_mesh(domain)

46 xdmf.write_function(uh)

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.130)
	$\displaystyle u(0)=u_{0},\quad u^{\prime}(L)=\alpha,$		(1.131)

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.130), obtemos

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

(1.132)

Aplicando a integração por partes, temos

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

(1.133)

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

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

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

(1.135)

e $L(v)$ é a forma linear

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

(1.136)

Exemplo 1.3.3.

Consideramos o problema

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

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

com

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

(1.140)

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

(1.141)

Então, consideramos o seguinte problema de elementos finitos associado: 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.142)

Código 6: ex_mef1d_neumann_nh.py

⬇

1from mpi4py import MPI

3# malha

4from dolfinx import mesh

5domain = mesh.create_unit_interval(MPI.COMM_WORLD,

6 nx = 5)

7# espaço

8from dolfinx import fem

9V = fem.functionspace(domain, ('P', 1))

11# c.c. dirichlet

12import numpy as np

13from dolfinx.fem import dirichletbc, locate_dofs_geometrical

14uD = fem.Function(V)

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

17def boundary_D(x):

18 return np.isclose(x[0], 0.)

20dofs_D = locate_dofs_geometrical(V, boundary_D)

21bc = dirichletbc(uD, dofs_D)

23# problema mef

24import ufl

25from dolfinx import default_scalar_type

26from dolfinx.fem.petsc import LinearProblem

27u = ufl.TrialFunction(V)

28v = ufl.TestFunction(V)

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

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

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

34L = f * v * ufl.dx

35L += g * v * ufl.ds

37problem = LinearProblem(a, L, bcs=[bc])

38uh = problem.solve()

40# armazena para visualização (paraview)

41from dolfinx import io

42from pathlib import Path

43results_folder = Path("results")

44results_folder.mkdir(exist_ok=True, parents=True)

45filename = results_folder / "u"

46with io.XDMFFile(domain.comm, filename.with_suffix(".xdmf"), "w") as xdmf:

47 xdmf.write_mesh(domain)

48 xdmf.write_function(uh)

1.3.3 Condições de Robin

Em revisão

Consideramos o seguinte problema com condições de contorno de Robin⁸⁸8Victor 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.143)
	$\displaystyle u^{\prime}(0)=r_{0}(u(0)-s_{0}),\leavevmode\nobreak\ -u^{\prime}% (L)=r_{L}(u(L)-s_{L}),$		(1.144)

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.143), obtemos

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

(1.145)

Aplicando a integração por partes, temos

\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})}=\int_{I}fc\,dx.

(1.146)

ou, mais adequadamente,

\int_{I}u^{\prime}v^{\prime}\,dx+r_{L}u(L)v(L)+r_{0}u(0)v(0)=\int_{I}fc\,dx+r_% {L}s_{L}v(L)+r_{0}s_{0}v(0).

(1.147)

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

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

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

(1.149)

e $L(v)$ é a forma linear

L(v)=\int_{I}fv\,dx+r_{L}s_{L}v(L)+r_{0}s_{0}v(0).

(1.150)

Exemplo 1.3.4.

Consideramos o problema

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

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

onde

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

(1.154)

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

(1.155)

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

⬇

1from mpi4py import MPI

3# malha

4from dolfinx import mesh

5domain = mesh.create_unit_interval(MPI.COMM_WORLD,

6 nx = 5)

7# espaço

8from dolfinx import fem

9V = fem.functionspace(domain, ('P', 1))

11# boundary colors

12from dolfinx.mesh import locate_entities

13from dolfinx.mesh import meshtags

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

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

16facet_indices, facet_markers = [], []

17fdim = domain.topology.dim - 1

18for (marker, locator) in boundaries:

19 facets = locate_entities(domain, fdim, locator)

20 facet_indices.append(facets)

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

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

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

24sorted_facets = np.argsort(facet_indices)

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

27# problema mef

28import ufl

29from dolfinx import default_scalar_type

30from dolfinx.fem.petsc import LinearProblem

31u = ufl.TrialFunction(V)

32v = ufl.TestFunction(V)

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

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

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

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

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

41L = f * v * ufl.dx

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

44problem = LinearProblem(a, L, bcs=[])

45uh = problem.solve()

47# armazena para visualização (paraview)

48from dolfinx import io

49from pathlib import Path

50results_folder = Path("results")

51results_folder.mkdir(exist_ok=True, parents=True)

52filename = results_folder / "u"

53with io.XDMFFile(domain.comm, filename.with_suffix(".xdmf"), "w") as xdmf:

54 xdmf.write_mesh(domain)

55 xdmf.write_function(uh)

1.3.4 Exercícios

Em revisão

E. 1.3.1.

Considere o problema

	$\displaystyle-u^{\prime\prime}+u^{\prime}+2u=-\cos(x),\quad x\in(0,\pi/2),$		(1.157)
	$\displaystyle u(0)=-0,3,\quad u(\pi/2)=-0,1.$		(1.158)

Obtenha uma aproximação por elementos finitos para a solução deste problema, empregando o espaço de elementos finitos linear sobre uma malha uniforme com $10$ células. Então, compare a aproximação computada com sua solução analítica $u(x)=0,1(\operatorname{sen}(x)+3\cos(x))$ , bem como, compute o erro $\|u-u_{h}\|_{L^{2}}$ .

Código.

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.

Política de uso de dados

Política de uso de dados

Método dos Elementos Finitos

1.3 Condições de Contorno

1.3.1 Condições de Dirichlet

Exemplo 1.3.1.

1.3.2 Condições de Neumann

Exemplo 1.3.2.

Exemplo 1.3.3.

1.3.3 Condições de Robin

Exemplo 1.3.4.

1.3.4 Exercícios

E. 1.3.1.

Envie seu comentário