Ver el código de «Ecuación de calor (ADPP)»

{{ TrabajoED | Ecuación de calor (ADPP). | [[:Categoría:EDP|EDP]]|[[:Categoría:EDP24/25|2024-25]] | Alex Heredero Santamaría, Diego Moñino Vizmanos, Pablo Vidal Nacle, Pablo Maestro Fernández}}

=. Introducción=
La ecuación del calor clásica, formulada por Joseph Fourier, describe la difusión térmica en medios continuos y ha sido ampliamente utilizada en la física y la ingeniería. Sin embargo, un aspecto problemático de esta ecuación es que predice una velocidad de propagación infinita para las perturbaciones térmicas. Es decir, si se genera una variación de temperatura en un punto del sistema, esta afectaría instantáneamente a todo el dominio, lo cual es físicamente irreal.

Este problema, conocido como la paradoja de la velocidad infinita de propagación, entra en conflicto con la relatividad especial, ya que ninguna perturbación física debería propagarse más rápido que la velocidad de la luz. Para corregir esta deficiencia, se han desarrollado modelos alternativos, entre ellos la ecuación de Cattaneo-Vernotte, que introduce un término de relajación en la ecuación del calor, limitando así la velocidad de propagación de las señales térmicas.

Resolver la ecuación de Cattaneo-Vernotte es fundamental para modelar procesos térmicos en situaciones donde la ecuación clásica de Fourier no es adecuada. Esto ocurre, por ejemplo, en la propagación del calor en materiales con estructura micro o nanoescala, en plasmas, y en medios donde la conducción térmica no es instantánea. Además, su estudio permite una mejor comprensión de los procesos de difusión modificados y su relación con principios relativistas.

Desde un punto de vista práctico, el análisis de esta ecuación puede mejorar modelos de transferencia de calor en la industria, la ingeniería de materiales y la física de semiconductores. Nuestros objetivos en este artículo serán introducir la ecuación de Cattaneo-Vernotte como una alternativa físicamente más realista, resolverla en un caso partícular, interpretar la solución obtenida y compararla con la ecuación clásica, y finalmente discutir las implicaciones físicas del modelo.

=. Problema concreto=
Planteamos el siguiente problema: Se considera una varilla metálica que ocupa el intervalo [0, 1] y que se encuentra aislada por su superficie lateral, de manera que la conducción de calor sólo se produce en la
dirección longitudinal. En el extremo derecho
se consigue mantener la temperatura a 10°C mientras que en el izquierdo la temperatura
es siempre de 1°C. Además, la temperatura en el instante inicial viene dada por la función <math> u_0 (x)= 10-10 \cdot 1_{[1/3,2/3]}(x)</math> . Así, el sistema que modeliza este problema es el siguiente:

<center><math> \begin{cases} 
      u_t-u_{xx} = 0 & x \in [0,1], t>0 \\
      u(0,t)=10 & t>0 \\
      u(1,t)=1 & t>0 \\
      u(x,0)=u_0(x) & x \in [0,1]
\end{cases} </math></center>

Su solución estacionaria es:

<center><math> v(x)=9x+1 </math></center>

y la solución general es:

<center><math> u(x,t)= 9x+1+  \sum_{n=1}^{\infty}b_n sen(n \pi x)e^{n^2\pi ^2 t}</math>,</center>

donde <math> b_n </math> viene dado por:

<center><math> b_n = \int_{-1}^{1}g(x)sen(n\pi x)dx </math>,</center>
con g(x):

<center><math> \begin{cases} 
      -9x-9 & x \in (-1,2/3)\cup (-1/3,0) \\
      -9x+9 & x \in (0,1/3)\cup (2/3,1) \\
      -9x+1 & x \in [-2/3,-1/3] \\
      -9x-1 & x \in [1/3,2/3]
\end{cases} </math></center>

Si visualizamos estas soluciones con <math> t \in [0,1] </math> y <math> x \in [0,1] </math>, representándola a través de matlab:

[[Archivo:sol2.jpeg|400px|thumb|right|Solución de la ecuación del calor]]

[[Archivo:sol11.jpeg|400px|thumb|right|Solución de la ecuación del calor]]

[[Archivo:Animación solución.gif|400px|thumb|right|Convergencia a la solución estacionaria]] 

<syntaxhighlight lang="matlab">

syms n
syms x
syms t
g=sin(n*pi*x);
f1=(-9*x-9)*g;
f2=(-9*x+9)*g;
f3=(-9*x+1)*g;
f4=(-9*x-1)*g;

F1=int(f1,-1,-2/3);
F2=int(f3,-2/3,-1/3);
F3=int(f1,-1/3,0);
F4=int(f2,0,1/3);
F5=int(f4,1/3,2/3);
F6=int(f2,2/3,1);

F(n)=F1+F2+F3+F4+F5+F6;

