Calor Placa Anillo (18B)

Trabajo realizado por estudiantes
Título	Ecuación del calor en una placa en forma de anillo (Grupo 18)
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2013-14
Autores	• Arantxa Abascal Colomar • Patricia Fernández Aibar • Paula Lacanal Cuadrado • David Ortiz Liriano • Álvaro Pintor Sousa • Alberto Rodríguez Fernández
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

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.

1 Introducción

Anillo circular entre los radios \(\rho=1\) y \(\rho=6\)

El problema nos pide considerar una placa plana en forma de anillo o corona circular, comprendida entre los radios \(\rho=1\) y \(\rho=6\), dividida en diferentes franjas circulares (sectores), cuya temperatura varía según la dirección radial.

Disponemos de la función que determina la temperatura inicial de la placa, que abarca distintos valores del radio \(\rho\):

\[ \begin{array}{c} \\u(\rho,0)=\begin{cases} 100(\rho - 1) & \text{ si } \rho \epsilon (1,2) \\ 100 & \text{ si } \rho \epsilon (2,5) \\ 90(6-\rho)+10 & \text{ si } \rho \epsilon (5,6) \end{cases} &&&(1) \end{array} \]

Los extremos de la placa están en contacto con objetos que, por diversas razones físicas, hacen que el límite interior de la placa se encuentre a \(0ºC\), mientras que el exterior está, de un modo similar, a \(10ºC\). Estas temperaturas se mantienen constantes con el paso del tiempo, es decir, durante \(t>0\).

Analizando la situación, podemos ver que se trata de un problema que se puede representar mediante la ecuación del calor y su distribución a lo largo del tiempo para los diferentes intervalos representados anteriormente \((1)\).

OBSERVACIÓN: Tal como indica el enunciado, se hace preciso trabajar con la ayuda de una parametrización en coordenadas polares por la facilidad práctica que supone. Sin embargo, la temperatura no varía en función de \(\theta\) (siendo \(\theta\) el ángulo girado con respecto al eje de abcisas y centro en el origen de coordenadas), razón por la cual asemejaremos este problema al de una varilla de longitud comprendida entre \(x\epsilon [1,6]\), cuya función de temperatura depende del tiempo y de la distancia \(\rho\) al extremo interior \(\rho=1\).

Considerando \(u(\rho,t)\) como la función que determinará el calor de la placa, las condiciones iniciales y de frontera, serán las siguientes:

• Para \(\rho=1\), tenemos una temperatura constante de \(0ºC\).
• Para \(\rho=6\), tenemos una temperatura constante de \(10ºC\).

2 Planteamiento del sistema de ecuaciones

Suponemos que la temperatura u de la placa en forma de anillo depende solo de la cordenada radial \(\rho\) y del tiempo t es decir \[

u = u(\rho,t) \]

y que satisface la ecuacion del calor \[

u_t - \Delta u =0 \]

El sistema de ecuaciones que satisface \( u = u(\rho,t)\) es el siguiente: \[

\begin{array}{c} u_t - u_{\rho\rho} =0 &  \rho  \epsilon (1,6)& y & t>0 \\ u(1,t)=0 ,  u(6,t)=10 & t>0 \\u(\rho,0)=\begin{cases} 100(\rho - 1) && \text{ si }  \rho   \epsilon (1,2) \\ 100 && \text{ si } \rho \epsilon (2,5) \\ 90(6-\rho) && \text{ si }  \rho \epsilon (5,6) \end{cases} \end{array}

\]

Resolvemos esta ecuación diferencial mediante el método del trapecio, y los resultados gráficos son los siguientes:

Podemos observar que en el entorno de \(\rho=2\), la gráfica realiza ciertas oscilaciones. Podría ser debido a que en ese punto, la solución de nuestra ecuación diferencial no tiene límite definido, haciendo que nuestra herramienta de cálculo matemático no sepa fijar un valor concreto para cada punto perteneciente a este entorno. Podríamos poner el ejemplo de la función [math]sin (1/x) [/math] , cuyo límite no se define completamente, y la gráfica obtenida es la siguiente:

El código de Matlab utilizado para la obtención de resultados, es el siguiente:

%
clear all 
a=1; b=6; 
h=0.1; 
x=a:h:b; %El vector x tiene N+1 elementos 
N=(b-a)/h;
% matriz K 
K=(1/h^2)*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1)); 
% término F y valor inicial 
rho=x(2:N); %quitamos primer y último elemento de x 
F=(0.*rho); %término independiente 
%valor inicial 
U=(0.*rho);

rho=1;
for i=1:N-1
    if i>1 & i<(2-1)/h
        rho=rho+h;
        U(1,i)=98*(rho-1);
elseif i>(2-1)/h & i<(5-2)/h
    rho=rho+h;
    U(1,i)=102-(2*rho);
