Ver el código de «Transmisión de calor de una varilla metálica (15B)»

[[Categoría:ED13/14]]
==Introducción==
Consideramos una varilla delgada de longitud L que ocupa el intervalo x ϵ (0;L) con L = 3 m., y cuyos extremos están colocados sobre objetos que mantienen una temperatura constante de 0 y 10 grados respectivamente. 


Inicialmente la temperatura de la varilla viene dada por u0(x) = 10x/3 salvo en su tercio central donde la temperatura ha subido hasta los 100 grados.

===Sistema de ecuaciones===

Suponiendo que la temperatura u(x, t) de la varilla satisface la ecuación del calor ut-uxx = 0, se ha planteado un sistema completo de ecuaciones que satisface u(x, t). Para desarrollar el caso, se ha modelizado una varilla de 3 metros de longitud, calor específico 1, densidad 1 y conductividad térmica 1. Designamos por u(x,t) la temperatura en la posición x de la en un tiempo t. Sus extremos mantienen la temperatura constante a 0 y 10 grados respectivamente y en el instante inicial la temperatura de la misma es 10x/3 salvo en el tercio central donde la temperatura inicial es 100ºC.  El sistema de ecuaciones que modeliza el problema será el siguiente:
Ecuación (ecuación diferencial en derivadas parciales):
[[Image:Trabajo15BFoto1.jpg|center|400px|b]]
Condiciones de frontera (son de tipo Dirichlet):
[[Image:Trabajo15BFoto2.jpg|center|400px|b]]
La condición inicial:
[[Image:Trabajo15BFoto3.jpg|center|400px|Método de Diferencias finitas y Trapecio]]
Este desarrollo resulta en un problema con condiciones de contorno no homogéneas, por lo que se ha realizado un cambio de variable de la forma:
[[Image:Trabajo15BFoto4.jpg|center|300px|Método de Diferencias finitas y Trapecio]]
Imponiendo que v(x,t) sea de la forma
[[Image:Trabajo15BFoto5.jpg|center|300px|Método de Diferencias finitas y Trapecio]]
y que cumpla las condiciones de contorno del problema, con lo que se obtiene:
[[Image:Trabajo15BFoto6.jpg|center|100px|Método de Diferencias finitas y Trapecio]]
Mediante el cambio, la función w(x,t) nos define el siguiente sistema, que sí tiene las condiciones de contorno homogéneas:
[[Image:Trabajo15BFoto7.jpg|center|400px|Método de Diferencias finitas y Trapecio]]

Una vez resuelto w(x,t), la solución u(x,t) buscada será:
[[Image:Trabajo15BFoto8.jpg|center|300px|Método de Diferencias finitas y Trapecio]]

==Métodos numéricos==
===Trapecio===
Una vez obtenido el sistema que define el desarrollo del calor, se ha resuelto por el método numérico de diferencias finitas para aproximar una solución. En este caso, se ha aplicado con Δx = 0,1. Asimismo, se ha resuelto mediante el método del trapecio en tiempo tomando Δt = Δx.

El código de Matlab® utilizado para resolver el problema por el método de las diferencias finitas y el trapecio es el siguiente:

{{matlab|codigo=
%Apartado 2 Método Diferencias finitas y del trapecio
%Discretizacion espacial
L=3; %Longitud de la varrilla
x0=0;xN=L;
dx=0.1;
N=(xN-x0)/dx;
xi=dx:dx:L-dx;
x=0:dx:L;
U01=0;
U02=100*ones(1,29)-(10/3).*xi;
d=1;c=1;k=1;nn=k/(d*c);
aa=1;bb=0;
cc=1;dd=0;
KK=-2*diag(ones(1,N-1))+diag(ones(1,N-2),1)+diag(ones(1,N-2),-1);
p=(2*aa*dx+bb)/(aa*dx+bb);
q=(2*cc*dx+dd)/(cc*dx+dd);
KK(1,1)=-p;
KK(N-1,N-1)=-q;
K=(nn/dx^2)*KK;
%Discretizacion temporal
t0=0;
tM=1; %Tiempo de estudio
dt=0.1;
M=(tM-t0)/dt;
t=0:dt:tM;
theta=0.5;
I=eye(N-1);
P=(I-dt*theta*K)\(I+dt*(1-theta)*K);
U=zeros(M+1,N-1);
U(1,11:21)=U02(1,11:21);
for n=1:M
    U(n+1,:)=U(n,:)*P';
end
u0=(bb/aa*dx+bb)*U(:,1);
uN=(dd/(cc*dx+dd))*U(:,N-1);
UU=[u0,U,uN];
[X,T]=meshgrid(x,t);
V=(10/3*ones(M+1,31)).*(X);
figure(100)
S=UU+V;
surf(X,T,S) %Representamos la solución
title('Método del trapecio')
xlabel('Longitud')
ylabel('Tiempo')
zlabel('Temperatura')
}}
El resultado es la siguiente gráfica:
[[Image:Trabajo15BApartado2.jpg|center|thumb|500px|Método de Diferencias finitas y Trapecio]]


