Política de uso de dados

Ao navegar por este site, você concorda com a política de uso de dados.

Política de uso de dados

Ao navegar por este site, você concorda com a política de uso de dados.

| | | |

Método de Elementos Finitos

1 Problemas Unidimensionais 1.6 Aplicação: EDP de Advecção-Difusão 1.8 Seleção de Aplicações

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

1.7 Aplicação: EDP Não-Linear

Em construção

Como exemplo de aplicação do MEF na solução de equações diferenciais parciais não-lineares, consideramos a equação de Fisher¹⁰¹⁰endnote: ¹⁰Ronald Aylmer Fisher, 1890-1962, biólogo inglês. Fonte: Wikipédia: Ronald Fisher. com dadas condição inicial e condições de contorno de Neumann¹¹¹¹endnote: ¹¹Carl Gottfried Neumann, 1832 - 1925, matemático alemão. Fonte: Wikipédia: Carl Neumann.

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

1.7.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)$ , o problema (1.7) pode ser aproximado pela iteração

	$\displaystyle\frac{u^{(k+1)}-u^{(k)}}{h_{t}}=\theta\left[u^{(k+1)}_{xx}+u^{(k+% 1)}\left(1-u^{(k+1})\right)\right]$
	$\displaystyle\qquad\qquad\qquad\;\;(1-\theta)\left[u^{(k)}_{xx}++u^{(k)}\left(% 1-u^{(k)}\right)\right],$		(1.183a)
	$\displaystyle u_{x}^{(k+1)}(0)=u_{x}^{(k+1)}(1)=0,$		(1.183b)

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

Observação 1.7.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}$ . Com isso e rearranjando os termos, cada iteração (1.7.1) se resume ao seguinte problema de valores de contorno

	$\displaystyle\frac{1}{h_{t}}u-\frac{1}{h_{t}}u^{0}-\theta\left[u_{x}x+u(1-u)\right]$
	$\displaystyle\qquad\qquad-(1-\theta)\left[u^{0}_{x}x+u^{0}(1-u^{0})\right],$		(1.184a)
	$\displaystyle u_{x}(0)=u_{x}(1)=0.$		(1.184b)

1.7.2 Formulação de Elementos Finitos

Em revisão

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

F(u;v)=0,\leavevmode\nobreak\ \forall v\in V,

(1.185)

onde

		$\displaystyle F(u;v):=\int_{0}^{1}\frac{1}{h_{t}}u\,dx-\int_{0}^{1}\frac{1}{h_% {t}}u^{0}\,dx$		(1.186)
		$\displaystyle\qquad+\theta\int_{0}^{1}u_{x}v_{x}\,dx-\theta\int_{0}^{1}u(1-u)v% \,dx$
		$\displaystyle\qquad+(1-\theta)\int_{0}^{1}u^{0}_{x}v_{x}\,dx-(1-\theta)\int_{0% }^{1}u^{0}(1-u^{0})v\,dx.$

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

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

(1.187)

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

F(u_{h};v)=0,\leavevmode\nobreak\ \forall v_{h}\in V_{h}.

(1.188)

Observação 1.7.2.

O problema (1.188) consiste em um sistema de equações não-lineares.

Exemplo 1.7.1.

Consideramos a equação de Fisher com condições inicial e de contorno

	$\displaystyle u_{t}=u_{xx}+u(1-u),\leavevmode\nobreak\ t\in(0,t_{f})\times(0,1),$		(1.189a)
	$\displaystyle u(0,x)=\cos^{2}(\pi x),\leavevmode\nobreak\ x\in[0,1],$		(1.189b)
	$\displaystyle u_{x}(t,0)=u_{x}(t,1)=0,\leavevmode\nobreak\ t\in[0,t_{f}],$		(1.189c)

com $tf=5$ .

Código 7: ex_mef1d_fisher.py

⬇

1 from mpi4py import MPI

2 import numpy as np

