Método dos Elementos Finitos

3 Aplicações 2D 3 Aplicações 2D 3.2 Equação de advecção-difusão-reação

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

3.1 Equação de advecção-difusão

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

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

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

(3.1)

onde $\kappa>0$ é o coeficiente de difusão, $\boldsymbol{a}$ é o campo de velocidade, $f$ é a fonte e $\Gamma$ é a fronteira do domínio $\Omega$ . Por simplicidade, consideramos que $\boldsymbol{a}$ é um campo de velocidade constante; os desenvolvimentos a seguir podem ser estendidos a campos com divergência nula, i.e.

\displaystyle\nabla\cdot\boldsymbol{a}=0.

(3.2)

Além disso, assumimos condições de contorno de Dirichlet homogêneas

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

(3.3)

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% )=\left(f,v\right)\quad\forall v\in H_{0}^{1}(\Omega).

(3.4)

Tomando $v=u$ em (3.4), a parte de advecção pode ser reescrita usando a identidade de Gauss¹¹1Johann Carl Friedrich Gauss, 1777 - 1855, matemático alemão. Fonte: Wikipédia: Carl Friedrich Gauss.:

	$\displaystyle\left(\boldsymbol{a}\cdot\nabla u,u\right)=\int_{\Omega}\left(% \boldsymbol{a}\cdot\nabla u\right)\,u\,dx$
	$\displaystyle\text{}\quad=\frac{1}{2}\int_{\Omega}\nabla\cdot(\boldsymbol{a}u^% {2})\,dx$
	$\displaystyle\text{}\quad=\frac{1}{2}\int_{\Gamma}\boldsymbol{a}\cdot% \boldsymbol{n}u^{2}\,ds$
	$\displaystyle\text{}\qquad-\frac{1}{2}\int_{\Omega}(\nabla\cdot\boldsymbol{a})% u^{2}\,dx=0.$

Com isso, a expressão em (3.4) pode ser limitada, à esquerda, por Poincaré²²2Henri Poincaré, 1854 - 1912, matemático francês. Fonte: Wikipédia: Henri Poincaré.–Friedrichs³³3Kurt Otto Friedrichs, 1901 - 1982, matemático estadunidense-alemão. Fonte: Wikipédia: Kurt Otto Friedrichs. e, à direita, por Cauchy⁴⁴4Augustin-Louis Cauchy, 1789-1857, matemático francês. Fonte: Wikipédia: Augustin-Louis Cauchy.–Schwarz⁵⁵5Karl Hermann Amandus Schwarz, 1843-1921, matemático alemão. Fonte: Wikipédia: Hermann Amandus Schwarz., resultando em:

\displaystyle\kappa C\|u\|_{H^{1}_{0}}^{2}\leq\kappa\|\nabla u\|^{2}=\left(f,u% \right)\leq\|f\|\|u\|\leq\|f\|\|u\|_{H^{1}_{0}}.

(3.5)

Assim, temos a estimativa de estabilidade

\displaystyle\|u\|_{H^{1}_{0}}\leq\frac{1}{\kappa C}\|f\|.

(3.6)

Esta estimativa mostra que a constante de estabilidade depende inversamente do coeficiente de difusão $\kappa$ . Como consequência, temos a indicação de que para valores pequenos de $\kappa$ , pequenas perturbações em $f$ podem levar a grandes variações em $u$ (instabilidade).

Classicamente, o problema de elementos finitos associado é: 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)=\left(f,v_{h}\right)\quad\forall v_{h}\in V_{h}.

(3.7)

Este problema herda a instabilidade do problema fraco, o que pode levar a soluções oscilatórias, especialmente para valores pequenos de $\kappa$ (regime de advecção-dominante).

Exemplo 3.1.1.

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

	$\displaystyle-\kappa\Delta u+\boldsymbol{a}\cdot\nabla u=0,\text{em }(0,1)^{2},$		(3.8)
	$\displaystyle u(0,y)=1,\leavevmode\nobreak\ y\in(0.8,1),$		(3.9)
	$\displaystyle u(0,y)=0,\leavevmode\nobreak\ y\in(0,0.8),$		(3.10)
	$\displaystyle u(x,1)=1,\leavevmode\nobreak\ x\in(0,1),$		(3.11)
	$\displaystyle u_{y}(x,0)=0,\leavevmode\nobreak\ x\in(0,1),$		(3.12)
	$\displaystyle u_{x}(1,y)=0,\leavevmode\nobreak\ y\in(0,1),$		(3.13)

onde $\boldsymbol{a}=(\sqrt{2}/2,-\sqrt{2}/2)$ e $\kappa=10^{-3}$ . A solução deste problema apresenta uma camada de fronteira ao longo da diagonal do quadrado, o que pode levar a oscilações na solução numérica se o método de elementos finitos clássico for utilizado.

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

Refer to caption — Figura 3.1: Solução numérica do problema de advecção-difusão utilizando o método de elementos finitos clássico.

Código 18: ad_mef2d.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 SpatialCoordinate,

16 TestFunction,

17 TrialFunction,

18 dot,

19 ds,

20 dx,

21 grad

22)

24import numpy as np

26# espaço

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

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

29u = TrialFunction(V)

30v = TestFunction(V)

32# forma bilinear

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

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

35a_form = kappa * dot(grad(u), grad(v)) * dx + dot(a, grad(u)) * v * dx

37# forma linear

38zero = Constant(mesh, 0.0)

39L_form = zero * v * dx

41# condições de Dirichlet

42x = SpatialCoordinate(mesh)

44def boundary_left_top(x):

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

46dofs_bc1 = locate_dofs_geometrical(V, boundary_left_top)

47one = Constant(mesh, 1.0)

48bc1 = dirichletbc(one, dofs_bc1, V)

