| | | |

Método dos Elementos Finitos

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

1.4 Seleção de aplicações

1.4.1 Problema não-linear

Como um exemplo de aplicação do método de elementos finitos na solução de equações diferenciais parciais não-lineares, consideramos o problema de Troesch

u′′=μsinh(μu),x(0,1), (1.186)
u(0)=0,u(1)=1, (1.187)

onde μ>0 é um parâmetro.

A formulação fraca deste problema consiste em: encontrar uV:={vH1(I);v(0)=0,v(1)=1} tal que

F(u,v):=L(v)a(u,v)=0,vV0, (1.188)

onde, V0={vH1(I);v(0)=0,v(1)=0} e as formas bilinear a(,) e linear L() são dadas por

a(u,v):=01uv𝑑x, (1.189)
L(v):=01μsinh(μu)v𝑑x. (1.190)

A formulação de elementos finitos deste problema consiste em: encontrar uhVh:={vhP1(I);vh(0)=0,vh(1)=1} tal que

a(uh,vh)=L(vh),vhVh,0, (1.191)

onde Vh,0={vhP1(I);vh(0)=0,vh(1)=0}.

Código 10: ex_mef1d_nlinear.py
1from mpi4py import MPI
2
3# malha
4from dolfinx import mesh
5domain = mesh.create_unit_interval(MPI.COMM_WORLD,
6 nx = 10)
7# espaço
8from dolfinx import fem
9V = fem.functionspace(domain, ('P', 1))
10
11# condição de contorno
12import numpy as np
13uD = fem.Function(V)
14# uD(0) = 0, uD(1) = 1
15def dirichlet_bc(x):
16 y = np.full(x.shape[1], 1.0)
17 y[np.isclose(x[0], 0.)] = 0.0
18 return y
19uD.interpolate(dirichlet_bc)
20
21tdim = domain.topology.dim
22fdim = tdim - 1
23domain.topology.create_connectivity(fdim, tdim)
24boundary_facets = mesh.exterior_facet_indices(domain.topology)
25boundary_dofs = fem.locate_dofs_topological(V, fdim,
26 boundary_facets)
27bc = fem.dirichletbc(uD, boundary_dofs)
28
29# problema mef
30import ufl
31from dolfinx import default_scalar_type
32from dolfinx.fem.petsc import NonlinearProblem
33from dolfinx.nls.petsc import NewtonSolver
34
35# unknown / test functions
36uh = fem.Function(V)
37v = ufl.TestFunction(V)
38
39# parameter
40mu = fem.Constant(domain, default_scalar_type(5))
41
42F = ufl.dot(ufl.grad(uh), ufl.grad(v)) * ufl.dx + (mu * ufl.sinh(mu * uh)) * v * ufl.dx
43
44## manually add Jacobian
45du = ufl.TrialFunction(V)
46J = ufl.dot(ufl.grad(du), ufl.grad(v)) * ufl.dx + (mu**2 * ufl.cosh(mu * uh) * du) * v * ufl.dx
47
48petsc_options = {
49 "snes_type": "newtonls",
50 "snes_linesearch_type": "none",
51 "snes_atol": 1e-6,
52 "snes_rtol": 1e-6,
53 "snes_monitor": None,
54 "ksp_error_if_not_converged": True,
55 "ksp_type": "gmres",
56 "ksp_rtol": 1e-8,
57 "ksp_monitor": None,
58}
59
60problem = NonlinearProblem(
61 F,
62 uh,
63 bcs=[bc],
64 petsc_options=petsc_options,
65 petsc_options_prefix="nonlinpoisson",
66)
67
68problem.solve()
69converged = problem.solver.getConvergedReason()
70num_iter = problem.solver.getIterationNumber()
71assert converged > 0, f"Solver did not converge, got {converged}."
72print(
73 f"Solver converged after {num_iter} iterations with converged reason {converged}."
74)

1.4.2 EDP evolutiva

Vamos estudar a aplicação do método de elementos finitos na solução de equações diferenciais parciais evolutivas no tempo. Para isso, consideramos como problema modelo a equação do calor com dadas condição inicial e condições de contorno de Dirichlet homogêneas

ut=αuxx+f,(t,x)(0,tf]×(a,b), (1.192)
u(0,x)=u0(x),x[a,b], (1.193)
u(t,a)=u(t,b)=0,t[0,tf], (1.194)

