Método dos Elementos Finitos

3 Aplicações 2D 3.1 Equação de advecção-difusão Referências

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

3.2 Equação de advecção-difusão-reação

Vamos considerar a equação de advecção-difusão-reação em 2D, dada por:

\displaystyle-\kappa\Delta u+\boldsymbol{a}\cdot\nabla u+\sigma u=f\quad\text{% em }\Omega,

\displaystyle u=0\quad\text{em }\Gamma,

(3.20)

onde $\sigma\geq 0$ é o coeficiente de reação. A formulação variacional fraca para este problema é: encontrar $u\in H_{0}^{1}(\Omega)$ tal que

\kappa\left(\nabla u,\nabla v\right)+\left(\boldsymbol{a}\cdot\nabla u,v\right% )+\sigma\left(u,v\right)=\left(f,v\right)\quad\forall v\in H_{0}^{1}(\Omega).

(3.21)

A presença do termo de reação $\sigma u$ pode melhorar a estabilidade do problema, especialmente para valores grandes de $\sigma$ . No entanto, para valores pequenos de $\kappa$ e $\sigma$ , o problema pode apresentar instabilidades semelhantes ao caso de advecção-difusão (Seção 3.1).

Uma formulação de elementos finitos estabilizada para o problema de advecção-difusão-reação é a formulação SUPG⁸⁸8Consulte o E.3.2.2. Uma formulação alternativa é a formulação GLS que estudaremos a seguir.

3.2.1 GLS

O método GLS (Galerkin/Least-Squares) é uma formulação estabilizada de elementos finitos baseada na adição de um termo de penalização que minimiza o resíduo do problema diferencial. A formulação GLS para o problema de advecção-difusão-reação (3.20) é dada por: encontrar $u_{h}\in V_{h}$ tal que

	$\displaystyle\kappa\left(\nabla u_{h},\nabla v_{h}\right)+\left(\boldsymbol{a}% \cdot\nabla u_{h},v_{h}\right)+\sigma\left(u_{h},v_{h}\right)$
	$\displaystyle\text{}\qquad+\sum_{K}\tau_{K}\left(-\kappa\Delta u_{h}+% \boldsymbol{a}\cdot\nabla u_{h}+\sigma u_{h},-\kappa\Delta v_{h}+\boldsymbol{a% }\cdot\nabla v_{h}+\sigma v_{h}\right)_{K}$
	$\displaystyle\text{}\quad=\left(f,v_{h}\right)+\sum_{K}\tau_{K}\left(f,-\kappa% \Delta v_{h}+\boldsymbol{a}\cdot\nabla v_{h}+\sigma v_{h}\right)_{K},$		(3.22)

onde $\tau_{K}$ é um parâmetro de estabilização que depende do tamanho do elemento $K$ . Alternativas comumente utilizadas para a escolha de $\tau_{K}$ para este problema são as seguintes [5, Subseção 2.4.3]:

	$\displaystyle\tau_{K}=\left(\left(\frac{2\\|\boldsymbol{a}\\|}{h_{K}}\right)^{2}% +9\left(\frac{4\kappa}{h_{K}^{2}}\right)^{2}+\sigma^{2}\right)^{-\frac{1}{2}},$		(3.23)
	$\displaystyle\tau_{K}=\left(\frac{2\\|\boldsymbol{a}\\|}{h_{K}}+\frac{4\kappa}{h% _{K}^{2}}+\sigma\right)^{-1}.$		(3.24)

Exemplo 3.2.1.

Vamos considerar o seguinte problema de advecção-difusão-reação em 2D:

	$\displaystyle-\kappa\Delta u+\boldsymbol{a}\cdot\nabla u+\sigma u=0,\text{em }% (0,1)^{2},$		(3.25)
	$\displaystyle u(0,y)=1,\leavevmode\nobreak\ y\in(0.8,1),$		(3.26)
	$\displaystyle u(0,y)=0,\leavevmode\nobreak\ y\in(0,0.8),$		(3.27)
	$\displaystyle u(x,1)=1,\leavevmode\nobreak\ x\in(0,1),$		(3.28)
	$\displaystyle u_{y}(x,0)=0,\leavevmode\nobreak\ x\in(0,1),$		(3.29)
	$\displaystyle u_{x}(1,y)=0,\leavevmode\nobreak\ y\in(0,1),$		(3.30)

onde $\boldsymbol{a}=(\sqrt{2}/2,-\sqrt{2}/2)$ , $\kappa=10^{-3}$ e $\sigma=1$ . Observe que este problema é semelhante ao do Exemplo 3.1.1, mas agora com a adição do termo de reação $\sigma u$ .

A Figura 3.3 mostra a solução numérica obtida utilizando o método de elementos finitos clássico implementado no Código 20. A malha utilizada é uma triangularização uniforme sob uma partição de $10\times 10$ células. Verifique!

Refer to caption — Figura 3.3: Solução numérica do problema de advecção-difusão-reação utilizando a formulação GLS do método de elementos finitos.

Código 20: adr_mef2d_gls.py

⬇

1from dolfinx import default_scalar_type

