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].

Solución del problema

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('u(x,t)')

En la gráfica podemos observar que la solución del poblema alcanza la solución estacionaria.

2.3 Flujo en los extremos

Por la ley de Fourier sabemos que [math] q=-k u_x [/math], siendo q el flujo. Por lo que para interpretar como es el flujo en nuestro problema, vamos a graficar como es en los extremos suponiendo que [math] k=1 [/math].

right

%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)')

Al observar la gráfica en el extremo izquierdo vemos que al principio no hay flujo debido a que tanto la barra como el extremo están a 0º. Una vez avanza el tiempo el flujo va siendo cada vez más negativo lo que implica que cada vez sale más flujo por el extremo izquierdo. Esto tiene sentido ya que como la barra se va calentando y el extremo se mantiene a 0º el flujo tiende ir de donde hay mayor temperatura a donde hay una temperatura menor.

Por el contrario, en la gráfica del extremos derecho podemos observar que inicialmente....

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')