Ecuación del calor

Resumen: En este ejemplo vamos a considerar un problema parabólico dependiente de tiempo.

Problema considerado

$\Omega\subset \mathbb{R}^2$ $T>0$ $u_0: \Omega \to \mathbb{R}$ $u_0\in L^2(\Omega)$ $f: \Omega\times (0,T) \to \mathbb{R}$ $f\in C^0(0,T; L^2(\Omega))$ $u:\Omega\times(0,T) \to \mathbb{R}$ tal que

\begin{matrix} (P) & {\begin{cases} \partial_{t} u - Δ u = f & en Ω \times (0, T), \\ u = 0 & sobre \partial Ω \times (0, T), \\ u (x, y, 0) = u_{0} & en Ω . \end{cases} \end{matrix}

$u = u(x,y,t)$ $(x,y)$ $t$ $u_0= u_0(x,y)$ $t=0$ $f = f(x,y,t)$ es la fuente externa del calor.

$\eqref{eq.calor}$ , se puede combinar el método de elementos finitos (distretización en espacio) con las técnicas usadas para aproximar el problema de Cauchy (distretización en tiempo) para EDO's. Más concretamente, procedemos como sigue:

$t$ $t$ $t$ $\Omega$ .
Cada problema estacionario se resuelve mediante el método de los elementos finitos sobre una triangulación dada, como en los ejemplos anteriormente considerados.

Distretización en tiempo

$[0,T]$ $M+1$ $0=t_0 <\dots < t_M=T$ $[0,T]$ $M$ $\Delta t$ (que se llama el paso en tiempo):

\begin{matrix} Δ t = \frac{T}{M}, t_{n} = n Δ t, n = 0, \dots, M . \end{matrix}

$n = 0$ , se tiene

\begin{matrix} u (x, y, t_{0}) = u (x, y, 0) = u_{0} (x, y) . \end{matrix}

$n = 1,\dots, M$ $\eqref{eq.calor}$ $t_n$ :

\begin{matrix} (1) & {\begin{cases} \partial_{t} u (x, y, t_{n}) - Δ u (x, y, t_{n}) = f (x, y, t_{n}) & en Ω, \\ u (x, y, t_{n}) = 0 & sobre \partial Ω . \end{cases} \end{matrix}

Pongamos

\begin{matrix} (2) & u^{n} (x, y) := u (x, y, t_{n}), f^{n} (x, y) := f (x, y, t_{n}) . \end{matrix}

$u$ $t$ por un cociente incremental (método de Euler explícito), tenemos

\begin{matrix} \partial_{t} u (x, y, t_{n}) \approx \frac{u (x, y, t_{n}) - u (x, y, t_{n - 1})}{Δ t} = \frac{u^{n} (x, y) - u^{n - 1} (x, y)}{Δ t} . \end{matrix}

$\eqref{eq.calor1}$ , tenemos la ecuación aproximada

\begin{matrix} (3) & \frac{u^{n} - u^{n - 1}}{Δ t} - Δ u^{n} = f^{n} en Ω . \end{matrix}

$n$ $u^n$ $u^{n-1}$ $f^n$ es un dato.

$u^0(x,y) = u(x,y,0)$ $n=1,\dots, M$ $u^n:\Omega \to \mathbb{R}$ tal que

\begin{matrix} (4) & {\begin{cases} \frac{1}{Δ t} u^{n} - Δ u^{n} = \frac{1}{Δ t} u^{n - 1} + f^{n} & en Ω, \\ u^{n} = 0 & sobre \partial Ω . \end{cases} \end{matrix}

$n$ $\eqref{eq.calorpdt}$ es un problema que no depende del tiempo, que podemos aproximar usando el método de los elementos finitos razonando como en los ejemplos anteriores.

Discretización en espacio

Formulación variacional

Observación: todos los cálculos que se hacen en esta sección son puramente formales, es decir, se supone siempre la regularidad suficiente para que todos ellos estén justificados.

$\eqref{eq.calorpdt}$ $H_0^1(\Omega)$ .

$\eqref{eq.calorpdt}$ $n$ $v\in H_0^1(\Omega)$ $\Omega$ $v$ $\partial\Omega$ , se obtiene:

\begin{matrix} \frac{1}{Δ t} \int_{Ω} u^{n} v d x d y + \int_{Ω} \nabla u^{n} \cdot \nabla v d x d y = \frac{1}{Δ t} \int_{Ω} u^{n - 1} v d x d y + \int_{Ω} f^{n} v d x d y \forall v \in H_{0}^{1} (Ω) . \end{matrix}

$n = 1,\dots, M$ $\eqref{eq.calorpdt}$ es

\begin{matrix} (PV) & {\begin{cases} Hallar u^{n} \in H_{0}^{1} (Ω) tal que \\ \frac{1}{Δ t} \int_{Ω} u^{n} v d x d y + \int_{Ω} \nabla u^{n} \cdot \nabla v d x d y = \frac{1}{Δ t} \int_{Ω} u^{n - 1} v d x d y + \int_{Ω} f^{n} v d x d y \\ \forall v \in H_{0}^{1} (Ω) . \end{cases} \end{matrix}

Problema aproximado

$P_1$ $\partial\Omega$ $\mathcal{T}_h$ $\overline{\Omega}$ . Consideramos

\begin{matrix} W_{h} = {v_{h} | v_{h} \in C^{0} (\overset{―}{Ω}), v_{h} |_{T} \in P_{1} (T) \forall T \in T_{h}} \subset H^{1} (Ω) \end{matrix}

$\text{dim}\, V_h = M$ $M$ $\mathcal{T}_h$ . Consideramos también

\begin{matrix} V_{h} = {v_{h} | v_{h} \in V_{h}, v_{h} |_{\partial Ω} = 0} \subset H_{0}^{1} (Ω), \end{matrix}

$\mathcal{T}_h$ no $\partial\Omega$ .

$n=1,\dots,M$ $\eqref{eq.calorPV}$ es

\begin{matrix} (PV_{h}) & {\begin{cases} Hallar u_{h}^{n} \in V_{h} tal que \\ \frac{1}{Δ t} \int_{Ω} u_{h}^{n} v_{h} d x d y + \int_{Ω} \nabla u_{h}^{n} \cdot \nabla v_{h} d x d y = \frac{1}{Δ t} \int_{Ω} u_{h}^{n - 1} v_{h} d x d y + \int_{Ω} f^{n} v_{h} d x d y \\ \forall v_{h} \in V_{h} . \end{cases} \end{matrix}

Resolución con FreeFEM

$\Omega$ $[0,1]\times[0,1]$ $u_0(x,y)= \sin(\pi x)$ $f(x,y,t) =$ $T = 1$ .

$[0,1]$ $\Delta t = 0.1$ $t_n = n \Delta t$ $n = 0,\dots, 10$ .

Construcción del mallado

Declaramos primero una variable tipo mesh y después construimos una triangulación del dominio usando square.


// Triangulacion
int[int] La = [1,1,1,1];
mesh Th = square(30, 30, flags = 3, label = La);
plot(Th, wait = 1);

Definición del espacio de elementos finitos

$\mathcal{T}_h$ $W_h$ $\overline{\Omega}$ $T\in\mathcal{T}_h$ es un polinomio de grado menor o igual que 1.


// Definicion del espacio de EF
fespace Wh(Th, P1);

Definición del problema variacional

$u^n=0$ $\partial \Omega$ bloqueo $W_h$ .

$u_h^0(x,y) = u(x,y,0)$ $n=1,\dots, M$ $\eqref{eq.calorPVh2}$ $u_h^{n-1}$ es calculado en la etapa anterior, que tenemos que conservar y actualizar de una etapa a otra. Para ello, usaremos una variable auxiliar que llamaremos uold.


// Declaracion de las variables del problema (elementos de Wh)
Wh u, v, uold, u0 = sin(x), f = x*y;

Obsérvese que en cada etapa de tiempo tenemos que resolver un problema variacional del mismo tipo. Para definir el problema usando problem , usaremos el parámetro opcional init que permite resolver varios problemas con la misma forma bilineal construyendo y factorizando la matriz del sistema una sola vez:

problem Calor(u, v, init = bool)

init = false (valor por defecto) la matriz se construye y se factoriza en cada problema,
init = true, la matriz se construye y se factoriza solo en el primer problema; en los siguientes se utiliza la matriz ya factorizada.


// Definicion del problema aproximado
problem Calor(u,v, init = 1) =
      int2d(Th)((dx(u)*dx(v)+dy(u)*dy(v)))
    + int2d(Th)(dt1*u*v)
    - int2d(Th)(dt1*uold*v)
    - int2d(Th)(f*v)
    + on(1, u = 0);

Iteraciones en tiempo, resolución del problema y representación

$n=1,\dots, M$ , resolveremos el correspondiente problema aproximado.

La estructura general del bucle indexado for es la siguiente:


  for (inicializacion; condicion; incremento) 
      { bloque de instucciones }

donde

inicializacion es la inicialización del índice del bucle,
condicion es una expresión con valor lógico: si toma valor verdadero (true), se repite el bloque de instrucciones, en caso contrario, no se repite,
incremento indica cómo cambia el índice del bucle de una iteración a otra,
bloque de instucciones es el conjunto de instrucciones a realizar.

Resolvemos el problema y representamos la solución en cada etapa de tiempo.


// Datos
int M = 10;
real T = 1., dt = T/M, dt1 = 1./dt;

// Iteraciones en tiempo 
u = u0;
for (int n = 1; n <= M; ++n){
  uold = u;
  Calor;
  plot(u, fill = true, dim = 3, wait = 1);
}

Obsérvese que el operador compuesto ++nn $t$ $t$ no aparece en el bucle.

$T$ .

Código completo del programa


// Datos 
real T = 1.;
int M = 10; 
real dt = T/M, dt1 = 1./dt;

// Triangulacion
int[int] La = [1,1,1,1];
mesh Th = square(30, 30, flags = 3, label = La);
plot(Th, wait = 1, ps = "MalladoCalor.eps"); 

// Espacio elementos finitos
fespace Wh(Th, P1);
Wh u, v, uold, u0 = sin(x), f = x*y;

// Definicion del problema
problem Calor(u,v, init = 1) =
        int2d(Th)((dx(u)*dx(v)+dy(u)*dy(v)))
      + int2d(Th)(dt1*u*v)
      - int2d(Th)(dt1*uold*v)
      - int2d(Th)(f*v)
      + on(1, u = 0);

// Iteraciones en tiempo
u = u0;
for (int n = 1; n <= M; ++n){
  uold = u;
  Calor;
  plot(u, fill = true, dim = 3, wait = 1);
}

Ejercicios

$\eqref{eq.calor}$ realizando la resolución manual con varf.
$\Omega = [0,1]^2$ $T = 3$ $\mu = 0.01$ , resolver el siguiente problema

\begin{matrix} (5) & {\begin{cases} \partial_{t} u - μ Δ u = x^{4} - 12 μ t x^{2} & en Ω \times (0, T), \\ u = t x^{4} & sobre \partial Ω \times (0, T), \\ u (x, y, 0) = 0 & en Ω . \end{cases} \end{matrix}

$\eqref{eq.calor3}$ usando la resolución manual con varf.
$\Omega = [0,L_x]\times [0,L_y]$ $L_x = 6$ $L_y = 1$ $T = 5$ $\alpha = 0.25$ $u_e = 25$ $a_0 = 10$ $a_1 = 90/Lx$ , resolver el siguiente problema

\begin{matrix} (6) & {\begin{cases} \partial_{t} u - \nabla \cdot (κ \nabla u) = 0 & en Ω \times (0, T), \\ κ \partial_{n} u + α (u - u_{e}) = 0 & sobre (Γ_{1} \cup Γ_{3}) \times (0, T), \\ u = a_{0} + L_{x} a_{1} & sobre Γ_{2} \times (0, T), \\ u = a_{0} & sobre Γ_{4} \times (0, T), \\ u (x, y, 0) = a_{0} + a_{1} x & en Ω . \end{cases} \end{matrix}

$\kappa$ $\Gamma_i$ $i = 1,2,3,4$ vienen dados como sigue:

\begin{matrix} \begin{matrix} κ = {\begin{cases} 2 & para y < L_{y} / 2, x \in [0, L_{x}], \\ 0.2 & si no, \end{cases} \end{matrix} \end{matrix}

$\eqref{eq.calor4}$ $u_e = 25 + (T_e-25)t/T$ $t$ $T_e = 200$ .
$\eqref{eq.calor4}$ usando la resolución manual con varf.

Anna Doubova - Rosa Echevarría - Dpto. EDAN - Universidad de Sevilla