elseif i>(5-2)/h & i<(6-5)/h
    rho=rho+h;
    U(1,i)=92*(6-rho);
    end
end

F=0; %AÑADIMOS CONDICIONES CONTORNO
 
 
% discretización en t 
j=h/4; % paso en t. 
t=0:j:3; 
M=length(t); %número de puntos de t 
 
sol(1,:)=[0,U,10]; 
U=U'; 
for k=1:M-1 
 Z=U+j/2*(-K*U); 
 U=(eye(N-1)+K*j/2)\Z; 
 sol(k+1,:)=[0,U',10]; 
end 
 
 
[Mx,Mt]=meshgrid(x,t); 
mesh(Mx,Mt,sol) ;

3 Representación de la temperatura en el tiempo para p=3

Vamos a dibujar el comportamiento de la temperatura para un radio constante \(\rho\)=3 en una gráfica 2D (temperatura/tiempo) para ello utilizamos el siguiente código MATLAB:

clear all 
h=0.1; 
a=2; b=5; 
x=a:h:b 
N=(b-a)/h; 
% Matriz K 
K=1/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1)); 
% término F 
xx=x(2:N); %quitamos primer y último elemento de x 
F=0*xx; 
U=100*ones(size(xx)); %valor inicial para un /rho=3
% discretización en t 
j= 0.5*h^2; % paso en t. Menor h2 para estabilidad con Euler. 
t=0:j:10; 
M=length(t); %número de puntos de t 
sol(1,:)=[0,U,10]; 
U=U'; 
for k=1:M-1 
 U= U + j*(-K*U +F'); 
 sol(k+1,:)=[0,U',10]; 
 end 
plot(t,sol(:,11)) %la columna 11 es la correspondiente a p=3

Varriacion de la temperatura en los 10 primeros segundos para \(\rho=3\)

Se observa que la temperatura inicial (t=0), es de 100 grados, condición impuesta en las condiciones iniciales. Con el paso del tiempo, la tendencia natural de la temperatura a crear un flujo de calor que viaja desde las zonas de temperatura más altas, a las zonas de temperatura más bajas, influye en la estabilización y homogeneización de esta a lo largo de la dirección radial del disco. Es decir, que a medida que pasa el tiempo, el flujo de calor va de las zonas con más temperatura a las de menor, por ello que en el instante inicial su temperatura sea máxima y luego baje, porque está transmitiendo calor. Como puede observarse una vez alcanzados los 10 primeros segundos la placa en ese radio \(\rho=3\) se estabilizada en 0ºC.

4 Resolución del sistema por diferentes métodos de discretizacion

OBSERVACIÓN: El método de Euler necesita un paso de tiempo del orden de [math]h^2/2[/math] para ser estable, siendo h el paso espacial.

4.1 Euler explicito

Resolvemos el sistema de ecuaciones diferenciales esta vez utilizando el método de Euler explícito, con las mismas condiciones que en el método de las diferencias finitas y del trapecio.

El código MATLAB utilizado es el siguiente:

%euler explicito 
clear all 
h=0.1; 
a=1; d=6;
N=(d-a)/h;
x=a:h:d;
%Matriz K
K=(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
K=1/h^2*K;
xx=x(2:N);
F=0*xx; 
%Valores iniciales
xx1=x(2:10);
xx2=x(11:41);
xx3=x(42:50);
U1 = 100*(xx1-1);
U2 =100*ones(size(xx2)); %valor inicial 
U3 =((90*(6-xx3))+10);
%Discretizacion de t
j= 0.5*h^2; % paso en t. Menor h2 para conseguir estabilidad con Euler. 
t=0:j:10; 
M=length(t); %número de puntos de t 
U=[U1,U2,U3];
sol(1,:)=[0,U,10];
U=U'; 
for k=2:M
 U= U + j*(-K*U +F'); 
 sol(k,:)=[0,U',10]; 
end
[Mx,Mt]=meshgrid(x,t);
mesh(Mx,Mt,sol)

Como podemos observar, la temperatura es máxima entre los radios 2 y 5 en el instante inicial. A partir de ese momento a medida que avanzan los segundos la temperatura desciende debido a la tendencia natural de la placa a estabilizarse, de manera que el calor se transmite desde las zonas de mayor temperatura a las de menor, hasta alcanzar una temperatura homogénea.

4.2 Euler implícito

Realizamos el mismo ejercicio con el método de Euler implícito utilizando el siguiente código MATLAB:

%euler implicito (incondicionalmente estable)
clear all 
h=0.1; 
a=1; d=6;
N=(d-a)/h;
x=a:h:d;
%Matriz K
K=(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
K=1/h^2*K;
xx=x(2:N);
F=0*xx; 
xx1=x(2:10);
xx2=x(11:41);
xx3=x(42:50);
%Valores iniciales
U1 = 100*(xx1-1);
U2 =100*ones(size(xx2)); 
U3 =((90*(6-xx3))+10);
%Discretizacion de t
j= 0.5*h^2;  
t=0:j:10; 
M=length(t); %número de puntos de t 
U=[U1,U2,U3];
sol(1,:)=[0,U,10];
U=U'; 
for k=2:M
 U=(eye(N-1)+j*K)\(U+F');
 sol(k,:)=[0,U',10]; 
end
[Mx,Mt]=meshgrid(x,t);
mesh(Mx,Mt,sol)

En este gráfico se observa lo mismo que lo dicho en el gráfico anterior. La única diferencia es el método utilizado para conseguirlo, es decir en el código MATLAB utilizado. Las conclusiones que se pueden deducir son las mismas que con el método de Euler explícito.

5 Apartado 5

5.1 Solución Analítica

Resolvemos la ecuación del calor en la placa y representamos la solución imponiendo diferentes intervalos de tiempo. El programa de representación es el siguiente.

clear all
N=50;
a=1; b=2; c=5; d=6;
h1=(b-a)/N;
h2=(c-b)/N;
x1=a:h1:b;
x2=b:h2:c;
x3=c:h1:d;
M=51;
k=10/M;
for m=1:N
uaprox =zeros(1,M);
for n=1:M
n1=2*(n-1)+1;
fac1=(-100/n1*pi)*exp(-(n1^2)*(pi^2)*(m-1)*k);
fac2=(200/n1*pi)*exp((-(n1^2)*(pi^2)*(m-1)*k)/9);
fac3=(-110/n1*pi)*exp(-(n1^2)*(pi^2)*(m-1)*k);
q1=sin(n1*pi*x1);
q2=sin(n1*pi*x2);
q3=sin(n1*pi*x3);
uaprox = uaprox + (fac1*q1)+(fac2*q2)+(fac3*q3);
sol(n,:)=[uaprox];
end 
Z=length(uaprox);
h=(d-a)/(Z-1);
x=[a:h:d];
hold on
plot(x,uaprox)
end

Podemos ver en las siguientes gráficas, representando la solución para t=10, 100 y 1000 segundos, como la temperatura tiende a estabilizarse en 0 siguiendo una función senoidal.

6 Apartado 6

Ahora colocamos en la frontera exterior de la placa (\rho = 6) 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 ujo de temperatura en la dirección radial es nulo. >Cual es el valor estacionario de la temperatura de la placa? >Cuanto tarda la temperatura en alcanzar el estado estacionario con un error del 5%?

El sistema en este caso es análogo al anterior. La diferencia que se impone es que el flujo de calor en la direccion radial sea nulo en la placa.(Condición de Neumann)

El sistema de ecuaciones que satisface \( u = u(\rho,t)\) es el siguiente: \[

\begin{array}{c} u_t - u_{\rho\rho} =0 &  \rho  \epsilon (1,6)& y & t>0 \\ u(1,t)=0 ,  ux(6,t)=0 & t>0 \\u(\rho,0)=\begin{cases} 100(\rho - 1) && \text{ si }  \rho   \epsilon (1,2) \\ 100 && \text{ si } \rho \epsilon (2,5) \\ 90(6-\rho) && \text{ si }  \rho \epsilon (5,6) \end{cases} \end{array}

\]

clear all 
h=0.1; 
a=1; b=2; c=5; d=6;
N=(d-a)/h;
x=a:h:d;
K=(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
K(N-1,N-1)=2/3; 
K(N-1,N-2)=-2/3; 
K=1/h^2*K;
xx=x(2:N)
F=0*xx; 
xx1=x(2:10);
xx2=x(11:41);
xx3=x(42:50);
U1 = 100*(xx1-1)
U2 =100*ones(size(xx2)) %valor inicial 
U3 =((90*(6-xx3))+10)
j= 0.5*h^2; % paso en t. Menor h2 para estabilidad con Euler. 
t=0:j:10; 
M=length(t); %número de puntos de t 
U=[U1,U2,U3]
sol(1,:)=[0,U,((-U(4)+U(5))/3)]
U=U'; 
for k=2:M
 U= U + j*(-K*U +F'); 
 sol(k,:)=[0,U',((-U(4)+U(5))/3)]; 
end
[Mx,Mt]=meshgrid(x,t);
mesh(Mx,Mt,sol)

Gráfica que describe el comportamiento de la temperatura a lo largo del tiempo según el radio, de acuerdo a las condiciones del enunciado.

7 Transformación del problema en disco

Disco circular de radio \(\rho<6\)

Considerando que la placa ocupa todo el disco \(\rho<6\), aproximaremos las soluciones usando el método de Fourier. De nuevo, la solución dependerá sólo de \(\rho\) y \(t\), y tomaremos ahora como condición frontera \(u(6;t) = 0\).