Ecuación del calor (Raúl, Sofía, Jaime)

Trabajo realizado por estudiantes
Título	Ecuación del calor
Asignatura	EDP
Curso	2023-24
Autores	Raúl Ortega Sofía Gómez Jaime Sáenz de Miera
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

En este artículo vamos a estudiar la ecuación del calor en una dimensión mediante el estudio de una varilla metálica de longitud 1 y que se encuentra aislada por su superficie lateral, es decir, la conducción de calor solo se produce en la dirección longitudinal. Además, vamos a ver cómo varía la solución en función de la variación de distintos parámatros como las condiciones iniciales o el coeficiente de difusión.

1 Preliminares

Antes de empezar a desarrollar el problema necesitamos conocer algunos conceptos que nos van a permitir modelizarlo y sacar sus ecuaciones.

• La ley de Fourier establece que el flujo de transferencia de calor es proporcional y de sentido contrario al gradiente de temperatura en esa dirección.

[math] q=-k\nabla u [/math]

Siendo k la conductividad térmica.

• Conservación de la energía.

2 Problema original

En este primer apartado, vamos a suponer que la temperatura inicial de la varilla es 0ºC y que en el extremo izquierdo se consigue mantener la temperatura a 0ºC mientras que en el derecho la temperatura es siempre de 1ºC. Teniendo esto en cuenta y los conceptos desarrollados en el apartado anterior, el sistema de ecuaciones que representa el problema es:

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

2.1 Homogeneización

Antes de empezar con la resolución del problema vamos a homogeneizar las condiciones de contorno. Para hacer esto vamos a restarle al problema la solución estacionaria, es decir, la que permanece invariante respecto al tiempo. Esto implica que [math] u_t\sim 0 [/math] y que [math] u(x,t)\sim v(x) [/math]. Por lo que el problema queda de la siguiente forma:

[math] \begin{cases} v_{xx}=0 & x\in(0,1)\\ v(0)=0\\ v(1)=1\\ \end{cases} [/math]

La solución de este problema, y por tanto la solución estacionaria es: [math] v(x)=x [/math].

Solución estacionaria

v=@(x,t) x;
colormap jet
fsurf(v, [0 1 0 1], 'EdgeColor','interp')
xlabel('x')
ylabel('t')
zlabel('v(x)')

Una vez tenemos la solución estacionaria, tomamos [math] w(x,t)=u(x,t)-v(x) [/math] y vemos como queda el problema homogeneizado para [math] w(x,t)[/math]:

[math] \begin{cases} w_t-w_{xx}=0 & x\in(0,1),t\gt0\\ w(0,t)=0&t\gt0\\ w(1,t)=0&t\gt0\\ w(x,0)=-x&x\in(0,1) \end{cases} [/math]

2.2 Resolución