Se puede observar en la imagen los cambios bruscos producidos por la escasa precisión del método del trapecio. En los siguientes métodos se observará cómo se suavizan estos cambios

===Euler===
Buscando un análisis comparativo de los métodos numéricos, se ha desarrollado el mismo problema aplicando el método de Euler explícito, así como el método de Euler implícito.
El sistema a resolver es el mismo que antes, si bien los códigos utilizados ahora serán:
Método de Euler explícito:
{{matlab|codigo=
%Apartado 4_1 Método Diferencias finitas y Euler explicito
%Discretizacion espacial
L=3; %Longitud de la varrilla
x0=0;xN=L;
dx=0.1;
N=(xN-x0)/dx;
xi=dx:dx:L-dx;
x=0:dx:L;
U01=0;
U02=100*ones(1,29)-(10/3).*xi;
d=1;c=1;k=1;nn=k/(d*c);
aa=1;bb=0;
cc=1;dd=0;
KK=-2*diag(ones(1,N-1))+diag(ones(1,N-2),1)+diag(ones(1,N-2),-1);
p=(2*aa*dx+bb)/(aa*dx+bb);
q=(2*cc*dx+dd)/(cc*dx+dd);
KK(1,1)=-p;
KK(N-1,N-1)=-q;
K=(nn/dx^2)*KK;
%Discretizacion temporal
t0=0;
tM=1; %Tiempo de estudio
dt=0.1;
M=(tM-t0)/dt;
t=0:dt:tM;
theta=0;
estab=dt*nn/dx^2;
I=eye(N-1);
P=(I-dt*theta*K)\(I+dt*(1-theta)*K);
U=zeros(M+1,N-1);
U(1,11:21)=U02(1,11:21);
for n=1:M
    U(n+1,:)=U(n,:)*P';
end
u0=(bb/aa*dx+bb)*U(:,1);
uN=(dd/(cc*dx+dd))*U(:,N-1);
UU=[u0,U,uN];
[X,T]=meshgrid(x,t);
V=(10/3*ones(M+1,31)).*(X);
figure(100)
S=UU+V;
surf(X,T,S) %Representamos la solución
title('Método de Euler explícito')
xlabel('Longitud')
ylabel('Tiempo')
zlabel('Temperatura')
}}
La gráfica que obtenemos hoy es la siguiente:
[[Image:Trabajo15BApartado4_1a.jpg|center|thumb|500px|Método de Euler explícito]]

Como se puede observar, la solución no tiene nada que ver con la anterior, 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:Δx/(Δt)^2 (sube la foto)

Y en este caso no se cumple por lo tanto basta con buscar un Δt que haga que el cociente sea ½, tomando un Δt=0,005. 
El nuevo código sería:
{{matlab|codigo=
%Apartado 4_1 Método Diferencias finitas y Euler explicito
%Discretizacion espacial
L=3; %Longitud de la varrilla
x0=0;xN=L;
dx=0.1;
N=(xN-x0)/dx;
xi=dx:dx:L-dx;
x=0:dx:L;
U01=0;
U02=100*ones(1,29)-(10/3).*xi;
d=1;c=1;k=1;nn=k/(d*c);
aa=1;bb=0;
cc=1;dd=0;
KK=-2*diag(ones(1,N-1))+diag(ones(1,N-2),1)+diag(ones(1,N-2),-1);
p=(2*aa*dx+bb)/(aa*dx+bb);
q=(2*cc*dx+dd)/(cc*dx+dd);
KK(1,1)=-p;
KK(N-1,N-1)=-q;
K=(nn/dx^2)*KK;
%Discretizacion temporal
t0=0;
tM=1; %Tiempo de estudio
dt=0.005;
M=(tM-t0)/dt;
t=0:dt:tM;
theta=0;
estab=dt*nn/dx^2;
I=eye(N-1);
P=(I-dt*theta*K)\(I+dt*(1-theta)*K);
U=zeros(M+1,N-1);
U(1,11:21)=U02(1,11:21);
for n=1:M
    U(n+1,:)=U(n,:)*P';
end
u0=(bb/aa*dx+bb)*U(:,1);
uN=(dd/(cc*dx+dd))*U(:,N-1);
UU=[u0,U,uN];
[X,T]=meshgrid(x,t);
V=(10/3*ones(M+1,31)).*(X);
figure(100)
S=UU+V;
surf(X,T,S) %Representamos la solución
title('Método de Euler explícito')
xlabel('Longitud')
ylabel('Tiempo')
zlabel('Temperatura')
}}
Y la gráfica

