Exibir código-fonte para Método dos Elementos Finitos - Equação do Calor

== Introdução ==

O Método dos Elementos Finitos (FEM) é uma ferramenta numérica amplamente utilizada na resolução de problemas de engenharia e física que envolvem equações diferenciais parciais. 

Dentre suas principais vantagens, o FEM permite abordar problemas complexos, mesmo em situações com geometrias irregulares, condições de contorno variadas e materiais com propriedades heterogêneas.

Neste trabalho, utilizamos o FEM para estudar a equação do calor, que descreve a distribuição de temperatura em um meio ao longo do tempo e do espaço. A escolha da equação do calor em uma dimensão se deve à existência de uma solução analítica exata, o que possibilita validar os resultados numéricos. Além disso, o método pode ser facilmente generalizado para problemas em dimensões superiores e configurações mais complexas, demonstrando sua versatilidade e eficácia em diferentes aplicações.

== Equação do Calor (1D) ==

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

=== Forma Forte ===
<math>\frac{\partial u(x,t)}{\partial t} = \alpha \frac{\partial^2 u(x,t)}{\partial x^2}</math>

==== Condição Inicial ====

A condição inicial escolhida em todas as simulações é:

<math>u(t=0,x)=sin(\pi x)</math>


[[Arquivo:initial condition.png|500px|center|thumb|Figura 1: Condição inicial para a equação do calor.]]

==== Condição de Contorno ====
Foram escolhidas as condições de contorno homogêneas de Dirichlet, ou seja, temperatura nula nas extremidades da barra duranto todo o tempo.

<math>u(t,x=0)=0=u(t,x=L)</math>

===Solução Analítica ===

A solução analítica é dada por:<ref name="analytical">Douglas, J. (1988). Introduction to partial differential equations (Vol. 2). New York: Academic Press. Chapter 4.</ref>

<math>u(x,t) = \sum_{n=1}^\infty b_n e^{-\alpha \left(\frac{n\pi}{L}\right)^2 t} \sin\left(\frac{n\pi x}{L}\right),</math>

Para <math>L=1</math>, <math>\alpha = 1</math> e as condições iniciais acima, todos os coeficientes são nulos exceto <math>b_1=1</math>

Obtemos então:

<math>u(x, t) = e^{-\pi^2 \alpha t} \sin(\pi x),</math>

== 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 <math>\Delta x</math>, 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 <math>v</math> do mesmo espaço de funções definido nos elementos.

<math>\frac{\partial u(x,t)}{\partial t}v = \alpha \frac{\partial^2 u(x,t)}{\partial x^2}v</math>

<math>\int_0^L \frac{\partial u(x,t)}{\partial t} v(x) \, dx = \alpha \int_0^L \frac{\partial^2 u(x,t)}{\partial x^2} v(x) \, dx</math>

<math>\int_0^1 \frac{\partial^2 u}{\partial x^2} v \, dx = 
\left[ \frac{\partial u}{\partial x} v \right]_0^1 - \int_0^1 \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} \, dx</math>

<math>\int_0^1 \frac{\partial u}{\partial t} v \, dx + \int_0^1 \alpha \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} \, dx = \left[ \alpha \frac{\partial u}{\partial x} v \right]_0^1</math>

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

<math>\int_0^1 \frac{\partial u}{\partial t} v \, dx + \int_0^1 \alpha \frac{\partial u}{\partial x} \frac{\partial v}{\partial x} \, dx = 0</math>

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

=== Fenics ===

O FEniCS (Finite Element Computational Software) é uma plataforma de software de código aberto projetada para resolver problemas numéricos usando o Método dos Elementos Finitos (FEM). Ele é amplamente utilizado para resolver equações diferenciais parciais (EDPs). 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:

<math>\frac{\partial u}{\partial t} \approx \frac{u^{t+1} - u^t}{\Delta t}</math>

Substituindo na forma fraca:
<math>\int_0^L \frac{u^{t+1} - u^t}{\Delta t} v \, dx + \int_0^L \alpha \frac{\partial u^{t+1}}{\partial x} \frac{\partial v}{\partial x} \, dx = 0</math>

