Ecuación de Calor en varilla(Grupo 16-C)

Trabajo realizado por estudiantes
Título	Ecuación del Calor en varilla. Grupo 16-C
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Laura Ramos Sangrós Luis Moreno López Álvaro Martínez de Andrés Gonzalo Olmos de la Cruz
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción y planteamiento del sistema

Se nos presenta una varilla delgada, homogénea, de sección constante y longitud L ocupando el intervalo x € (0,L) con L = 4 m. Las condiciones físicas de la varilla son densidad=1 y conductividad térmica(k)=2. En los extremos de la barra están colocados objetos a una temperatura constante de 0 y 5 grados respectivamente, e inicialmente la temperatura de la varilla viene dada por la función: u0 = e^(-6(x-2)^2) + (5-5x/4) Suponiendo que la temperatura de la varilla u(x,t) satisface la ecuación de calor ut-2uxx = 0, El problema (P) con su sistema de ecuaciones que satisface u(x,t) será:

    [math] \left \{\begin{matrix}    Ut – 2Uxx = 0   ;      t\gt0, x € (0,4) \\  U(0,t) = 5 ; U(4,t) = 0  ;      t\gt0 \\ U(x,0) = e^(-6(x-2)^2) + (5-5x/4)  ;  t\gt0, x € (0,4) \end{matrix} \right . [/math]

2 Métodos de resolución

En segundo lugar se nos pide resolverlo por el método de diferencias finitas con un tamaño de paso y un intervalo de tiempo de 0.1 por el método del trapecio.

El código en Matlab y su correspondiente gráfico será:

%Trabajo ejercicio 2 
clear all , clc , clf

%Resolverlo por trapecio
%datos iniciales
a=0; b=4; h=0.1;
x=a:h:b;
N=(b-a)/h; %el vector x tiene n+1 elementos
q=2; %término que acompaña a Uxx
%nodos interiores
xx=x(2:N);
xx=xx';
%cond contorno
ua=5; ub=0;
%cond inicial
u0=exp(-6*(xx-2).^2)+5-5*xx/4;
% matriz K
K=q/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
% término F 
F=0*xx;
F(1)=F(1)+q*ua/h^2;
F(end)=F(end)+q*ub/h^2;
% discretización en t
t0=0;tM=1;
k=0.1; % paso temporal.
t=t0:k:tM;
M=length(t)-1; %número de subintervalos temporales
sol(:,1)=u0;
for i=1:M
%Trapecio
sol(:,1+i)=(eye(size(K))+(k/2)*K)\(sol(:,i)+(k/2)*(-K*sol(:,i)+2*F));
end
%Añado frontera
Ua=q*ua*ones(1,length(t));
Ub=q*ub*ones(1,length(t));
sol=[Ua;sol;Ub];

figure(1)
[Mt,Mx]=meshgrid(t,x); 
surf(Mx,Mt,sol)
xlabel('Longitud');
ylabel('Tiempo');
zlabel('Temperatura');

A continuación se nos pide ver cómo se comporta la temperatura en el punto medio de la barra(X=2). Para ello congelaremos en un gráfico 2D (temperatura/tiempo) en X=2 a lo largo del tiempo.Usando el método de diferencias finitas y trapecio,el código será similar al apartado anterior salvo unos cambios al final:

%Trabajo ejercicio 3 
clear all , clc , clf

%Resolverlo por trapecio

%datos iniciales
a=0; b=4; h=0.1;
x=a:h:b;
N=(b-a)/h; %el vector x tiene n+1 elementos
q=2; %término que acompaña a Uxx
%nodos interiores
xx=x(2:N);
xx=xx';
%cond contorno
ua=5; ub=0;
%cond inicial
u0=exp(-6*(xx-2).^2)+5-5*xx/4;
% matriz K
K=q/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
% término F 
F=0*xx;
F(1)=F(1)+q*ua/h^2;
F(end)=F(end)+q*ub/h^2;
% discretización en t
t0=0;tM=1;
k=0.1; % paso temporal.
t=t0:k:tM;
M=length(t)-1; %número de subintervalos temporales
sol(:,1)=u0;
for i=1:M
%Trapecio
sol(:,1+i)=(eye(size(K))+(k/2)*K)\(sol(:,i)+(k/2)*(-K*sol(:,i)+2*F));
end
%Añado frontera
Ua=q*ua*ones(1,length(t));
Ub=q*ub*ones(1,length(t));
sol=[Ua;sol;Ub];

figure(2)
%Comportamiento de la temperatura en el punto medio de la barra
hold on
long=(b-a)/2;
fila=round((long-a)/h+1);
plot(t,sol(fila,:))
xlabel('Tiempo');
ylabel('Temperatura');
hold off

Es el punto donde más variación de temperatura se observa.

Terminamos la primera parte resolviendo este problema (P) con las mismas condiciones que las anteriores y por el método de diferencias finitas, mediante la resolución de los métodos de Euler eplícito, Euler implícito y Runge-Kutta de orden 4, para más tarde comparar sus gráficas como se pide en el enunciado. A primera vista las gráficas obtenidas serán:

%Trabajo ejercicio 4 
clear all , clc , clf

%datos iniciales
a=0; b=4; h=0.1;
x=a:h:b;
N=(b-a)/h; %el vector x tiene n+1 elementos
q=2; %término que acompaña a Uxx
%nodos interiores
xx=x(2:N);
xx=xx';
%cond contorno
ua=5; ub=0;
%cond inicial
u0=exp(-6*(xx-2).^2)+5-5*xx/4;
% matriz K
K=q/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
% término F
F=0*xx;
F(1)=F(1)+q*ua/h^2;
F(end)=F(end)+q*ub/h^2;
% discretización en t
t0=0;tM=1;
k=0.1 % paso temporal.
t=t0:k:tM;
M=length(t)-1; %número de subintervalos temporales
sol(:,1)=u0;
sol1(:,1)=u0;
sol2(:,1)=u0;
for i=1:M
%Euler explicito
sol(:,1+i)=sol(:,i)+k*(-K*sol(:,i)+F);
%Euler implicito
sol1(:,1+i)=(eye(size(K))+k*K)\(sol1(:,i)+k*F);
%Kutta
k1=-K*sol2(:,i)+F;
k2=-K*(sol2(:,i)+(1/2)*k1*k)+F;
k3=-K*(sol2(:,i)+(1/2)*k2*k)+F;
k4=-K*(sol2(:,i)+k3*k)+F;
sol2(:,i+1)=sol2(:,i)+(k/6)*(k1+k2*2+k3*2+k4);
end
%incluimos cond contorno
Ua=q*ua*ones(1,length(t));
Ub=q*ub*ones(1,length(t));
sol=[Ua;sol;Ub];
sol1=[Ua;sol1;Ub];
sol2=[Ua;sol2;Ub];
%grafico para comparar
figure(1)
[Mt,Mx]=meshgrid(t,x); 
mesh(Mx,Mt,sol) 
title('Euler explícito')
xlabel('Longitud');
ylabel('Tiempo');
zlabel('Temperatura');
figure(2)
mesh(Mx,Mt,sol1)
title(' Euler implícito')
xlabel('Longitud');
ylabel('Tiempo');
zlabel('Temperatura');
figure(3)
mesh(Mx,Mt,sol2)
title('Runge-Kutta')
xlabel('Longitud');
ylabel('Tiempo');
zlabel('Temperatura');

Como se puede observar, la solución no tiene nada que ver con la realidad,solo concuerda en el caso de Euler Implícito,lo cual induce a pensar que hay un error en el desarrollo de la modelización. Si se revisa el programa, se advierte que la distorsión se debe a que con las condiciones Δx y Δt dadas el problema no cumple la condición de estabilidad necesaria: Δt/(Δx)^2 Y en este caso no se cumple por lo tanto basta con buscar un Δt que haga que cumpla: Δt/(Δx)^2 = k/h^2 = 1/2Δt , con lo que Δt = 0.0025. El nuevo código sería:

%Trabajo ejercicio 4 
clear all , clc , clf

%datos iniciales
a=0; b=4; h=0.1;
x=a:h:b;
N=(b-a)/h; %el vector x tiene n+1 elementos
q=2; %término que acompaña a Uxx
%nodos interiores
xx=x(2:N);
xx=xx';
%cond contorno
ua=5; ub=0;
%cond inicial
u0=exp(-6*(xx-2).^2)+5-5*xx/4;
% matriz K
K=q/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
% término F
F=0*xx;
F(1)=F(1)+q*ua/h^2;
F(end)=F(end)+q*ub/h^2;
% discretización en t
t0=0;tM=1;
k=0.0025; % paso temporal.
t=t0:k:tM;
M=length(t)-1; %número de subintervalos temporales
sol(:,1)=u0;
sol1(:,1)=u0;
sol2(:,1)=u0;
for i=1:M
%Euler explicito
sol(:,1+i)=sol(:,i)+k*(-K*sol(:,i)+F);
%Euler implicito
sol1(:,1+i)=(eye(size(K))+k*K)\(sol1(:,i)+k*F);
%Kutta
k1=-K*sol2(:,i)+F;
k2=-K*(sol2(:,i)+(1/2)*k1*k)+F;
k3=-K*(sol2(:,i)+(1/2)*k2*k)+F;
k4=-K*(sol2(:,i)+k3*k)+F;
sol2(:,i+1)=sol2(:,i)+(k/6)*(k1+k2*2+k3*2+k4);
end
%incluimos cond contorno
Ua=q*ua*ones(1,length(t));
Ub=q*ub*ones(1,length(t));
sol=[Ua;sol;Ub];
sol1=[Ua;sol1;Ub];
sol2=[Ua;sol2;Ub];
%grafico para comparar
figure(1)
[Mt,Mx]=meshgrid(t,x); 
mesh(Mx,Mt,sol) 
xlabel('Longitud');
ylabel('Tiempo');
zlabel('Temperatura');
figure(2)
mesh(Mx,Mt,sol1)
xlabel('Longitud');
ylabel('Tiempo');
zlabel('Temperatura');
figure(3)
mesh(Mx,Mt,sol2)
title('Kutta');
xlabel('Longitud');
ylabel('Tiempo');
zlabel('Temperatura');

3 Variación de la temperatura en el tiempo

Otra conclusión interesante que se extrae de este problema es que para tiempos grandes la temperatura apenas cambia en el tiempo. En este caso el término ut se puede despreciar y la temperatura toma un valor estacionario en cada punto de la varilla. Esta ecuación estacionaria será la siguiente al despreciar este término: 2Uxx = 0 Para demostrarlo gráficamente, basta con utilizar el método de Euler implícito, tanteando con distintos tiempos. El código es el mismo que antes pero al variar el valor de t(t=0,1,2,10), se obtienen las siguientes gráficas.

Para saber cuándo empieza a ser estacionaria la temperatura probamos a editar un gráfico en 2D con t=1 y efectivamente a partir de este instante la solución parece ser estable,es decir,estacionaria. Lo demostramos en estos gráficos,que en el primero se observa como cambia en los puntos centrales de la barra con un paso de 0.1 hasta que se convierte en estacionaria, osea en t=1. De ahí que las gráficas se vayan acercando a una única hasta llegar al punto de estar superpuestas en una misma e única, que será la estacionaria y se cumple a partir de t=1.

Este segundo gráfico mostrará cómo se comporta la función estacionaria a partir de ese instante t=1 para cualquier tiempo más grande, que seguirá el mismo patrón para cualquier tiempo mayor de t=1.

La función a la que se parece la solución u(x,t) para tiempos grandes es: ye=5-(5/4)*x, resuelta analíticamente a partir de la ecuación estacionaria antes dicha (2Uxx= 0). La función gráficamente es(es idéntica para cualquier tiempo):

Así la distancia entre la solución transitoria y estacionaria será la resta entre ambas soluciones --> err=abs(ye-sol1) que en matlab dará una tabla con las distancias en cada instante t por paso(0.1).

4 Extremo con pieza aislante

Continuando con la manipulación, al colocar en el extremo derecho una pieza aislante (en lugar de un objeto con temperatura constante 5) se provoca que no haya pérdida de calor en ese extremo, es decir, que el flujo de temperatura sea nulo. Esto provoca un cambio en las condiciones de frontera para el extremo derecho de la barra por lo que si queremos encontrar el nuevo valor estacionario el problema será:

%Trabajo ejercicio 6
clear all , clc , clf

%datos iniciales
a=0; b=4; h=0.1;
x=a:h:b;
N=(b-a)/h; %el vector x tiene n+1 elementos
q=2; %término que acompaña a Uxx
%nodos interiores
xx=x(1:N);
xx=xx';
%cond contorno
ua=0; ub=0;
%cond inicial
u0=exp(-6*(xx-2).^2)+5-5*xx/4;
% matriz K
K=q/h^2*(2*diag(ones(1,N))-diag(ones(1,N-1),-1) -diag(ones(1,N-1),1));
K(1,2)=-4/h^2;
% término F 
F=0*xx;
F(1)=F(1)-4*ua/h^2;
F(end)=F(end)+q*ub/h^2;
% discretización en t
t0=0;tM=50;
k=0.1; % paso temporal.
t=t0:k:tM;
M=length(t)-1; %número de subintervalos temporales
sol1(:,1)=u0;
for i=1:M
  %Euler implicito
  sol1(:,1+i)=(eye(size(K))+k*K)\(sol1(:,i)+k*F);
end
%Añado frontera(solo Dirichlet)
Ub=q*ub*ones(1,length(t));
sol1=[sol1;Ub];

%gráfico
[Mt,Mx]=meshgrid(t,x);
figure(1)
mesh(Mx,Mt,sol1); 
xlabel('Longitud');
ylabel('Tiempo');
zlabel('Temperatura');

Se sacan 2 conclusiones claras: La primera es que el valor estacionario de la temperatura de la barra es 0(se puede demostrar analíticamente también), porque se ve que a partir de un determinado tiempo la gráfica se mantiene constante a temperatura 0.

La segunda que el tiempo que tarda la temperatura en alcanzar ese estado estacionario está alrededor de los 20 segundos. Para comprobar exactamente cuanto tiempo tarda en alcanzar ese estado estacionario con un error del 5 % se procede a dibujar un gráfico 2D para ver exactamente el tiempo pedido:

CONCLUSIÓN: Se alcanza el estado estacionario con un error menos del 5% a partir de los 21 segundos.

Se solicita a continuación, resolver la misma situación anterior pero usando el método de Fourier con distintos términos de la serie(1,3,5,10 y 20). Dibujamos las aproximaciones en una misma gráfica para t=0.5.

%Ejercicio 7
clear all, clc , clf
%datos del problema
a=0;b=4; h=0.1;
x=a:h:b;
Q=[1,3,5,10,20];
u0=exp(-6*(x-2).^2)+5-5*x/4; %valor inicial
t=0.5;
hold on
plot(x,u0,'m')
for i=1:length(Q)
 U=0;
 for j=1:Q(i)
 p=cos((j-(1/2))*pi/4*x);
 aj=trapz(x,u0.*p)/trapz(x,p.*p);
 U=U+aj*exp(-2*(j-1/2)^2*(pi^2)*t/16)*cos((j-(1/2))*pi*x/4);
end
plot(x,U)
xlabel('Longitud');
ylabel('Temperatura');
end
legend('u0','Q=1','Q=3','Q=5','Q=10','Q=20','Location','best')
hold off

El objetivo de esta gráfica es ver la temperatura a lo largo de la barra cuando han pasado 0.5 unidades de tiempo. Hemos hecho la representación usando q= 1,3,5,10 y 20, aunque en la gráfica solo se observan la de q=1 y q=20, puesto que las demás se encuentran superpuestas. Teóricamente la más exacta es la de q=20, pero en ambas podemos observar como la temperatura en el inicio de la barra parte de 3,7 grados y a medida que nos acercamos al final de esta va decreciendo hasta llegar a 0, cumpliendo así nuestra condición de contorno u(4,t)=0. También hemos añadido la temperatura inicial de la barra para poder compararlas, observando así que hasta x=2.5 la temperatura ha decrecido respecto a la inicial y a partir de 2.5 ha subido levemente.

5 Pérdida de calor de la varilla

Por último se supone que se produce una pérdida de calor de la varilla a través del aire. Esta pérdida es constante y de 16 ºC. ¿Cómo se traduce esta pérdida de calor a la ecuación?

La nueva ecuación teniendo en cuenta este nuevo factor quedará: Ut-Uxx+U = 16 ; que manteniendo las condiciones de contorno de la primera situación del trabajo se nos pide encontrar el tiempo para alcanzar el estado estacionario:

%Trabajo ejercicio 8 
clear all , clc , clf

%datos iniciales
a=0; b=4; h=0.1;
x=a:h:b;
N=(b-a)/h; %el vector x tiene n+1 elementos
%nodos interiores
xx=x(2:N);
xx=xx';
%cond contorno
ua=5; ub=0;
%cond inicial
u0=exp(-6*(xx-2).^2)+5-5*xx/4;
% matriz K
K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
I=eye(size(K));
% término F
F=16*xx;
F(1)=F(1)+16+ua/h^2;
F(end)=F(end)+16+ub/h^2;
% discretización en t
t0=0;tM=10;
k=0.0005; % paso temporal.
t=t0:k:tM;
M=length(t)-1; %número de subintervalos temporales
sol(:,1)=u0;
for i=1:M
%Euler explicito
sol(:,1+i)=sol(:,i)+k*(-(K+I)*sol(:,i)+F);
end
%incluimos cond contorno
Ua=ua*ones(1,length(t));
Ub=ub*ones(1,length(t));
sol=[Ua;sol;Ub];

[Mt,Mx]=meshgrid(t,x); 
mesh(Mx,Mt,sol)
xlabel('Longitud');
ylabel('Tiempo');
zlabel('Temperatura');

En este caso de pérdida de calor parece que la solución estacionaria se obtiene antes,alrededor de los 4-5 segundos. Veamos en que momento exacto con error de 10^-3 se alcanza este estado obteniendo la gráfica en 2D:

Así se ve que la temperatura estacionaria es 4.655 que se mantiene constante(estacionaria) a partir de 4.485 segundos. Con esto se obtiene lo pedido; el tiempo que se alcanza el estado estacionario con un error de 10^-3

6 Condición de flujo saliente

El último experimento del trabajo, es resolver la ecuación anterior pero con una condición de frontera diferente:

                                       Ux(0,t) = 2, U(1,t)= -2  ; t(0,10) x(0,1)

La interpretación física de estas últimas condiciones es que existe un flujo de calor saliente por el extremo izquierdo(X=0) de valor constante 2ºC y en el extremo derecho(X=1) se ha colocado un objeto a temperatura constante de -2ºC. La longitud de la varilla ahora es de 1 m.

Resolvemos esta cuestión:

%Trabajo apartado 9
%Ecuacion Ut - Uxx + u = 16

clear all , clc , clf

%datos iniciales
a=0; b=1; h=0.1;
x=a:h:b;
N=(b-a)/h; %el vector x tiene n+1 elementos

%nodos interiores
xx=x(1:N);
xx=xx';
%cond contorno
ua=-2; ub=2;
%cond inicial
u0=exp(-6*(xx-2).^2)+5-5*xx/4;
% matriz K
K=1/h^2*(2*diag(ones(1,N))-diag(ones(1,N-1),-1) -diag(ones(1,N-1),1));
K(1,2)=-2/h^2;
I = eye(size(K));
% término F 
F=16*xx;
F(1)=F(1)-(2*ua/h^2)+16;
F(end)=F(end)+16+ub/h^2;
% discretización en t
t0=0;tM=10;
k=0.0005; % paso temporal.
t=t0:k:tM;
M=length(t)-1; %número de subintervalos temporales
sol(:,1)=u0;
for i=1:M
  %Euler explicito
  sol(:,1+i)=sol(:,i)+k*(-(K+I)*sol(:,i)+F);
end
%Añado frontera(solo Dirichlet)
Ub=ub*ones(1,length(t));
sol=[sol;Ub];

%gráfico
[Mt,Mx]=meshgrid(t,x);
mesh(Mx,Mt,sol); 
xlabel('Longitud');
ylabel('Tiempo');
zlabel('Temperatura');