50def boundary_left_bottom(x):

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

52dofs_bc2 = locate_dofs_geometrical(V, boundary_left_bottom)

53bc2 = dirichletbc(zero, dofs_bc2, V)

55def boundary_top(x):

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

57dofs_bc3 = locate_dofs_geometrical(V, boundary_top)

58bc3 = dirichletbc(one, dofs_bc3, V)

60bcs = [bc1, bc2, bc3]

62# problema linear

63solver_options = {

64 "ksp_type": "gmres",

65 "ksp_rtol": 1.0e-8,

66 "ksp_atol": 1.0e-10,

67 "ksp_max_it": 1000,

68 "ksp_monitor": None,

69 "ksp_converged_reason": None,

70 "pc_type": "jacobi",

71}

72problem = LinearProblem(

73 a_form,

74 L_form,

75 bcs=bcs,

76 petsc_options_prefix="ad_mef2d",

77 petsc_options=solver_options

78)

79uh = problem.solve()

81# plot 3D surface

82import pyvista

83from dolfinx.plot import vtk_mesh

85topology, cell_types, geometry = vtk_mesh(V)

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

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

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

91plotter = pyvista.Plotter()

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

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

94plotter.view_vector((1, 1, 1))

95plotter.add_axes()

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

97plotter.show()

3.1.1 SUPG

Uma formulação de elementos finitos estabilizada para o problema de advecção-difusão (3.1) é a formulação SUPG (Streamline Upwind Petrov⁶⁶6Georgi Iwanowitsch Petrov, 1912 - 1987, engenheiro soviético. Fonte: Wikipedia: Georgi Iwanowitsch Petrov.-Galerkin⁷⁷7Boris Galerkin, 1871 - 1945, engenheiro e matemático soviético. Fonte: Wikipédia: Boris Galerkin.). A formulação SUPG é 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)+\sum_{K}\tau_{K}\left(-\kappa\Delta u_{h}+% \boldsymbol{a}\cdot\nabla u_{h},\boldsymbol{a}\cdot\nabla v_{h}\right)_{K}$
		$\displaystyle\text{}\quad=\left(f,v_{h}\right)+\sum_{K}\tau_{K}\left(f,% \boldsymbol{a}\cdot\nabla v_{h}\right)_{K},$		(3.14)

onde $\tau_{K}$ é um parâmetro de estabilização que depende do tamanho do elemento $K$ . Os termos adicionais na formulação SUPG atuam como um mecanismo de estabilização, reduzindo as oscilações na solução numérica pela adição de um termo de difusividade artificial ao longo das linhas de fluxo do campo de velocidade $\boldsymbol{a}$ . Existem várias formas de escolha do parâmetro $\tau_{K}$ , sendo uma escolha comum a seguinte [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}\right)^{-\frac{1}{2}},

(3.15)

onde $h_{K}$ é o tamanho do elemento $K$ . Para a justificativa desta escolha e de outras alternativas, consulte [5].

Exemplo 3.1.2.

Vamos considerar o mesmo problema do Exemplo 3.1.1, mas agora utilizando a formulação SUPG. A Figura 3.2 mostra a solução numérica obtida utilizando a formulação SUPG com uma triangularização sob uma malha uniforme de $10\times 10$ células.

O Código 18 pode ser facilmente modificado pela inclusão do termo de estabilização SUPG. Consulte a implementação do termo de estabilização no Código 19 e verifique a solução implementando-a no Código 18.

Código 19: ad_mef2d_supg.py

⬇

1from ufl import CellDiameter, sqrt

2# SUPG stabilization

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

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

5tau0 = Constant(mesh, 1.0)

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

7a_form += tau * dot(a, grad(u)) * dot(a, grad(v)) * dx

3.1.2 Exercícios

Em revisão

E. 3.1.1.

Obtenha aproximações de elementos finitos da solução do seguinte problema de advecção-difusão em 1D:

	$\displaystyle-\kappa u^{\prime\prime}+au^{\prime}=10e^{-5x}-4e^{-x},\text{em }% (0,1),$		(3.16)
	$\displaystyle u(0)=0,\leavevmode\nobreak\ u(1)=1.$		(3.17)

Verifique a estabilidade da solução numérica para diferentes valores do número de Péclet $Pe_{K}=\frac{ah_{K}}{2\kappa}=0.25,\dotsc,5$ . Use a formulação SUPG (com escolha de $\tau_{K}$ análoga à (3.15)) para obter soluções estáveis e compare os resultados com a formulação clássica de elementos finitos. Por fim, faça uma análise de convergência de malha para ambos os métodos.

E. 3.1.2.

Para problemas de advecção-difusão em 1D, uma escolha apropriada do parâmetro de estabilização $\tau_{K}$ na formulação SUPG é dada por

\displaystyle\tau_{K}=\frac{h_{K}}{2a}\left(\coth(Pe_{K})-\frac{1}{Pe_{K}}% \right),

(3.18)

onde $Pe_{K}=\frac{ah_{K}}{2\kappa}$ é o número de Péclet local. Implemente esta escolha de $\tau_{K}$ para o problema do E.3.1.1 e compare os resultados com a escolha de $\tau_{K}$ análoga à (3.15).

E. 3.1.3.

Faça uma análise de convergência de malha para o problema de advecção-difusão em 2D apresentado no Exemplo 3.1.1, utilizando tanto a formulação clássica de elementos finitos quanto a formulação SUPG. Uma escolha alternativa para o parâmetro de estabilização $\tau_{K}$ na formulação SUPG é dada por

\displaystyle\tau_{K}=\left(\frac{2a}{h}+\frac{4\kappa}{h^{2}}\right)^{-1}.

(3.19)

Para o problema em questão, esta escolha é melhor que a sugerida em (3.15)? 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