onde f=f(t,x) denota uma dada fonte, α>0 é a difusividade e u0 é a condição inicial e tf>0 é o tempo final.

Antes de elaborarmos a formulação fraca, vamos assumir uma 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 θ e denotando u~(k)u(t(k),x) e f(k)=f(t(k),x), o problema pode ser aproximado pela iteração

u~(k+1)u~(k)ht=θ(αu~xx(k+1)+f(k+1))
(1θ)(αu~xx(k)+f(k)), (1.195)
u~(k+1)(a)=u~(k+1)(b)=0, (1.196)

para k=0,1,2,,nt e u~(0)=u0. Ou seja, cada iteração k se resume ao seguinte problema de valores de contorno

αθu~xx+1htu~=1htu~0+(1θ)αu~xx0
(1θ)f0+θf, (1.197)
u~(a)=u~(b)=0, (1.198)

onde u~0=u~(k) e u~=u~(k+1) para fins de simplificação da notação121212Análogo para f(k) e f(k+1)..

A formulação fraca deste problema de valores de contorno consiste em: encontrar u~V0:=H01(a,b) tal que

a(u~,v)=L(v),vV0, (1.199)

onde

a(u~,v):=θabαu~v𝑑x+ab1htu~v𝑑x, (1.200)
L(v):=(θ1)abαu~0v𝑑x+ab1htu~0v𝑑x
θabfv𝑑x+(1θ)abf0v𝑑x. (1.201)

Logo, a formulação de elementos finitos P1 deste problema consiste em: encontrar u~hVh,0:={vhP1(a,b);vh(a)=0,vh(b)=0} tal que

a(u~h,vh)=L(vh),vhVh,0. (1.202)
Exemplo 1.4.1.

Consideremos o problema de calor

ut=uxx+f,(t,x)(0,1]×(0,1), (1.203)
u(0,x)=sen(πx),x[0,1], (1.204)
u(t,0)=u(t,1)=0,t[0,1], (1.205)

onde f(t,x)=(π21)etsen(πx). A solução analítica deste problema é u(t,x)=etsen(πx). O Código 11 a seguir implementa a aproximação por elementos finitos deste problema. A Figura 1.6 apresenta uma comparação da solução analítica com a aproximação de elementos finitos P1 computada com n=40 células. Verifique!

Refer to caption
Figura 1.6: Solução analítica e aproximação por elementos finitos do problema de calor do Exemplo 1.4.1.
Código 11: ex_mef1d_calor.py
1from mpi4py import MPI
2
3# esquema temporal
4thta = 0.5
5nt = 10
6t0 = 0.
7tf = 1.
8dt = (tf - t0) / nt
9
10# malha espacial
11from dolfinx import mesh
12nx = 40
13domain = mesh.create_unit_interval(MPI.COMM_WORLD,
14 nx = nx)
15
16# espaço
17from dolfinx import fem
18V = fem.functionspace(domain, ('P', 1))
19
20# condição inicial
21import numpy as np
22u0 = fem.Function(V)
23u0.interpolate(lambda x: np.sin(np.pi * x[0]))
24
25# condição de contorno
26uD = fem.Function(V)
27uD.interpolate(lambda x: np.full(x.shape[1], 0.))
28
29tdim = domain.topology.dim
30fdim = tdim - 1
31domain.topology.create_connectivity(fdim, tdim)
32boundary_facets = mesh.exterior_facet_indices(domain.topology)
33boundary_dofs = fem.locate_dofs_topological(V, fdim,
34 boundary_facets)
35bc = fem.dirichletbc(uD, boundary_dofs)
36
37# fonte
38import ufl
39x = ufl.SpatialCoordinate(domain)
40def f(t, x):
41 return (ufl.pi**2 - 1) * ufl.exp(-t) * ufl.sin(ufl.pi * x[0])
42
43# iteração temporal
44from dolfinx.fem.petsc import LinearProblem
45u = ufl.TrialFunction(V)
46v = ufl.TestFunction(V)
47
48t = t0
49for k in range(nt):
50 a = u * v * ufl.dx + thta * dt * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
51 L = u0 * v * ufl.dx - (1 - thta) * dt * ufl.dot(ufl.grad(u0), ufl.grad(v)) * ufl.dx + \
52 thta * dt * f(t, x) * v * ufl.dx + (1 - thta) * dt * f(t + dt, x) * v * ufl.dx
53 t += dt
54 print(f'Iteração {k+1}/{nt}, t={t:.3f}')
55 problem = LinearProblem(a, L, bcs=[bc],
56 petsc_options_prefix="calor1d")
57 uh = problem.solve()
58 u0.x.array[:] = uh.x.array[:]

