| | | | |

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

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

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 Fisher1010endnote: 10Ronald Aylmer Fisher, 1890-1962, biólogo inglês. Fonte: Wikipédia: Ronald Fisher. com dadas condição inicial e condições de contorno de Neumann1111endnote: 11Carl Gottfried Neumann, 1832 - 1925, matemático alemão. Fonte: Wikipédia: Carl Neumann.

ut=uxx+u(1u),(t,x)(0,tf]×(0,1), (1.182a)
u(0,x)=u0(x),x[0,1], (1.182b)
ux(t,0)=ux(t,1)=0,t[0,tf]. (1.182c)

1.7.1 Discretização do Tempo

Consideramos os nt+1 tempos discretos t(k)=kht, passo no tempo ht=tf/nt, k=0,1,2,,nt. Seguindo esquema θ denotando u(k)u(t(k),x), o problema (1.7) pode ser aproximado pela iteração

u(k+1)u(k)ht=θ[uxx(k+1)+u(k+1)(1u(k+1))]
(1θ)[uxx(k)++u(k)(1u(k))], (1.183a)
ux(k+1)(0)=ux(k+1)(1)=0, (1.183b)

onde u(0)=u0.

Observação 1.7.1.

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

  • θ=0.: Euler explícito.

  • θ=1.: Euler implícito.

  • θ=0.5: Crank-Nicolson.

Por simplificação da notação, vamos suprimir o super-índice k, denotando u(k+1):=u, u(k)=u0. Com isso e rearranjando os termos, cada iteração (1.7.1) se resume ao seguinte problema de valores de contorno

1htu1htu0θ[uxx+u(1u)]
(1θ)[ux0x+u0(1u0)], (1.184a)
ux(0)=ux(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 uV:=H1[0,1] tal que

F(u;v)=0,vV, (1.185)

onde

F(u;v):=011htu𝑑x011htu0𝑑x (1.186)
+θ01uxvx𝑑xθ01u(1u)v𝑑x
+(1θ)01ux0vx𝑑x(1θ)01u0(1u0)v𝑑x.

Então, assumindo uma malha com nx células Ii=[xi,xi+1] de tamanho hx=1/nx e nodos xi=(i1)hx, i=0,1,2,,nx, escolhemos o espaço de elementos finitos

Vh:={vC0([a,b]):v|IiP1(Ii),i=0,1,,nx}. (1.187)

Com isso, a formulação de elementos finitos do problema (1.7.1) consiste em: encontrar uhVh tal que

F(uh;v)=0,vhVh. (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

ut=uxx+u(1u),t(0,tf)×(0,1), (1.189a)
u(0,x)=cos2(πx),x[0,1], (1.189b)
ux(t,0)=ux(t,1)=0,t[0,tf], (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
9
10  # parâmetros
11  tf = 5.
12
13  # esquema theta
14  theta = 0.5
15
16  # discretização no tempo
17  nt = 100
18  ht = tf/nt
19
20  # malha
21  domain = mesh.create_unit_interval(MPI.COMM_WORLD,
22                                     nx = 5)
23  x = ufl.SpatialCoordinate(domain)
24
25  # espaço
26  V = fem.functionspace(domain, ('P', 1))
27
28  # mef fun.s
29  v = ufl.TestFunction(V)
30  u = fem.Function(V)
31
32  # condição inicial
33  t = 0.
34  u0 = fem.Function(V)
35  u0.interpolate(lambda x: np.cos(np.pi*x[0])**2)
36
37  # inicialização
38  u.x.array[:] = u0.x.array[:]
39
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)
45
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.)
51
52
53  # iteração no tempo
54  for k in range(nt):
55
56      t += ht
57      print(f"{k+1}: t = {t:.4g}")
58
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
69
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)
75
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)
81
82      u0.x.array[:] = u.x.array[:]

1.7.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!