Colabore! Saiba mais screen_rotation

Ao navegar por este site , você concorda com a Política de Uso de Dados.

| | | | |

1 Problemas Unidimensionais 1.4 Malhas Auto-Adaptativas 1.6 Aplicação: EDP de Advecção-Difusão

1.5 Aplicação: EDP Evolutiva

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

Em revisão

Como exemplo de aplicação do método de elementos finitos (MEF) na solução de equações diferenciais parciais evolutivas no tempo, consideramos a equação do calor com dadas condição inicial e condições de contorno de Dirichlet homogêneas

	$\displaystyle u_{t}=\alpha u_{xx}+f,\leavevmode\nobreak\ (t,x)\in(0,t_{f}]% \times(a,b),$		(1.165a)
	$\displaystyle u(0,x)=u_{0}(x),x\in[a,b],$		(1.165b)
	$\displaystyle u(t,a)=u(t,b)=0,\leavevmode\nobreak\ t\in[0,t_{f}],$		(1.165c)

onde $f=f(t,x)$ denota uma dada fonte.

1.5.1 Discretização do Tempo

Consideramos os $n_{t}+1$ tempos discretos $t^{(k)}=kh_{t}$ , passo no tempo $h_{t}=t_{f}/n_{t}$ , $k=0,1,2,\dotsc,n_{t}$ . Seguindo esquema $\theta$ denotando $u^{(k)}\approx u\left(t^{(k)},x\right)$ e $f^{(k)}=f\left(t^{(k)},x\right)$ , o problema (1.5) pode ser aproximado pela iteração

	$\displaystyle\frac{u^{(k+1)}-u^{(k)}}{h_{t}}=\theta\left(\alpha u^{(k+1)}_{xx}% +f^{(k+1)}\right)$
	$\displaystyle\qquad\qquad\qquad\;\;(1-\theta)\left(\alpha u^{(k)}_{xx}+f^{(k)}% \right),$		(1.166a)
	$\displaystyle u^{(k+1)}(a)=u^{(k+1)}(b)=0,$		(1.166b)

onde $u^{(0)}=u_{0}$ .

Observação 1.5.1.

(Esquema $\theta$ .) O esquema $\theta$ e um forma robusta de escrever diferentes esquemas de discretização em uma única expressão:

•

$\theta=0.$ : Euler explícito.
•

$\theta=1.$ : Euler implícito.
•

$\theta=0.5$ : Crank-Nicolson.

Por simplificação da notação, vamos suprimir o super-índice $k$ , denotando $u^{(k+1)}:=u$ , $u^{(k)}=u^{0}$ e similar para $f^{(k)}$ . Com isso e rearranjando os termos, cada iteração (1.5.1) se resume ao seguinte problema de valores de contorno

	$\displaystyle-\alpha\theta u_{xx}+\frac{1}{h_{t}}u=\frac{1}{h_{t}}u^{0}+(1-% \theta)\alpha u^{0}_{xx}$
	$\displaystyle\qquad\qquad\qquad\quad+(1-\theta)f^{0}+\theta f,$		(1.167a)
	$\displaystyle u(a)=u(b)=0.$		(1.167b)

1.5.2 Formulação de Elementos Finitos

A formulação fraca do problema (1.5.1) consiste em: encontrar $u\in V:=H^{1}_{0}(a,b)$ tal que

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

(1.168)

onde

	$\displaystyle a(u,v):=\int_{a}^{b}\theta\alpha u^{\prime}v^{\prime}\,dx+\int_{% a}^{b}\frac{1}{h_{t}}uv\,dx,$		(1.169)
	$\displaystyle L(v):=(1-\theta)\int_{a}^{b}\alpha u^{0^{\prime}}v^{\prime}\,dx+% \int_{a}^{b}\frac{1}{h_{t}}u^{0}v\,dx$
	$\displaystyle\qquad\theta\int_{a}^{b}fv\,dx+(1-\theta)\int_{a}^{b}f^{0}v\,dx$		(1.170)