[[Image:Trabajo15BApartado4_1b.jpg|center|thumb|500px|Método de Euler explícito, ahora corregido]]

Método de Euler Implícito:
{{matlab|codigo=
%Apartado 4_2 Método Diferencias finitas y Euler implicito
%Discretizacion espacial
L=3; %Longitud de la varrilla
x0=0;xN=L;
dx=0.1;
N=(xN-x0)/dx;
xi=dx:dx:L-dx;
x=0:dx:L;
U01=0;
U02=100*ones(1,29)-(10/3).*xi;
d=1;c=1;k=1;nn=k/(d*c);
aa=1;bb=0;
cc=1;dd=0;
KK=-2*diag(ones(1,N-1))+diag(ones(1,N-2),1)+diag(ones(1,N-2),-1);
p=(2*aa*dx+bb)/(aa*dx+bb);
q=(2*cc*dx+dd)/(cc*dx+dd);
KK(1,1)=-p;
KK(N-1,N-1)=-q;
K=(nn/dx^2)*KK;
%Discretizacion temporal
t0=0;
tM=1; %Tiempo de estudio
dt=0.1;
M=(tM-t0)/dt;
t=0:dt:tM;
theta=1;
I=eye(N-1);
P=(I-dt*theta*K)\(I+dt*(1-theta)*K);
U=zeros(M+1,N-1);
U(1,11:21)=U02(1,11:21);
for n=1:M
    U(n+1,:)=U(n,:)*P';
end
u0=(bb/aa*dx+bb)*U(:,1);
uN=(dd/(cc*dx+dd))*U(:,N-1);
UU=[u0,U,uN];
[X,T]=meshgrid(x,t);
V=(10/3*ones(M+1,31)).*(X);
figure(100)
S=UU+V;
surf(X,T,S) %Representamos la solución
title('Método de Euler implicito')
xlabel('Longitud')
ylabel('Tiempo')
zlabel('Temperatura')
}}
[[Image:Trabajo15BApartado4_2imp.jpg|center|thumb|500px|Método de Euler implícito]]

==Desarrollo práctico del modelo==
===Variación del modelo en un punto determinado===
Habiendo desarrollado una solución numérica, se puede usar esta para definir el comportamiento de la temperatura en el punto medio de la barra en una gráfica 2D temperatura/tiempo. Usando el método de diferencias finitas y trapecio, ya desarrollado, bastará con suprimir la línea de la función surf (representar superficie) y añadir una con la función plot para obtener la gráfica deseada.
El código para este apartado es el siguiente:
{{matlab|codigo=
%Apartado 3: Gráfica temperatura/tiempo
%Discretizacion espacial
L=3; %Longitud de la varilla
x0=0;xN=L;
dx=0.1;
N=(xN-x0)/dx;
xi=dx:dx:L-dx;
x=0:dx:L;
U01=0;
U02=100*ones(1,29)-(10/3).*xi;
d=1;c=1;k=1;nn=k/(d*c);
aa=1;bb=0;
cc=1;dd=0;
KK=-2*diag(ones(1,N-1))+diag(ones(1,N-2),1)+diag(ones(1,N-2),-1);
p=(2*aa*dx+bb)/(aa*dx+bb);
q=(2*cc*dx+dd)/(cc*dx+dd);
KK(1,1)=-p;
KK(N-1,N-1)=-q;
K=(nn/dx^2)*KK;
%Discretizacion temporal
t0=0;
tM=1; %Tiempo de estudio
dt=0.1;
M=(tM-t0)/dt;
t=0:dt:tM;
theta=0.5;
I=eye(N-1);
P=(I-dt*theta*K)\(I+dt*(1-theta)*K);
U=zeros(M+1,N-1);
U(1,11:21)=U02(1,11:21);
for n=1:M
    U(n+1,:)=U(n,:)*P';
end
u0=(bb/aa*dx+bb)*U(:,1);
uN=(dd/(cc*dx+dd))*U(:,N-1);
UU=[u0,U,uN];
[X,T]=meshgrid(x,t);
V=(10/3*ones(M+1,31)).*(X);
figure(100)
S=UU+V;
plot(t,S(:,16)) %Representamos la gráfica temperatura/tiempo
title('Gráfico temperatura/tiempo')
xlabel('Tiempo')
ylabel('Temperatura')
}}
La gráfica obtenida:

[[Image:Trabajo15BApartado3.jpg|center|thumb|500px|Evolución de la temperatura]]

===Desarrollo del modelo 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. Para demostrarlo gráficamente, basta con utilizar el método de Euler explícito, tanteando con distintos tiempos. El código es el mismo que antes pero al variar el valor de t, se obtienen las siguientes gráficas:

Con t=0
[[Image:Trabajo15BApartado5_0.jpg|center|thumb|400px|Desarrollo en t 0]]
Con t=1
[[Image:Trabajo15BApartado4_1b.jpg|center|thumb|400px|Desarrollo en t 1]]
Con t=2
[[Image:Trabajo15BApartado5_2.jpg|center|thumb|400px|Desarrollo en t 2]]
Con t=5
[[Image:Trabajo15BApartado5_5.jpg|center|thumb|400px|Desarrollo en t 5]]
Con t=10
[[Image:Trabajo15BApartado5_10.jpg|center|thumb|400px|Desarrollo en t 10]]
Con t=50
[[Image:Trabajo15BApartado5_50.jpg|center|thumb|400px|Desarrollo en t 50]]

Como se puede observar, a partir de t=5 la solución toma un valor estacionario de 0. La interpretación física de este problema indica que la varilla está totalmente fría y aislada.

===Modelización con elementos aislantes===
Continuando con la experimentación, al colocar en el extremo derecho una pieza aislante (en lugar de un objeto con temperatura constante 10) 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á:
{{matlab|codigo=
%Apartado 6_1 Método Diferencias finitas y Euler explicito
%Discretizacion espacial
L=3; %Longitud de la varrilla
x0=0;xN=L;
dx=0.1;
N=(xN-x0)/dx;
xi=dx:dx:L-dx;
x=0:dx:L;
U01=(10/3).*xi;
U02=100*ones(1,29);
d=1;c=1;k=1;nn=k/(d*c);
aa=1;bb=0;
cc=0;dd=1;
KK=-2*diag(ones(1,N-1))+diag(ones(1,N-2),1)+diag(ones(1,N-2),-1);
p=(2*aa*dx+bb)/(aa*dx+bb);
q=(2*cc*dx+dd)/(cc*dx+dd);
KK(1,1)=-p;
KK(N-1,N-1)=-q;
K=(nn/dx^2)*KK;
%Discretizacion temporal
t0=0;
tM=30; %Tiempo de estudio
dt=0.005;
M=(tM-t0)/dt;
t=0:dt:tM;
theta=0;
estab=dt*nn/dx^2;
I=eye(N-1);
P=(I-dt*theta*K)\(I+dt*(1-theta)*K);
U1=zeros(M+1,N-1);
U2=zeros(M+1,N-1);
U3=zeros(M+1,N-1);
U1(1,2:10)=U01(1,2:10);
U2(1,11:20)=U02(1,11:20);
U3(1,21:29)=U01(1,21:29);
U=U1+U2+U3;
for n=1:M
    U(n+1,:)=U(n,:)*P';
end
u0=(bb/aa*dx+bb)*U(:,1);
uN=(dd/(cc*dx+dd))*U(:,N-1);
UU=[u0,U,uN];
[X,T]=meshgrid(x,t);
figure(100)
surf(X,T,UU) %Representamos la solución
title('Método de Euler explícito')
xlabel('Longitud')
ylabel('Tiempo')
zlabel('Temperatura')
pause
%Representamos la x=3
plot(t,UU(:,31))
xlabel('Longitud')
ylabel('Temperatura')
}}
[[Image:Trabajo15BApartado6.jpg|center|thumb|500px|Flujo con aislante en un extremo]]


[[Image:Trabajo15BApartado6b.jpg|center|thumb|500px|Variación a lo largo de la barrra]]


Vemos que el valor estacionario de la temperatura es cero grados. Para ver con mayor precisión en que instante se alcanza la estacionalidad, dibujamos la gráfica Temperatura-Tiempo del extremo x=3. En dicha gráfica observamos que a partir de x=22.22 la temperatura es cero con un error despreciable (se dispone en dicho punto de una temperatura de 0.1º), desde luego mucho menor que el 5% deseado.