Transmisión de calor de una varilla metálica (15B)
 Este artículo está en versión beta. El autor de este artículo no lo ha terminado todavía, por favor no lo edites hasta que elimine este mensaje.
|
Contenido
1 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.
1.1 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):
Condiciones de frontera (son de tipo Dirichlet):
La condición inicial:
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:
Imponiendo que v(x,t) sea de la forma
y que cumpla las condiciones de contorno del problema, con lo que se obtiene:
Mediante el cambio, la función w(x,t) nos define el siguiente sistema, que sí tiene las condiciones de contorno homogéneas:
Una vez resuelto w(x,t), la solución u(x,t) buscada será:
2 Métodos numéricos
2.1 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:
%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:
2.2 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: Apartado4_1.m
%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:
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.
Y en este caso:
No se, cumple por lo tanto basta con buscar un Δt que haga que el cociente sea ½, este Δt=0,005. El nuevo código sería: Códigos\Apartado4_1b.m
%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
Método de Euler Implícito:
Apartado4_2.m
%Apartado 4_2 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=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 explicito')
xlabel('Longitud')
ylabel('Tiempo')
zlabel('Temperatura')
3 Ap3
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: Apartado3.m La gráfica obtenida:
4 Ap5
5. Comprobar gráficamente 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. ¿Qué ecuaciones crees que debería verificar este estado estacionario? ¿A qué función se debe parecer la solución u(x; t) para tiempos grandes? Calcular la distancia entre esa solución estacionaria y la solución u(x; t) para t = 0; 1; 2; 10 Utilizando el método de Euler explícito tanteamos con distintos t, el código es el mismo que antes pero cambiamos el valor de t y tenemos las siguientes gráficas: Con t=0
Con t=1
Con t=2
Con t=5
Con t=10
Con t=50
Como vemos a partir de t=5 la solución toma un valor estacionario de 0. Esto indica que la varilla está totalmente fría y aislada.
5 Ap6
6. Ahora colocamos en el extremo derecho una pieza aislante (en lugar de un objeto con temperatura constante 10). El aislante hace que no haya pérdida de calor en ese extremo, es decir, el flujo de temperatura es nulo. ¿Cuál es el valor estacionario de la temperatura de la barra? ¿Cuánto tarda la temperatura en alcanzar el estado estacionario con un error del 5%? En este caso se produce un cambio en las condiciones de frontera para el extremo derecho de la barra por lo que nuestro problema ahora sería:
6 Ap7
7. Resolver el apartado anterior usando el método de Fourier con 1, 3, 5, 10 y 20 términos de la serie. Dibujar las aproximaciones en una misma gráfica para t = 0:5.
7 Ap8
8. Supongamos ahora que se produce una pérdida de calor de la varilla a través del aire que tiene una temperatura constante de 16 grados. En este caso la ecuación es ut - uxx + (u-16) = 0. Si mantenemos las condiciones de contorno del apartado primero, calcular el estado estacionario y el tiempo en que se alcanza con un error de 10-3. En este caso el sistema que modeliza nuestro problema es el siguiente:
8 Ap9
9. Resolver la ecuación anterior con las condiciones de contorno u(0; t) = 10 sin(t), ux(1; t) = 3. Dibujar la solución en t 2 [0; 10]. ¿Qué interpretación física tienen estas condiciones de frontera?