Então, assumindo uma malha com $n_{x}$ células $I_{i}=[x_{i},x_{i+1}]$ de tamanho $h_{x}=(b-a)/n_{x}$ e nodos $x_{i}=a+(i-1)h_{x}$ , $i=0,1,2,\dotsc,n_{x}$ , escolhemos o espaço de elementos finitos

V_{h,0}:=\{v\in C^{0}([a,b]):\leavevmode\nobreak\ v|_{I_{i}}\in P_{1}(I_{i}),% \leavevmode\nobreak\ i=0,1,\dotsc,n_{x},v(a)=v(b)=0\}.

(1.171)

Com isso, a formulação de elementos finitos do problema (1.5.1) consiste em: encontrar $u_{h}\in V_{h,0}$ tal que

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

(1.172)

Exemplo 1.5.1.

Consideramos o seguinte problema de calor

	$\displaystyle u_{t}=u_{xx}+(\pi^{2}-1)e^{-t}\operatorname{sen}(\pi x),% \leavevmode\nobreak\ (t,x)\in(0,1]\times(0,1),$		(1.173a)
	$\displaystyle u(0,x)=\operatorname{sen}(\pi x),\leavevmode\nobreak\ x\in[0,1],$		(1.173b)
	$\displaystyle u(t,0)=u(t,1)=0.$		(1.173c)

⬇

1from mpi4py import MPI

2import ufl

3from dolfinx import mesh

4from dolfinx import fem

5from dolfinx import default_scalar_type

6from dolfinx.fem.petsc import LinearProblem

8# parâmetros

9tf = 1.

10alpha = 1.

12# esquema theta

13theta = 0.5

15# discretização no tempo

16nt = 10

17ht = tf/nt

19# malha

20domain = mesh.create_unit_interval(MPI.COMM_WORLD,

21 nx = 5)

22x = ufl.SpatialCoordinate(domain)

24# espaço

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

27# fonte

28f = fem.Function(V)

29def f(t,x):

30 return (ufl.pi**2-1.)*ufl.exp(-t)*ufl.sin(ufl.pi*x[0])

32# condição de contorno

33import numpy as np

34uD = fem.Function(V)

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

37def boundary_D(x):

38 return np.logical_or(np.isclose(x[0], 0.),

39 np.isclose(x[0], 1.))

41dofs_D = fem.locate_dofs_geometrical(V, boundary_D)

42bc = fem.dirichletbc(uD, dofs_D)

44# mef fun.s

45u = ufl.TrialFunction(V)

46v = ufl.TestFunction(V)

48# condição inicial

49t = 0.

50u0 = fem.Function(V)

51u0.interpolate(lambda x: np.sin(np.pi*x[0]))

53# fonte

54def f(t, x):

55 return (ufl.pi**2-1.)*ufl.exp(-t)*ufl.sin(ufl.pi*x[0])

57# visualização (paraview)

58from dolfinx import io

59from pathlib import Path

60results_folder = Path("results")

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

63# iteração no tempo

64for k in range(nt):

65 t += ht

67 # forma bilinear

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

69 a += u * v / ht * ufl.dx

71 # forma linear

72 L = (theta-1.) * ufl.dot(ufl.grad(u0), ufl.grad(v)) * ufl.dx

73 L += u0 * v / ht * ufl.dx

74 L += theta * f(t, x) * v * ufl.dx

75 L += (1.-theta) * f(t-ht, x) * v * ufl.dx

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

78 uh = problem.solve()

80 # armazena para visualização (paraview)

81 filename = results_folder / f"u_{k:0>6}"

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

83 xdmf.write_mesh(domain)

84 xdmf.write_function(uh, t)

86 u0.x.array[:] = uh.x.array[:]

1.5.3 Exercícios

Em construção

Envie seu comentário

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. Aproveito para agradecer a todas/os que de forma assídua ou esporádica contribuem enviando correções, sugestões e críticas!

| | | |