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

Revisión del 21:55 19 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 \in (1,2) \\ 100 & \text{ si } \rho \in (2,5) \\ 90(6-\rho)+10 & \text{ si } \rho \in (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\in [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 sólo de la coordenada 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 \in (1,6)& t>0 \end{array}\\

\\ \begin{array}{c} \\\\u(1,t)=0 & u(6,t)=10 & t>0 \end{array}\\\\ \\ \begin{array}{c}\\ \\u(\rho,0)=\begin{cases} 100(\rho - 1) & \text{ si } \rho \in (1,2) \\ 100 & \text{ si } \rho \in (2,5) \\ 90(6-\rho)+10 & \text{ si } \rho \in (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:

left

right

left

right

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:

center

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);
% U=[1,N-2];
rho=1;
for i=1:N-1
    if i>1 & i<(2-1)/h
        rho=rho+h;
        U(1,i)=100*(rho-1);
elseif i>(2-1)/h & i<(5-1)/h
    rho=rho+h;
    U(1,i)=100;
elseif i>(5-1)/h & i<(6-1)/h
    rho=rho+h;
    U(1,i)=90*(6-rho)+10;
    end
end

 
% 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 \(\rho=3\)

Vamos a dibujar el comportamiento de la temperatura para un radio constante \(\rho=3\) en una gráfica 2-D (eje de abcisas, temperatura y eje de ordenadas, 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 \rho=3

Variación de la temperatura en los \(10\) primeros segundos para \(\rho=3\)

Se puede observar que la temperatura inicial (\(t=0\)), es de \(100ºC\), 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 alta, a las zonas de temperatura más baja, 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, es por eso que en el instante inicial su temperatura es máxima y luego disminuye, 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]\frac{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)

Euler explícito

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)

Euler Implícito

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.

4.3 Euler modificado

Por último, probamos con el método de Euler modificado. Para ello adjuntamos el siguiente código MATLAB:

Código

Imagen

Discusión

5 Variación de la temperatura para tiempos grandes

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

\begin{array}{c} u_t - u_{\rho\rho} =0 &  \rho \in (1,6)& t>0 \end{array}\\

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

\]

La solución numérica de este sistema consiste en dar como solución aproximada una función en serie de Fourier. \[

\begin{array}{c} U( \rho,t)= U_0( \rho,t) + W( \rho,t) \end{array}\\

\]

\begin{equation}

U( \rho,t)=2x-2+\sum_{i=1}^k{\frac{-230-2.5(-1)^k}{\pi\cdot{k}}}\cdot{e^{-((k\cdot{\pi}^2)\cdot{t}}}\cdot{\sin({k\cdot{\pi}\cdot{x})}} \end{equation}

El programa utilizado para representar la solución se expone a continuación. Utilizaremos k=10 para realizar la aproximacion, ya que conseguimos un error pequeño.

clear all 
a=1; b=6; c=0; d=1;
h=0.1;  N=(b-a)/h;
x=a:h:b;
M=50;
k=(d-c)/M;
t=c:k:d;
[X,T]=meshgrid(x,t);
xx1=x(2:10);
xx2=x(11:41);
xx3=x(42:50);
U1 = 100*(xx1-1);
U2 =100*ones(size(xx2)); 
U3 =((90*(6-xx3))+10);
U0=2*X-2;
Q=10;
U=U0;
for j=1:Q
 pj=sin(pi*x*k);
 pj1=sin(pi*xx1*k);
 pj2=sin(pi*xx2*k);
 pj3=sin(pi*xx3*k);
 fourier(j)=(trapz(xx1,U1.*pj1)+trapz(xx2,U2.*pj2)+trapz(xx3,U3.*pj3))/trapz(x,pj.^2);
 U=U+fourier(j)*(exp(-((pi^2)*j*T))).*(sin(pi*X*j));
 end
figure
mesh(X,T,U)

Representando la superficie para diferentes tiempos t=0,2 y 10 segundos observamos como la temperatura tiende a estabilizarse a un valor estacionario para cada punto, estableciéndose un incremento de temperatura lineal en el anillo, en la dirección del radio. La disminución de la temperatura en el tiempo se plasma dentro de la función \( u = u(\rho,t)\) = XXXXX desarrollada en serie de Fourier, en forma de \(e^-t=1/e^t\). Es precisamente este operador de la función el que impone un valor estacionario de la temperatura en la placa. Este valor estacionario seguirá la función \(U_0( \rho,t)= 2x - 2\), que es la función que nos permite homogeneizar las condiciones de contorno de nuestro problema inicial. De esta forma, para tiempo infinito la temperatura en \(\rho=1\) será 0ºC y en \(\rho=6\) será 10ºC

Solución para un t=0 segundos y diferencia con la solución estacionaria

Solución para un t=1 segundos y diferencia con la solución estacionaria

Vemos que para t=1 segundo la solución ya coincide prácticamente con al estacionaria. Esto significa que la temperatura de la placa se estabiliza muy rápidamente. En la siguiente gráfica se observa como la variación de la temperatura se estabiliza para tiempos pequeños.

Solución para un t=1 segundos y diferencia con la solución estacionaria

6 Variacion de las condiciones de frontera

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 flujo de temperatura en la dirección radial es nulo, no entra ni sale calor en ese extremo de la placa. El valor estacionario de la temperatura de la placa será 0ºC. La temperatura tarda en alcanzar el estado estacionario XXXXXXX 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 \in (1,6)& t>0 \end{array}\\

\\ \begin{array}{c} \\\\u(1,t)=0 & u_\rho(6,t)=0 & t>0 \end{array}\\\\ \\ \begin{array}{c}\\ \\u(\rho,0)=\begin{cases} 100(\rho - 1) & \text{ si } \rho \in (1,2) \\ 100 & \text{ si } \rho \in (2,5) \\ 90(6-\rho)+10 & \text{ si } \rho \in (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\). \[

\begin{array}{c} u_t - u_{\rho\rho} =0 &  \rho \in (0,6)& t>0 \end{array}\\

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

\]

7.1 Problema de autovalores

7.2 Función de Bessel

Mediante la

7.3 Aproximación de la solución