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.
Contenido
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 por conducción en un medio isótropo es proporcional y de sentido contrario al gradiente de temperatura en esa dirección.
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, 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].
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]. Esto da 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)=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]. Con esta igualdad deducimos que si extendemos la función [math] -x [/math] de forma ... 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]
Graficamos los 10 primeros términos de la serie en [math] t \in [0,1][/math] y podemos observar que se aproxima al estado estacionario.
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)')
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')