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

Revisión actual del 23:59 12 ago 2019

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

En las gráficas se puede observar cómo para tiempo 0, los valores de temperatura para cada radio son los indicados en las condiciones iniciales. A partir de ahí y a medida que va pasando el tiempo el flujo de calor se expande a las zonas de menor temperatura, esto es, de los radios interiores (2 a 5) a los exteriores, de manera que para un tiempo igual a 10 segundos la placa pasa a tener menor temperatura en los radios interiores, llegando a tener todos los radios una misma temperatura. Sin embargo se aprecia que en el extremo del anillo con radio igual a 6 tenemos unas condiciones físicas que hacen que en ese extremo la temperatura valga 10. Esto hace que para tiempo igual a 10 segundo no toda la placa este homogeneizada en temperatura,puesto que en ese extremo la temperatura está obligada a ser 10.

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:10; 
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

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, y por tanto, se obtienen las mismas conclusiones.

4.3 Euler modificado

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

%euler modificado (euler+trapecio)
clear all 
h=0.1; 
a=1; b=2; c=5; d=6;
N=(d-a)/h;
x=a:h:d;
K=(1/(h^2))*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1) -diag(ones(1,N-2),1));
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,10];
U=U'; 
for k=2:M 
   k1=(K)*U+F';
   k2=((j/2)+(K))*U+F';
   U=U-(j/2)*(k1+k2);
   sol(k,:)=[0,U',10]; 
end
[X,T]=meshgrid(x,t);
mesh(X,T,sol)

Euler modificado

Vemos que gráficamente no se aprecian variaciones, siendo prácticamente iguales todas las gráficas. La única diferencia es el método utilizado para conseguirlo, es decir, el código MATLAB utilizado; mientras que los métodos de Euler explícito y modificado precisan de un paso temporal pequeño (del orden de [math]h^2/2[/math] , siendo \(h\) el paso de discretización espacial) para ser estables y conseguir una buena aproximación, no es así con el método de Euler implícito, que no necesita un paso de tan pequeño orden y la aproximación conseguida es bastante acertada. Hay que añadir que el inconveniente del uso de un paso pequeño, no es más que para intentar el menor colapso del programa matemático utilizado (MATLAB en este caso), y conseguir el resultado en el menor tiempo posible.

Después de este razonamiento, vemos que el método más óptimo es de Euler implícito.

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 por Fourier se expone a continuación. Utilizaremos k=10 para realizar la aproximación, 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 temperatura a lo largo de la placa para diferentes tiempos t=0,1,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( \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})}}$ desarrollada en serie de Fourier, en forma de $e^{-((k\cdot{\pi}^2)\cdot{t}} $. Esta funcion tiende a 0 para tiempos grandes. 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

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

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

Vemos que para t=10 segundo la solución ya coincide prácticamente con al 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. El sistema en este caso es análogo al anterior. La diferencia que se impone es que el flujo de calor en la dirección 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.

La variación de la temperatura a lo largo del disco no es uniforme, si no que el extremo\(\rho=6\) tiene una temperatura mayor que \(\rho=1\). Esto es debido a la condición inicial de colocar un elemento aislante en el extremo. En las gráficas expuestas a continuación se puede ver la temperatura en la placa para diferentes tiempos. Para un tiempo grande de 50 segundos podemos ver que la temperatura ya es prácticamente 0 en toda la placa.

Temperaturas en la placa para t=1 seg

Temperaturas en la placa para t=2 seg

Temperaturas en la placa para t=50 seg

7 Transformación del problema en disco

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}

\]

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

7.1 Problema de autovalores

7.2 Función de Bessel

7.3 Aproximación de la solución