sol=9*x+1;
k=20; %numero de elementos
for i=1:k
    sol=sol+F(i)*sin(i*pi*x)*exp(-(i^2)*pi^2*t);
end
u(x,t)=sol;
[X, T] = meshgrid(linspace(0,1,50), linspace(0,1,50));
Z=u(X,T);
Z=double(Z);
figure
surf(X, T, Z) 
xlabel('x'), ylabel('t'), zlabel('u(x,t)')
title('Gráfico de la solución')
colorbar
shading interp 

% Definir parámetros
xv = linspace(0,1,100); % Dominio de x (100 puntos entre 0 y 1)
t_values = linspace(0,0.1,60); % Valores de t para cada frame (60 frames)

% Crear la figura
fig = figure;
v = VideoWriter('animacion_ftejercicio2.mp4', 'MPEG-4'); % Guardar animación en video
v.FrameRate = 5; % FPS
open(v);

for t = t_values
    y = u(xv,t); % Evaluar la función en x para el tiempo t
    plot(xv, y, 'b', 'LineWidth', 2); % Graficar en 1D
    ylim([0, 12]); % Mantener el mismo rango en y
    xlabel('x');
    ylabel('u(x,t)');
    title(sprintf('t = %.2f', t)); % Mostrar el tiempo actual
    grid on;
    drawnow;
    
    % Capturar frame para el video
    frame = getframe(fig);
    writeVideo(v, frame);
end

close(v); % Cerrar archivo de video
close(fig); % Cerrar la figura

disp('Animación guardada como animacion_ft.mp4');
</syntaxhighlight>

Observamos que la solución u comienza en la temperatura correcta (salvo por oscilaciones) y para un tiempo suficientemente pequeño se aproxima considerablemente a la solución de equilibrio. No obstante, en el siguiente apartado se comentará un detalle de la solución que se debe estudiar en mayor profundidad.

=. Paradoja de la velocidad de propagación=

[[Archivo:Paradoja.jpeg|400px|thumb|derecha|Solución para t=0]]

Observamos que, en nuestra solución, para <math>t>0</math> la solución es estrictamente positiva. Sin embargo, en <math>t=0</math> observamos que para el intervalo <math>  x \in [1/3,2/3]</math> la temperatura es cero. Esto contradice la teoría de la relatividad de Einstein, ya que podemos tomar un valor <math>t>0</math> arbitrariamente pequeño de manera que la temperatura sea no nula, lo cual no es posible bajo la teoría de la relatividad. Este fenómeno se conoce como la paradoja de la velocidad de propagación del calor. 

==. Ecuación de Cattaneo-Vernotte==
Para solucionar esta aparente falla en nuestro modelo, Cattaneo propuso una solución modificando la Ley de Fourier de la siguiente manera:

<center><math>q=-k\nabla u \rightarrow q+ \tau \frac{\partial q}{\partial t} =k\nabla u</math></center>

Añadimos el factor <math> \tau \frac{\partial q}{\partial t} </math> a la Ley de Fourier, donde <math>\tau</math> se conoce como el término de relajación térmica. De esta manera, se modela un retardo en la respuesta de flujo de calor a los cambios de la temperatura. Si ahora aplicamos la conservación de la energía, nos queda la conocida ecuación hiperbólica del calor:

<center><math> \frac{\partial u}{\partial t} + \tau \frac{\partial^{2}u}{\partial t^{2}} = \alpha \Delta u </math></center>

donde <math> \alpha </math> es la difusividad térmica. Para nuestro problema, la ecuación queda de la siguiente forma:

<center><math> \begin{cases} 
      \tau u_{tt}+u_t-u_{xx} = 0 & x \in [0,1], t>0 \\
      u(0,t)=10 & t>0 \\
      u(1,t)=1 & t>0 \\
      u(x,0)=u_0(x) & x \in [0,1] \\
      u_t(x,0)=0 & x \in [0,1]
\end{cases} </math></center>

Hemos introducido en las condiciones iniciales el factor <math> u_t(x,0) </math>, que se conoce como la velocidad térmica inicial. En nuestro caso hemos establecido esta velocidad nula ya que se asume que en el estado inicial el sistema está en equilibrio.

Esta nueva ecuación garantiza que la velocidad de propagación del calor se produzca con velocidad finita, por lo que se utiliza en modelos donde el efecto de propagación instantánea tiene importancia. 

Homogeneizamos el problema con el cambio de variable <math> w(t,x)=u(t,x) -v(x)</math>. Si resolvemos esta ecuación por separación de variables, tenemos las siguientes soluciones para la función espacial <math>X(x)</math> y temporal <math> T(t)</math>:

La ecuación espacial con estas condiciones nos proporciona la colección de soluciones <math> X_{n}(x)=A_{n}sen(n\pi x) </math>