2from dolfinx.fem import (

3 Constant,

4 Function,

5 functionspace,

6 dirichletbc,

7 form,

8 locate_dofs_geometrical

10from dolfinx.fem.petsc import LinearProblem

11from dolfinx.mesh import create_unit_square

13from mpi4py import MPI

14from ufl import (

15 CellDiameter,

16 SpatialCoordinate,

17 TestFunction,

18 TrialFunction,

19 div,

20 dot,

21 dx,

22 grad,

23 sqrt

24)

26import numpy as np

28# espaço

29mesh = create_unit_square(MPI.COMM_WORLD, 10, 10)

30V = functionspace(mesh, ("Lagrange", 1))

31u = TrialFunction(V)

32v = TestFunction(V)

34# forma bilinear

35kappa = Constant(mesh, 1e-3)

36sigma = Constant(mesh, 1.0)

37a = Constant(mesh, np.array([np.sqrt(2)/2, -np.sqrt(2)/2], dtype=default_scalar_type))

38a_form = (

39 kappa * dot(grad(u), grad(v)) * dx

40 + dot(a, grad(u)) * v * dx

41 + sigma * u * v * dx

42)

44# GLS stabilization

45h = CellDiameter(mesh) # local cell size h_K

46a_norm = sqrt(dot(a, a))

47tau0 = Constant(mesh, 1.0)

48tau = tau0 / sqrt((2.0 * a_norm / h)**2 + 9.0 * (4.0 * kappa / h**2 + sigma)**2)

49res_u = -kappa * div(grad(u)) + dot(a, grad(u)) + sigma * u

50res_v = -kappa * div(grad(v)) + dot(a, grad(v)) + sigma * v

51a_form += tau * res_u * res_v * dx

53# forma linear

54zero = Constant(mesh, 0.0)

55L_form = zero * v * dx

57# condições de Dirichlet

58x = SpatialCoordinate(mesh)

60def boundary_left_top(x):

61 return np.logical_and(np.isclose(x[0], 0), x[1] >= 0.8)

62dofs_bc1 = locate_dofs_geometrical(V, boundary_left_top)

63one = Constant(mesh, 1.0)

64bc1 = dirichletbc(one, dofs_bc1, V)

66def boundary_left_bottom(x):

67 return np.logical_and(np.isclose(x[0], 0), x[1] < 0.8)

68dofs_bc2 = locate_dofs_geometrical(V, boundary_left_bottom)

69bc2 = dirichletbc(zero, dofs_bc2, V)

71def boundary_top(x):

72 return np.isclose(x[1], 1)

73dofs_bc3 = locate_dofs_geometrical(V, boundary_top)

74bc3 = dirichletbc(one, dofs_bc3, V)

76bcs = [bc1, bc2, bc3]

78# problema linear

79solver_options = {

80 "ksp_type": "gmres",

81 "ksp_rtol": 1.0e-8,

82 "ksp_atol": 1.0e-10,

83 "ksp_max_it": 1000,

84 "ksp_monitor": None,

85 "ksp_converged_reason": None,

86 "pc_type": "jacobi",

87}

88problem = LinearProblem(

89 a_form,

90 L_form,

91 bcs=bcs,

92 petsc_options_prefix="ad_mef2d",

93 petsc_options=solver_options

94)

95uh = problem.solve()

97# plot 3D surface

98import pyvista

99from dolfinx.plot import vtk_mesh

100

101topology, cell_types, geometry = vtk_mesh(V)

102

103grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

104grid.point_data["uh"] = uh.x.array.real

105surface = grid.warp_by_scalar("uh", factor=1.0)

106

107plotter = pyvista.Plotter()

108plotter.add_mesh(surface, show_edges=False, scalars="uh",

109 scalar_bar_args={"title": "uh", "vertical": True})

110plotter.view_vector((1, 1, 1))

111plotter.add_axes()

112plotter.view_vector((-2, 0, 1))

113plotter.show()

114# save as PNG image

115plotter.screenshot("fig.png")

3.2.2 Exercícios

E. 3.2.1.

Estude a estabilidade do problema de advecção-difusão-reação do Exemplo 3.2.1 para:

a)

um problema de convecção-difusão dominante, com $\kappa=10^{-4}$ e $\sigma=1$ ;
b)

um problema de reação dominante, com $\kappa=10^{-4}$ , $\left\|\boldsymbol{a}\right\|=10^{-3}$ e $\sigma=1$ .

Compare os resultados obtidos com a formulação estabilizada GLS para ambos os casos.

E. 3.2.2.

Implemente a formulação SUPG para o problema de advecção-difusão-reação do Exemplo 3.2.1. Uma escolha apropriada do parâmetro de estabilização $\tau_{K}$ para este problema é dada por

\displaystyle\tau_{K}=\left(\left(\frac{2\|\boldsymbol{a}\|}{h_{K}}\right)^{2}% +9\left(\frac{4\kappa}{h_{K}^{2}}\right)^{2}+\sigma^{2}\right)^{-\frac{1}{2}}.

(3.31)

Compare os resultados obtidos com a formulação clássica de elementos finitos para diferentes valores de $\kappa$ e $\sigma$ .

E. 3.2.3.

Refaça o E.3.2.1 utilizando a formulação SUPG para o problema de advecção-difusão-reação. Qual formulação estabilizada (SUPG ou GLS) fornece melhores resultados para cada um dos casos estudados? Justifique sua resposta.

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