1.4.3 Sistema de equações

Consideramos o seguinte problema de equações diferenciais ordinárias com valores de contorno

u0′′+u1=f0,xI=(iL,iR) (1.206)
u1′′u0=f1,xI (1.207)
u0(iL)=u0L,u0(iR)=u0R, (1.208)
u1(iL)=u1L,u1(iR)=u1R, (1.209)

onde f0, f1, u0L, u0R, u1L, u1R são dados.

Para construirmos uma aproximação por elementos finitos podemos tomar o seguinte problema fraco associado: encontrar 𝒖=(u0,u1)𝒰:=U0×U1 tal que

a(𝒖,𝒗)=L(𝒗),𝒗=(v0,v1)𝒱0:=V0×V0, (1.210)

onde U0={vH1(I);v(iL)=u0L,v(iR)=u0R}, U1={vH1(I);v(iL)=u1L,v(iR)=u1R}, V0={vH1(I);v(0)=v(L)=0}, a forma bilinear é

a(𝒖,𝒗)=Iu0v0𝑑x+Iu1v1𝑑x+Iu1v0𝑑xIu0v1𝑑x (1.211)

e a forma linear é

L(𝒗)=If0v0𝑑x+If1v1𝑑x. (1.212)

Então, o problema de elementos finitos associado no espaço das funções lineares por partes lê-se: encontrar 𝒖h=(uh0,uh1)𝒰h:=Uh0×Uh1 tal que

a(𝒖h,𝒗h)=L(𝒗h),𝒗h=(vh0,vh1)𝒱0:=Vh0×Vh0, (1.213)

onde Uh0={vhP1(I);vh(iL)=u0L,vh(iR)=u0R}, Uh1={vhP1(I);vh(iL)=u1L,vh(L)=u1R}, Vh0={vhP1(I);vh(0)=vh(L)=0}.

Exemplo 1.4.2.

Consideremos o seguinte problema de valor de contorno

u0′′+u1=sen(x)+cos(x),x(π,π) (1.214)
u1′′u0=cos(x)sin(x),x(π,π) (1.215)
u0(π)=0,u0(π)=0, (1.216)
u1(π)=1,u1(π)=1. (1.217)

Sua solução analítica é u0(x)=sen(x) e u1(x)=cos(x). O Código 12 a seguir implementa a aproximação por elementos finitos deste problema. A Figura 1.7 apresenta uma comparação das soluções analíticas e das aproximações de elementos finitos uh0 e uh1, estas construídas no espaço dos polinômios lineares por partes sobre uma malha uniforme de 10 células. Verifique!

Refer to caption
Figura 1.7: Solução analítica versus aproximação de elementos finitos para o problema do Exemplo 1.4.2.
Código 12: ex_mef1d_sistema.py
1from mpi4py import MPI
2
3# malha
4import numpy as np
5from dolfinx import mesh
6domain = mesh.create_interval(MPI.COMM_WORLD,
7 nx = 10,
8 points=np.array([-np.pi, np.pi]))
9# espaço
10from basix.ufl import element, mixed_element
11from dolfinx import fem
12P1 = element("Lagrange", domain.basix_cell(), 1)
13V = fem.functionspace(domain, mixed_element([P1, P1]))
14
15# condição de contorno
16import ufl
17uD = fem.Function(V)
18uD.sub(0).interpolate(lambda x: np.full(x.shape[1], 0.))
19uD.sub(1).interpolate(lambda x: np.full(x.shape[1], -1.))
20
21tdim = domain.topology.dim
22fdim = tdim - 1
23domain.topology.create_connectivity(fdim, tdim)
24boundary_facets = mesh.exterior_facet_indices(domain.topology)
25boundary_dofs = fem.locate_dofs_topological(V, fdim,
26 boundary_facets)
27bc = fem.dirichletbc(uD, boundary_dofs)
28
29# problema mef
30from dolfinx.fem.petsc import LinearProblem
31u = ufl.TrialFunction(V)
32u0, u1 = ufl.split(u)
33v = ufl.TestFunction(V)
34v0, v1 = ufl.split(v)
35
36# fonte
37x = ufl.SpatialCoordinate(domain)
38f0 = ufl.sin(x[0]) + ufl.cos(x[0])
39f1 = ufl.cos(x[0]) - ufl.sin(x[0])
40
41# rhs
42lhs = ufl.dot(ufl.grad(u0), ufl.grad(v0)) * ufl.dx + ufl.dot(ufl.grad(u1), ufl.grad(v1)) * ufl.dx + ufl.dot(u1, v0) * ufl.dx - ufl.dot(u0, v1) * ufl.dx
43
44rhs = f0 * v0 * ufl.dx + f1 * v1 * ufl.dx
45
46problem = LinearProblem(lhs, rhs, bcs=[bc],
47 petsc_options_prefix="ex_mef1d_sistema")
48uh = problem.solve()