Una vez homogeneizado el sistema podemos empezar a buscar la solución por el método de separación de variables, en el cual vamos a suponer que la solución es de la forma: [math] w(x,t)=X(x)T(t) [/math]. Como [math]w(x,t)[/math] es solución obtenemos que se tiene que cumplir [math]\frac{X''(x)}{X(x)}=\frac{T'(t)}{T(t)}=\lambda_k[/math] dando lugar a la resolución de dos problemas por separado, uno que depende de [math] x [/math] y otro que depende de [math] t[/math]. Al solucionar ambos problemas obtenemos el siguiente resultado: [math] w(x,t)=\sum_{k=1}^\infty c_k\sin(k\pi x)e^{-k^2\pi^2 t} [/math].

Para que [math] w(x,t) [/math] sea solución de nuestro problema, solo falta ver que cumpla la condición incial, es decir, que [math] w(x,0)=c_k\sin(k\pi x)=-x [/math]. De esta igualdad sacamos que si extendemos la función [math] -x [/math] de forma impar los coeficientes [math] c_k[/math] van a ser los coeficiente de la serie de Fourier. Por tanto, la solución queda de la siguiente manera: [math] w(x,t)=\sum_{k=1}^\infty\frac{2(-1)^k}{k\pi}\sin(k\pi x)e^{-k^2\pi^2 t}[/math]

Esta función es la solución del problema homogeneizado, por lo que para obtener la solución del problema original tenemos que deshacer el cambio de variable. La solución es:

[math] u(x,t)=\sum_{k=1}^\infty\frac{2(-1)^k}{k\pi}\sin(k\pi x)e^{-k^2\pi^2 t}+x [/math]

Ahora vamos a graficar los 10 primeros términos de la serie en [math] t \in [0,1][/math].

clear all
%Definimos las funciones
wi=@(x,t,k) ((2*(-1)^(k))/(k*pi)).*sin(k.*pi.*x).*exp(-k^2.*pi^2.*t);
v=@(x) x;

%Sumamos los 10 primeros términos
w=@(x,t) 0;
for n=1:10
    w=@(x,t) w(x,t) + wi(x,t,n);
end
u=@(x,t) w(x,t)+v(x);

% Gráfica en 3D
close all
colormap jet
fsurf(u,[0 1 0 1], 'EdgeColor','interp')
xlabel('x')
ylabel('t')
zlabel('w(x,t)')

En la gráfica podemos observar que la solución se aproxima a la solución estacionaria.

2.3 Flujo en los extremos

Por la ley de Fick sabemos que el flujo es proporcional a menos la derivada de la temperatura. Por lo que para interpretar como es el flujo en n uestro problema, vamos a graficar como es en los extremos.

%Definimos las derivadas 
dwi=@(x,t,k) (2*(-1)^(k+1)).*cos(k.*pi.*x).*exp(-k^2.*pi^2.*t);
dv=@(x) 1;

%Sumamos los 10 primeros términos
dw0=@(t) 0;
dw1=@(t) 0;
for n=1:10
    dw0=@(t) dw0(t) + dwi(0,t,n);
    dw1=@(t) dw1(t) + dwi(1,t,n);
end
du0=@(t) dw0(t)-dv(0);
du1=@(t) dw1(t)-dv(1);

% Gráfica 
close all
t=linspace(0,1,100);
subplot(2,1,1)
plot(t,du0(t))
xlabel('t')
ylabel('u´(0,t)')
subplot(2,1,2)
plot(t,du1(t))
xlabel('t')
ylabel('u´(1,t)')

3 Variación coeficiente de conductividad

En este apartado vamos a ver como varia la solución cuando variamos el coeficiente de conductividad térmica de 1 a 1/2.

[math] \begin{cases} u_t-\frac{1}{2}u_{xx}=0 & x\in(0,1),t\gt0\\ u(0,t)=0&t\gt0\\ u(1,t)=1&t\gt0\\ u(x,0)=0&x\in(0,1) \end{cases} [/math]

Diferencia con la estacionaria

%Definimos las funciones
wi1=@(x,t,k) ((2*(-1)^(k))/(k*pi)).*sin(k.*pi.*x).*exp(-k^2.*pi^2.*t);
wi2=@(x,t,k) ((2*(-1)^(k))/(k*pi)).*sin(k.*pi.*x).*exp(-k^2.*pi^2.*t./2);
v=@(x) x;

%Sumamos los 10 primeros términos
w1=@(x,t) 0;
w2=@(x,t) 0;
for n=1:10
    w1=@(x,t) w1(x,t) + wi1(x,t,n);
    w2=@(x,t) w2(x,t) + wi2(x,t,n);
end

u1=@(x,t) w1(x,t)+v(x);
u2=@(x,t) w2(x,t)+v(x);

%Hacemos la diferencia con la solución estacionaria
dif1=@(t) u1(0.5,t)-v(0.5);
dif2=@(t) u2(0.5,t)-v(0.5);

% Gráfica 
close all
t=linspace(0,1,100);
subplot(2,1,1)
plot(t,dif1(t))
xlabel('t')
ylabel('w(0.5,t)-v(0.5)')
subplot(2,1,2)
plot(t,dif2(t))
xlabel('t')
ylabel('w(0.5,t)-v(0.5)')

4 Variación de las condiciones de contorno y las iniciales

[math] \begin{cases} u_t-u_{xx}=0 & x\in(0,1),t\gt0\\ u(0,t)=0&t\gt0\\ u(1,t)=0&t\gt0\\ u(x,0)=max\{0,1-4|x-1/2|\}&x\in(0,1) \end{cases} [/math]

u0=@(x) max(0,1-4*abs(x-1/2));
uk=@(x,t,k) sin(k.*pi.*x).*exp(-k.^2*pi.^2.*t);
close all

a=0; b=1;                    
x=linspace(a,b,100);
t=linspace(a,b,100);
N=length(x)-1;
h=(b-a)/N;
w=ones(N+1,1);              
w(1)=1/2; w(N+1)=1/2;
f_n=zeros(length(x),length(t));
for i=1:length(t)
    for k=1:100
        g=u0(x)*sin(k*pi*x)'; 
        a_k=2*h*w'*g;
        f_n(:,i)= f_n(:,i)'+a_k.*uk(x,t(i),k);
    end
end
colormap jet 
surf(x,t,f_n,'EdgeColor','interp')