Método dos Elementos Finitos - Equação do Calor

Equação do Calor (1D)

A equação do calor em 1D descreve a evolução da temperatura $u(x,t)$ ao longo do tempo em uma barra unidimensional, sujeita a condições de contorno e a uma difusão térmica ($\alpha$). A equação diferencial que governa o comportamento do sistema é dada por:

Forma Forte

$\frac{\partial u(x,t)}{\partial t}=\alpha \frac{{\partial }^{2}u(x,t)}{\partial x^2}$

Condição Inicial

$u(t=0,x)=sin(\pi x)$

Figura 1: Condição inicial para a equação do calor.

Condição de Contorno

$u(t,x=0)=0=u(t,x=L)$

Método dos Elementos Finitos (FEM)

O Método dos Elementos Finitos (FEM, do inglês Finite Element Method) é uma técnica amplamente utilizada para resolver equações diferenciais parciais, como a equação do calor. Para aplicar o FEM à equação do calor em 1D, a barra [0,1] é dividida em um número finito de elementos (intervalos) de comprimento $\mathrm{\Delta }x$, e as variáveis contínuas são aproximadas por funções de base definidas localmente em cada elemento.

Forma Fraca

Apesar de o problema estar completo na forma forte, que é a mostrada acima, precisamos representá-lo na forma equivalente, chamada forma fraca. A forma fraca é preferível por diversos motivos, sendo um dos principais o fato de que ela não possui derivadas de segunda ordem, possibilitando o uso de funções lineares para fazer a aproximação, que na forma forte teriam essa derivada nula. Para obter a forma fraca a partir da forma forte, deve-se seguir os seguintes passos:

• Multiplica-se por uma função teste $v$ do mesmo espaço de funções definido nos elementos.

$\frac{\partial u(x,t)}{\partial t}v=\alpha \frac{{\partial }^{2}u(x,t)}{\partial x^2}v$

$\underset{0}{\overset{L}{\int }}\frac{\partial u(x,t)}{\partial t}v(x)dx=\alpha \underset{0}{\overset{L}{\int }}\frac{{\partial }^{2}u(x,t)}{\partial x^2}v(x)dx$

$\underset{0}{\overset{1}{\int }}\frac{{\partial }^{2}u}{\partial x^2}vdx={\left[\frac{\partial u}{\partial x}v\right]}_0^1-\underset{0}{\overset{1}{\int }}\frac{\partial u}{\partial x}\frac{\partial v}{\partial x}dx$

$\underset{0}{\overset{1}{\int }}\frac{\partial u}{\partial t}vdx+\underset{0}{\overset{1}{\int }}\alpha \frac{\partial u}{\partial x}\frac{\partial v}{\partial x}dx={\left[\alpha \frac{\partial u}{\partial x}v\right]}_0^1$

O lado direito zera devido às condições de contorno

$\underset{0}{\overset{1}{\int }}\frac{\partial u}{\partial t}vdx+\underset{0}{\overset{1}{\int }}\alpha \frac{\partial u}{\partial x}\frac{\partial v}{\partial x}dx=0$

Essa é a forma fraca, mas ainda precisamos discretizá-la tanto no espaço quanto no tempo para aplicar o FEM.

Fenics

O Fenics é uma plataforma de computação popular e de código aberto para resolver equações diferenciais parciais (PDEs) usando o método dos elementos finitos (FEM). Na prática, ajuda a discretizar a parte espacial e a resolver o sistema de equações lineares.

Discretização Temporal

Devemos discretizar a parte temporal para obter a forma fraca final para usar o Fenics. Para isso, usaremos o método de Euler Implícito.

Fazemos a aproximação:

$\frac{\partial u}{\partial t}\approx \frac{u^{t+1}-u^t}{\mathrm{\Delta }t}$

Substituindo na forma fraca: $\underset{0}{\overset{L}{\int }}\frac{u^{t+1}-u^t}{\mathrm{\Delta }t}vdx+\underset{0}{\overset{L}{\int }}\alpha \frac{\partial u^{t+1}}{\partial x}\frac{\partial v}{\partial x}dx=0$

Multiplicando por $\mathrm{\Delta }t$: $\underset{0}{\overset{L}{\int }}{u}^{t+1}vdx-\underset{0}{\overset{L}{\int }}{u}^{t}vdx+\mathrm{\Delta }t\underset{0}{\overset{L}{\int }}\alpha \frac{\partial u^{t+1}}{\partial x}\frac{\partial v}{\partial x}dx=0$

Reorganizando os termos, temos: $\underset{0}{\overset{L}{\int }}{u}^{t+1}vdx+\mathrm{\Delta }t\underset{0}{\overset{L}{\int }}\alpha \frac{\partial u^{t+1}}{\partial x}\frac{\partial v}{\partial x}dx=\underset{0}{\overset{L}{\int }}{u}^{t}vdx$

Para fazer a implementação no Fenics, devemos escrever na versão mais geral, com dimensão arbitrária, utilizando gradiente e produto escalar:

$\underset{\mathrm{\Omega }}{\int }{u}^{t+1}vd\mathrm{\Omega }+\mathrm{\Delta }t\underset{\mathrm{\Omega }}{\int }\alpha \mathrm{\nabla }{u}^{t+1}\cdot \mathrm{\nabla }vd\mathrm{\Omega }=\underset{\mathrm{\Omega }}{\int }{u}^{t}vd\mathrm{\Omega }$

Resultados

Como esperado pela solução analítica, a partir das condições iniciais e de contorno escolhidas, a amplitude do seno diminui exponencialmente com o tempo

Figura 2: Evolução temporal para a equação do calor utilizando o método dos elementos finitos, com dt=0.001.