3 import ufl

4 from dolfinx import mesh

5 from dolfinx import fem

6 from dolfinx import default_scalar_type

7 from dolfinx.fem.petsc import NonlinearProblem

8 from dolfinx.nls.petsc import NewtonSolver

10 # parâmetros

11 tf = 5.

13 # esquema theta

14 theta = 0.5

16 # discretização no tempo

17 nt = 100

18 ht = tf/nt

20 # malha

21 domain = mesh.create_unit_interval(MPI.COMM_WORLD,

22 nx = 5)

23 x = ufl.SpatialCoordinate(domain)

25 # espaço

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

28 # mef fun.s

29 v = ufl.TestFunction(V)

30 u = fem.Function(V)

32 # condição inicial

33 t = 0.

34 u0 = fem.Function(V)

35 u0.interpolate(lambda x: np.cos(np.pi*x[0])**2)

37 # inicialização

38 u.x.array[:] = u0.x.array[:]

40 # visualização (paraview)

41 from dolfinx import io

42 from pathlib import Path

43 results_folder = Path("results")

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

46 # armazena para visualização (paraview)

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

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

49 xdmf.write_mesh(domain)

50 xdmf.write_function(u, 0.)

53 # iteração no tempo

54 for k in range(nt):

56 t += ht

57 print(f"{k+1}: t = {t:.4g}")

59 # forma fraca

60 ## time term

61 F = 1./ht * u * v * ufl.dx

62 F -= 1./ht * u0 * v * ufl.dx

63 ## diffusion term

64 F += theta * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx

65 F += (1.-theta) * ufl.dot(ufl.grad(u0), ufl.grad(v)) * ufl.dx

66 ## reaction term

67 F -= theta * u * (1. - u) * v * ufl.dx

68 F -= (1.-theta) * u0 * (1. - u0) * v * ufl.dx

70 problem = NonlinearProblem(F, u)

71 solver = NewtonSolver(MPI.COMM_WORLD, problem)

72 n, converged = solver.solve(u)

73 print(f"\tNewton iterations: {n}")

74 assert(converged)

76 # armazena para visualização (paraview)

77 filename = results_folder / f"u_{k+1:0>6}"

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

79 xdmf.write_mesh(domain)

80 xdmf.write_function(u, t)

82 u0.x.array[:] = u.x.array[:]

1.7.3 Exercícios

Em construção

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.

Método de Elementos Finitos

1 Problemas Unidimensionais 1.6 Aplicação: EDP de Advecção-Difusão 1.8 Seleção de Aplicações

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

1.7 Aplicação: EDP Não-Linear

Em construção

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

1.7.1 Discretização do Tempo

	$\displaystyle\frac{u^{(k+1)}-u^{(k)}}{h_{t}}=\theta\left[u^{(k+1)}_{xx}+u^{(k+% 1)}\left(1-u^{(k+1})\right)\right]$
	$\displaystyle\qquad\qquad\qquad\;\;(1-\theta)\left[u^{(k)}_{xx}++u^{(k)}\left(% 1-u^{(k)}\right)\right],$		(1.183a)
	$\displaystyle u_{x}^{(k+1)}(0)=u_{x}^{(k+1)}(1)=0,$		(1.183b)

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

Observação 1.7.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.

	$\displaystyle\frac{1}{h_{t}}u-\frac{1}{h_{t}}u^{0}-\theta\left[u_{x}x+u(1-u)\right]$
	$\displaystyle\qquad\qquad-(1-\theta)\left[u^{0}_{x}x+u^{0}(1-u^{0})\right],$		(1.184a)
	$\displaystyle u_{x}(0)=u_{x}(1)=0.$		(1.184b)

1.7.2 Formulação de Elementos Finitos

Em revisão

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

F(u;v)=0,\leavevmode\nobreak\ \forall v\in V,

(1.185)