1.4.4 Exercícios

E. 1.4.1.

Obtenha uma aproximação por elementos finitos da solução do seguinte problema não-linear de Poisson131313Siméon Denis Poisson, 1781 - 1840, matemático francês. Fonte: Wikipédia:Siméon Denis Poisson.

(q(u)u)=π2(3sen3(πx)sen(πx)),x(0,1), (1.218)
u(0)=0,u(1)=0, (1.219)

onde q(u)=1+u2. Faça o esboço dos gráficos da solução analítica u(x)=sen(πx) e da sua aproximação de elementos finitos. Então, faça uma análise da convergência de malha e verifique as taxas de convergência dos erros uuhL2 e (uuh)L2.

E. 1.4.2.

Considere o problema de Fisher

ut=uxx+λu(1u),(t,x)(0,tf]×(0,1), (1.220)
u(0,x)=(1+eλ6x)2,x[0,1], (1.221)
ux(t,0)=(1+e56λt)2,t[0,tf], (1.222)
ux(t,1)=(1+eλ656λt)2,t[0,tf], (1.223)

onde λ=1 (taxa de crescimento populacional) e tf=1 (tempo final). A solução analítica deste problema é

u(t,x)=(1+eλ6x56λt)2. (1.224)

Faça três códigos de elementos finitos para este problema utilizando os seguintes esquemas de tempo:

  1. a)

    θ=0 (esquema de Euler explícito);

  2. b)

    θ=1 (esquema de Euler implícito);

  3. c)

    θ=0.5 (esquema de Crank-Nicolson).

Em cada caso, faça o esboço dos gráficos da solução analítica e da sua aproximação de elementos finitos. Então, faça uma análise da convergência das malhas espacial e temporal.

E. 1.4.3.

Considere o seguinte problema não-linear de Calor

ut=(q(u)ux)+f,(t,x)(0,tf]×(0,1), (1.225)
u(0,x)=sen(πx),x[0,1], (1.226)
u(t,0)=u(t,1)=0,t[0,tf], (1.227)

onde q(u)=1+u2, a fonte é

f(t,x)=(π21)ets+e3ts(3s22), (1.228)

com s:=sin(πx). A solução analítica deste problema é u(t,x)=etsen(πx). Faça três códigos de elementos finitos para este problema utilizando os seguintes esquemas de tempo:

  1. a)

    Runge-Kutta de ordem 4 (RK4 - explícito);

  2. b)

    Método de Adams-Moulton de ordem 2 (A-M-2);

  3. c)

    Método preditor-corretor de Adams-Bashforth-Moulton de ordem 2 (A-B-M-2);

Quais dos códigos mostrou-se mais eficiente?

E. 1.4.4.

Estude aproximações por elementos finitos para o seguinte problema de Burgers

ut+uux=νuxx,(t,x)(0,1]×(1,1), (1.229)
u(0,x)=2νπsen(πx)2+cos(πx),x[1,1], (1.230)
u(t,1)=u(t,1)=0,t[0,1], (1.231)

para ν=0.1,0.01,0.001. A solução analítica deste problema é

u(t,x)=2νπeνπ2tsen(πx)2+eνπ2tcos(πx). (1.232)
E. 1.4.5.

Estude aproximações por elementos finitos para o seguinte problema de Burgers

ut+uux=νuxx,(t,x)(0,1]×(1,1), (1.233)
u(0,x)=sen(πx),x[1,1], (1.234)
u(t,1)=u(t,1)=0,t[0,1], (1.235)

onde ν>0. Para ν=0.01/π, a solução analítica é

u(t,x)=sen(π(xη))f(xη)eη24νtdηf(xη)eη24νt𝑑η, (1.236)