Multiplicando por <math>\Delta t</math>:
<math>\int_0^L u^{t+1} v \, dx - \int_0^L u^t v \, dx + \Delta t \int_0^L \alpha \frac{\partial u^{t+1}}{\partial x} \frac{\partial v}{\partial x} \, dx = 0</math>

Reorganizando os termos, temos:
<math>\int_0^L u^{t+1} v \, dx + \Delta t \int_0^L \alpha \frac{\partial u^{t+1}}{\partial x} \frac{\partial v}{\partial x} \, dx = \int_0^L u^t v \, dx</math>

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

<math>\int_{\Omega} u^{t+1} v \, d\Omega + \Delta t \int_{\Omega} \alpha \nabla u^{t+1} \cdot \nabla v \, d\Omega = \int_{\Omega} u^t v \, d\Omega</math> <ref name=weakform>Fangohr, H., & Logg, A. (2020). Solving partial differential equations in Python (Capítulo 3). Springer.</ref>

== 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


[[Arquivo:evolution.gif]]

O intervalo de tempo escolhido tem grande influência na acurácia do resultado. Nas figuras 3, 4 e 5 são mostradas a solução obtida pelo FEM em comparação com a solução analítica para diferentes instantes de tempo.
É possível perceber que a solução, independente do intervalo de tempo (dt) escolhido, é muito boa para tempos curtos. Entretanto, quanto maior o tempo, maior será o erro para dt maiores.

[[Arquivo:comparing_t_0.02.png|500px|center|thumb|Figura 2: Comparação da solução analítica com a solução do método dos elementos finitos para diferentes intervalos de tempo e mantendo fixo t=0.02.]]

[[Arquivo:comparing_t_0.05.png|500px|center|thumb|Figura 3: Comparação da solução analítica com a solução do método dos elementos finitos para diferentes intervalos de tempo e mantendo fixo t=0.05.]]

[[Arquivo:comparing_t_0.1.png|500px|center|thumb|Figura 4: Comparação da solução analítica com a solução do método dos elementos finitos para diferentes intervalos de tempo e mantendo fixo t=0.01..]]

[[Arquivo:comparing_t_0.3.png|500px|center|thumb|Figura 5: Comparação da solução analítica com a solução do método dos elementos finitos para diferentes intervalos de tempo e mantendo fixo t=0.3.]]

Alem disso, para demonstrar a influência do número de elementos para a qualidade do resultado, são comparadas nas figuras 6, 7 e 8 a solução analítica com a solução obtida pelo FEM. Para um número de elementos muito pequeno, a curva é uma aproximação grosseira e claramente linear em partes. Quando aumentamos os elementos, entretanto, rapidamente a curva suaviza e obtemos uma ótima aproximação para o resultado analítico. Assim como discutido anteriormente, a aproximação piora com o passar do tempo.

[[Arquivo:comparing_mesh_0.02.png|500px|center|thumb|Figura 6: Comparação de diferentes números de elementos (mesh) com a solução analítica, mantendo fixos dt=0.005 e mantendo fixo t=0.02.]]

[[Arquivo:comparing_mesh_0.1.png|500px|center|thumb|Figura 7: Comparação de diferentes números de elementos (mesh) com a solução analítica, mantendo fixos dt=0.005 e mantendo fixo t=0.1.]]


[[Arquivo:comparing_mesh_0.3.png|500px|center|thumb|Figura 8: Comparação de diferentes números de elementos (mesh) com a solução analítica, mantendo fixos dt=0.005 e mantendo fixo t=0.3.]]

== Conclusão ==

Este trabalho aplicou o Método dos Elementos Finitos (FEM) para resolver a equação do calor unidimensional, utilizando a solução analítica como referência para validar os resultados numéricos. A comparação entre ambas as soluções revelou que o FEM oferece boas aproximações, com precisão que depende da escolha do passo temporal e do número de elementos.

Embora o erro numérico cresça com o tempo, principalmente para passos temporais maiores, o FEM se mostrou eficaz e estável, especialmente quando combinado com o método de Euler Implícito para a discretização temporal. Além disso, a generalização do método para problemas em dimensões superiores ou com geometrias mais complexas destaca sua versatilidade e aplicabilidade em diferentes contextos da engenharia e física, apesar de que talvez fosse necessário aplicar outros métodos de discretização temporal.

== Referências ==

<references />