onde

		$\displaystyle F(u;v):=\int_{0}^{1}\frac{1}{h_{t}}u\,dx-\int_{0}^{1}\frac{1}{h_% {t}}u^{0}\,dx$		(1.186)
		$\displaystyle\qquad+\theta\int_{0}^{1}u_{x}v_{x}\,dx-\theta\int_{0}^{1}u(1-u)v% \,dx$
		$\displaystyle\qquad+(1-\theta)\int_{0}^{1}u^{0}_{x}v_{x}\,dx-(1-\theta)\int_{0% }^{1}u^{0}(1-u^{0})v\,dx.$

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

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

(1.187)

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

F(u_{h};v)=0,\leavevmode\nobreak\ \forall v_{h}\in V_{h}.

(1.188)

Observação 1.7.2.

O problema (1.188) consiste em um sistema de equações não-lineares.

Exemplo 1.7.1.

Consideramos a equação de Fisher com condições inicial e de contorno

	$\displaystyle u_{t}=u_{xx}+u(1-u),\leavevmode\nobreak\ t\in(0,t_{f})\times(0,1),$		(1.189a)
	$\displaystyle u(0,x)=\cos^{2}(\pi x),\leavevmode\nobreak\ x\in[0,1],$		(1.189b)
	$\displaystyle u_{x}(t,0)=u_{x}(t,1)=0,\leavevmode\nobreak\ t\in[0,t_{f}],$		(1.189c)

com $tf=5$ .

Código 7: ex_mef1d_fisher.py

⬇

1 from mpi4py import MPI

2 import numpy as np

3 import ufl

4 from dolfinx import mesh

5 from dolfinx import fem

6 from dolfinx import default_scalar_type

7 from dolfinx.fem.petsc import NonlinearProblem

8 from dolfinx.nls.petsc import NewtonSolver

10 # parâmetros

11 tf = 5.

13 # esquema theta

14 theta = 0.5

16 # discretização no tempo

17 nt = 100

18 ht = tf/nt

20 # malha

21 domain = mesh.create_unit_interval(MPI.COMM_WORLD,

22 nx = 5)

23 x = ufl.SpatialCoordinate(domain)

25 # espaço

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

28 # mef fun.s

29 v = ufl.TestFunction(V)

30 u = fem.Function(V)

32 # condição inicial

33 t = 0.

34 u0 = fem.Function(V)

35 u0.interpolate(lambda x: np.cos(np.pi*x[0])**2)

37 # inicialização

38 u.x.array[:] = u0.x.array[:]

40 # visualização (paraview)

41 from dolfinx import io

42 from pathlib import Path

43 results_folder = Path("results")

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

46 # armazena para visualização (paraview)

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

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

49 xdmf.write_mesh(domain)

50 xdmf.write_function(u, 0.)

53 # iteração no tempo

54 for k in range(nt):

56 t += ht

57 print(f"{k+1}: t = {t:.4g}")

59 # forma fraca

60 ## time term

61 F = 1./ht * u * v * ufl.dx

62 F -= 1./ht * u0 * v * ufl.dx

63 ## diffusion term

64 F += theta * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx

65 F += (1.-theta) * ufl.dot(ufl.grad(u0), ufl.grad(v)) * ufl.dx

66 ## reaction term

67 F -= theta * u * (1. - u) * v * ufl.dx

68 F -= (1.-theta) * u0 * (1. - u0) * v * ufl.dx

70 problem = NonlinearProblem(F, u)

71 solver = NewtonSolver(MPI.COMM_WORLD, problem)

72 n, converged = solver.solve(u)

73 print(f"\tNewton iterations: {n}")

74 assert(converged)

76 # armazena para visualização (paraview)

77 filename = results_folder / f"u_{k+1:0>6}"

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

79 xdmf.write_mesh(domain)

80 xdmf.write_function(u, t)

82 u0.x.array[:] = u.x.array[:]

1.7.3 Exercícios

Em construção

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.

Pedro H A Konzen

| | | |