onde f(y)=ecos(πy)2πν. Estude também o comportamento das aproximações por elementos finitos para ν cada vez menor.

E. 1.4.6.

Estude aproximações por elementos finitos para o seguinte problema de Burgers

ut+uux=νuxx,(t,x)(0,1]×(1,1), (1.237)
u(0,x)=2sign(x),x[1,1], (1.238)
u(t,1)=ua(t,1),u(t,1)=ua(t,1),t[0,1], (1.239)

onde ν>0. Para ν=1, a solução analítica é

ua(t,x)=2G(t,x)G(t,x)G(t,x)+G(t,x), (1.240)

onde G(t,x)=12etxerfc(2tx2t). Estude também o comportamento das aproximações por elementos finitos para ν cada vez menor.

E. 1.4.7.

Estude aproximações por elementos finitos para o seguinte problema de Solow-Swan

Kt=s(x)AKϕL1ϕ(x)δK+d2Kx2,x(0,l),t>0, (1.241)
K(t,x)=K0(x),xΩ=[0,l],t=0, (1.242)
Kx=0,xΩ={0,l},t>0. (1.243)

para o caso em que l10, T10, ϕ1/3, δ(x)0.05, A(x)1, d(x)0.25, L(x)=1+0.3x2[1cos(4πx/l)], k0(x)=cos4(πx/(2l)π/2), and s(x)=0.2cos4(πx/(2l)).

E. 1.4.8.

Estude aproximações por elementos finitos para o seguinte sistema de equações de difusão-reação

ε0u0′′+au0+bu1=f0,x(0,1), (1.244)
ε1u1′′+cu0+du1=f1,x(0,1), (1.245)
u0(0)=u0(1)=0, (1.246)
u1(0)=u1(1)=0, (1.247)

onde ε0,ε1>0 são os coeficientes de difusão, a,b,c,d são os coeficientes de reação, e f0,f1 são as fontes. Como estudo de caso, considere ε0=1, ε1=2, a=1, b=1, c=1, d=2. Para a solução manufaturada u0(x)=sen(πx), u1(x)=cos(πx), as fontes são

f0(x)=π2sen(πx)2cos(π(x+1/4)), (1.248)
f1(x)=sen(πx)+2cos(πx)+2π2cos(πx). (1.249)

Estude a convergência de malha e verifique as taxas de convergência dos erros u0uh0L2, (u0uh0)L2, u1uh1L2 e (u1uh1)L2.

E. 1.4.9.

Estude aproximações por elementos finitos para o seguinte sistema de equações de difusão-reação não-linear

ε0u0′′u0+u02u1=f0,x(0,1), (1.250)
ε1u1′′u0u12=f1,x(0,1), (1.251)
u0(0)=π,u0(1)=π, (1.252)
u1(0)=1,u1(1)=1. (1.253)

onde ε0,ε1>0 são os coeficientes de difusão, e f0,f1 são as fontes. Como estudo de caso, considere ε0=1, ε1=10, e as seguintes fontes

f0(x)=(sen(πx)cos(2πx)1+π2)sen(πx), (1.254)
f1(x)=(sen(πx)cos(2πx)+40π2)cos(2πx). (1.255)

A solução analítica deste problema é u0(x)=sen(πx) e u1(x)=cos(2πx). Estude a convergência de malha e verifique as taxas de convergência dos erros u0uh0L2, (u0uh0)L2, u1uh1L2 e (u1uh1)L2.

E. 1.4.10.

Estude aproximações por elementos finitos para o seguinte problema de FitzHugh-Nagumo

ut=Duuxx+λuu3vκ,(t,x)(0,tf]×(L,L), (1.256)
vt=Dvvxx+uv,(t,x)(0,tf]×(L,L), (1.257)
ux(t,L)=ux(t,L)=0,t[0,tf], (1.258)
vx(t,L)=vx(t,L)=0,t[0,tf], (1.259)
u(0,x)=ex2,v(0,x)=0,x[L,L], (1.260)

onde Du,Dv>0 são os coeficientes de difusão, λ é o parâmetro de excitabilidade, e κ é o parâmetro de estímulo. Obtenha uma solução precisa para Du=0.1, Dv=0.01, λ=1, κ=0, L=10 e tf=10.


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!

Opcional. Preencha seu nome para que eu possa lhe contatar.
Opcional. Preencha seu e-mail para que eu possa lhe contatar.
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.

Licença Creative Commons
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
| | | |