Para la ecuación temporal, tenemos la solución <math> T_{n}(t)=e^{\frac{-t}{2\tau}}(c_{1}cos(\omega t) + c_{2}sen((\omega t)) </math>, con <math> \omega = \frac{\sqrt{4\tau n^{2}\pi^{2}-1}}{2\tau} </math>

Por tanto, la solución de la ecuación de Cattaneo-Vernotte en nuestro caso es:

<center><math> u(x,t) = 1+9x+ \sum_{n=1}^{\infty} A_{n}e^{\frac{-t}{2\tau}}(c_{1}cos(\omega t) + c_{2}sen((\omega t)) \cdot sen(n\pi x) </math></center>

Ahora, utilizando las condiciones frontera, obtenemos <math> c_{2}=0</math>  y <math> A_{n}c_{1}=b_{n} </math>, quedándonos la solución

<center> <math>u(x,t) = 1+9x+ \sum_{n=1}^{\infty} b_{n}e^{\frac{-t}{2\tau}}cos(\omega t)  \cdot sen(n\pi x) </math></center>
 
donde <math> b_{n} </math> son los mismos coeficientes que en la solución de la ecuación original.

[[Archivo:Solucion_cattaneo1.jpg|400px|thumb|right|Solución de la ecuación de Cattaneo-Vernotte]]

[[Archivo:Solucion_cattaneo2.jpg|400px|thumb|right|Solución de la ecuación de Cattaneo-Vernotte]]

[[Archivo:Animacioncattaneo.gif|400px|thumb|right|Convergencia a la solución estacionaria]] 

<syntaxhighlight lang="matlab">

syms n
syms x
syms t
g=sin(n*pi*x);
f1=(-9*x-9)*g;
f2=(-9*x+9)*g;
f3=(-9*x+1)*g;
f4=(-9*x-1)*g;

F1=int(f1,-1,-2/3);
F2=int(f3,-2/3,-1/3);
F3=int(f1,-1/3,0);
F4=int(f2,0,1/3);
F5=int(f4,1/3,2/3);
F6=int(f2,2/3,1);

F(n)=F1+F2+F3+F4+F5+F6;

sol=9*x+1;
k=40; %numero de elementos
for i=1:k
    sol=sol+F(i)*sin(i*pi*x)*exp(-t/(2*tau))*cos(t*(sqrt(4*tau*i^2*pi^2-1))/(2*tau));
end
u(x,t)=sol;

[X, T] = meshgrid(linspace(0,1,50), linspace(0,1,50));
Z=u(X,T);
Z=double(Z);
figure
surf(X, T, Z) 
xlabel('x'), ylabel('t'), zlabel('u(x,t)')
title('Gráfico de la solución')
colorbar
shading interp 

% Definir parámetros
xv = linspace(0,1,100); % Dominio de x (100 puntos entre 0 y 1)
t_values = linspace(0,1,60); % Valores de t para cada frame (60 frames)

% Crear la figura
fig = figure;
v = VideoWriter('animacion_ftejercicio3.mp4', 'MPEG-4'); % Guardar animación en video
v.FrameRate = 5; % FPS
open(v);

for t = t_values
    y = u(xv,t); % Evaluar la función en x para el tiempo t
    plot(xv, y, 'b', 'LineWidth', 2); % Graficar en 1D
    ylim([0, 12]); % Mantener el mismo rango en y
    xlabel('x');
    ylabel('u(x,t)');
    title(sprintf('t = %.2f', t)); % Mostrar el tiempo actual
    grid on;
    drawnow;
    
    % Capturar frame para el video
    frame = getframe(fig);
    writeVideo(v, frame);
end

close(v); % Cerrar archivo de video
close(fig); % Cerrar la figura

disp('Animación guardada como animacion_ft.mp4');

</syntaxhighlight>



Como se puede apreciar en las gráficas, la convergencia a la solución estacionaria es suave y no abrupta como en el problema original, logrando una modelización más fiel de la transferencia de calor tras una perturbación térmica, como queríamos comprobar al inicio del artículo.

=. Referencias =

* Amin Moosaie. Non-Fourier heat conduction in a finite medium with insulated boundaries and arbitrary initial conditions ([https://www.sciencedirect.com/science/article/pii/S0735193307001637 [1<nowiki>]</nowiki>]).

* Francisco R. Villatoro. La ciencia de la mula Francis: La velocidad de la propagación del calor, entre la paradoja y la entropía([https://francis.naukas.com/2008/10/22/la-velocidad-de-la-propagacion-del-calor-entre-la-paradoja-y-la-entropia/ [2<nowiki>]</nowiki>])

* Marc Calvo Schwarzwälder. Non-Fourier Heat Conduction. The Maxwell-Cattaneo Equations ([https://upcommons.upc.edu/bitstream/handle/2117/78480/memoria.pdf [3<nowiki>]</nowiki>]).

[[Categoría:EDP]]
[[Categoría:EDP24/25]]