Diferencia entre revisiones de «Calor Placa Anillo (18B)»

Revisión del 18:34 17 may 2014

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

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\)

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

%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)

4.2 Euler implícito

%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)

5 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 analogo 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.

6 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\).

6.1 Problema de autovalores

6.2 Función de Bessel

6.3 Aproximación de la solución