Ecuación del calor en una placa en forma de anillo

Trabajo realizado por estudiantes
Título	Ecuación del calor en una placa en forma de anillo. Grupo 6B
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2013-14
Autores	Elisabet Sánchez López, Ana Santos Martín, Isabel Sastre Furones, Ángel Antonio Villa Figueroa
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

Tenemos una placa plana en forma de anillo. Vamos a trabajar en coordenadas polares en el plano y sabemos que la placa está comprendida entre los radios ρ=1 y ρ=6 y su temperatura inicial viene definida por la función u(ρ,0):

[math]\left\{\begin{matrix}100(ρ-1) ……… ρ ∈ [1, 2] \\100 ……… ρ ∈ [2, 5]\\90(6-ρ)+10 ……… ρ ∈ [5, 6]\end{matrix}\right.[/math]

Si colocamos en los extremos interior ρ=1 y exterior ρ=6 objetos a una temperatura constante de 0 y 10 grados respectivamente, es decir, las condiciones de frontera serán:

[math]\left\{\begin{matrix}u(1,t)=0 \\u(6,t)=10\end{matrix}\right.[/math]

[math] u_t-\Delta u= u_t- \frac{\partial^2 u }{\partial ρ^2} -\frac{1 }{ρ} \frac{\partial u }{\partial ρ}-\frac{1 }{ρ^2} \frac{\partial^2 u }{\partial θ^2}= u_t- \frac{\partial^2 u }{\partial ρ^2} -\frac{1 }{ρ} \frac{\partial u }{\partial ρ} = 0 [/math]

Así la ecuación queda de la siguiente manera:

[math] u_t-\frac{1 }{ρ} u_ρ - u_(ρρ) = 0 [/math]

Suponiendo que la temperatura u de la placa únicamente depende de la coordenada radial y del tiempo, es decir, u=u(ρ,t) y sabiendo que satisface la ecuación del calor [math] u_t-\Delta u=0[/math] , vamos a plantear el sistema de ecuaciones que satisface u(ρ,t) en coordenadas polares:

2 Resolución por el método de diferencias finitas

2.1 Método del trapecio en el tiempo

Tomamos los pasos de discretización del espacio [math] \Delta ρ=0.1[/math] y del tiempo [math] \Delta t=(\Delta ρ/4)[/math] en t ∈ [0; 10]

%Datos
a=1;b=6;ya=0;yb=10;h=0.1;dt=h/4;
N=(b-a)/h;r=[a:h:b];rr=[a+h:h:b-h];
%Definimos K1
K1=1/h^2*(2*diag(ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1));
%Definimos K2
K2=1/(2*h)*( diag(ones(N-2,1),-1)-diag(ones(N-2,1),1) );
M=diag(1./rr);
K2=M*K2;

T0=0;TF=10;
t=T0:dt:TF;

%Definimos la malla
[R,T]=meshgrid(r,t);
%Definimos F
F=zeros(1,N-1);
F(end)=F(end)+10/h^2+10/(2*5.9*h);
%Condición inicial
for i=1:length(rr)
    if rr(i)<2
        u(i)=100*(rr(i)-1);
    elseif rr(i)<5
        u(i)=100;
    else
        u(i)=90*(6-rr(i))+10;
    end
end
sol(1,:)=[0,u,10];
%Método del trapecio
for m=1:length(t)-1
    Z=u'+(dt/2)*(-(K1+K2)*u'+F');
    W=eye(N-1)+(dt/2)*(K1+K2);
    W=W\Z;
    sol(m+1,:)=[0,W',10];
    u=W';
end
%Resultados gráficos
surf(R,T,sol)

Superficie temperatura en el eje radial ρ∈[1,6] en el intervalo de tiempo

Como se puede observar en la gráfica, la temperatura va disminuyendo hasta alcanzar valores casi constantes a lo largo del tiempo. En el centro (del eje radial) se ve un descenso más brusco, debido a que la temperatura inicial es más alta, y en los extremos se mantiene el mismo valor. Por lo tanto, existe un enfriamiento de la placa, sobre todo en la parte central, hasta unos valores constantes de temperatura: 0 en ρ=1, 10 en ρ=6 y entre 0-10 en ρ∈(1,6).

Temperatura en ρ=3

Partiendo del programa anterior, queremos ahora dibujar el comportamiento de la temperatura en los puntos para los cuales ρ=3 en una gráfica 2D temperatura/tiempo. Implementamos el siguiente código Matlab y obtenemos dicha gráfica.

%Cálculo del punto N que corresponde a ρ=3
N3=(3-a)/h+1;
%Dibujo
plot(sol(:,N3),t)

Temperatura en los puntos de la corona ρ=3

Con esto, se demuestra lo dicho en el anterior apartado, pero se puede ver de forma mucho más clara para ρ=3 que existe un enfriamiento en la corona, llegando a una temperatura entre los 0 y 10ºC.

i

Temperatura en los puntos de la corona ρ=3

Que demuestra lo comentado anteriormente que cuando se supera el intevalo de tiempo se llegan a valores constantes.

2.2 Método de Euler explícito en el tiempo

Puesto que los datos son los mismos que los usados en el método de diferencias finitas con el método del trapecio en el tiempo, únicamente incluimos en este apartado la implementación del nuevo método utilizado.

%Método de Euler explícito
for m=1:length(t)-1
    Z=U'+(dt)*(-(K1+K2)*U'+F');
    sol(m+1,:)=[0,Z',10];
    U=Z';
end       
surf(R,T,sol)

En esta gráfica resultante se puede observar que, como el método es explícito, el sistema es inestable. Esto significa que el paso tomado es grande, por eso tiene los límites muy elevados.

2.3 Método de Euler implícito en el tiempo

Análogamente, en este apartado se especifica la aplicación del método de Euler implícito, obviando los datos (que continúan siendo los mismos).

%Método de Euler implícito:
for m=1:length(t)-1 
    Z=U'+dt*(F');
    W=eye(N-1)+dt*(K1+K2);
    W=W\Z;
    sol(m+1,:)=[0,W',10]
    U=W';
end

Obtenemos una superficie muy parecida a la que salía con el método del trapecio.

2.4 Método de Euler modificado en el tiempo

Al utilizar el método de Euler modificado en el tiempo, observamos que éste es inestable, al igual que el método de Euler explícito.

%Método de Euler modificado:
for m=1:length(t)-1
    M1=-(K1+K2)*U'+F';
    M2=-(K1+K2)*(U'+M1*dt)+F';
    W=U'+(dt/2)*(M1+M2);
    sol(m+1,:)=[0,W',10];
    U=W';
end

3 Estado estacionario

3.1 Con las mismas condiciones

El sistema a resolver para tiempos grandes queda de la siguiente manera:

[math]\left\{\begin{matrix} -\frac{1 }{ρ} u_ρ - u_{ρρ} = 0 \\ u(1)=0 ; u(6)=10\end{matrix}\right.[/math]

El paso que se toma es de 0.1 y 0.01. Como se trata de un problema de contorno, lo resolvereos implementando el siguiente código:

%Datos
a=1; b=6; ya=0; yb=10;
h=[0.1,0.01];
color='rg'
%Para obtener ambas soluciones en la misma gráfica
hold on
for i=1:2
    N=(b-a)/h(i); r=[a:h(i):b]';rr=[a+h(i):h(i):b-h(i)]';
    K1=1/h(i)^2*(2*diag(ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1));
    K2=1/(2*h(i))*( diag(ones(N-2,1),-1)-diag(ones(N-2,1),1) );
    M=diag(1./rr);
    K2=M*K2;
    F=zeros(1,N-1);
    F(end)=F(end)+10/h(i)^2+10/(2*5.9*h(i));
    u=(K1+K2)\F';
    uu=[ya;u;yb];
    plot(r,uu,color(i));
end
hold off

Distribución de la temperatura para tiempos grandes con ambos pasos

Como podemos observar, casi no se aprecia la diferencia entre hacerlo con h=0.1 y h=0.01. Además, demuestra que da igual que sean tiempos grandes porque para nuestra placa se obtendrá un valor dentro del intervalo de tiempo dado y, como se había dicho antes, se llegan a valores constantes de temperatura.

Detalle de la imagen anterior: rojo(0.1) y verde(0.01)

En el zoom podemos comprobar como dicha diferencia es mínima.

3.2 Con las nuevas condiciones

Cuando colocamos en la frontera exterior de la placa una pieza aislante, no se produce pérdida de calor en dicho extremo, por lo que el flujo de temperatura en la dirección radial es nulo.

Una vez alcanzado el estado estacionario, el problema que modeliza la temperatura de la placa es:

[math]\left\{\begin{matrix} -\frac{1}{ρ} u_ρ - u_{ρρ} = 0 \\ u(1)=0 ; u_ρ(6)=0\end{matrix}\right.[/math]

Para obtener el valor estacionario de la temperatura de a placa, resolvemos el problema aplicando el cambio de variable [math] u_ρ = v [/math].

Así, tenemos que resolver una ecuación de variables separables cuya solución será: [math] u_ρ = v = \frac{c_1}{ρ} [/math]

Integrando, obtenemos [math] u(ρ)=c_1·log ρ + c_2 [/math].

Finalmente, obligando al cumplimiento de las condiciones, obtenemos el valor de las constantes: [math]\left\{\begin{matrix} u(1)=0 =\gt c_2=0 \\ u_ρ(6)=0 =\gt c_1=0\end{matrix}\right.[/math]

Por lo tanto, [math] u(ρ)=0 [/math].

Ahora comprobamos cuánto tarda la temperatura en alcanzar el estado estacionario con un error del 5%. Para ello, utilizamos un método cualquiera de los progamados anteriormente y añadimos al bucle del tiempo las siguientes líneas:

for m=1:length(t)-1
    ...
    if (max(abs(sol)))<0.05
         break
    end
end
m
t(m)

Así, obtenemos como resultado t=9.9750 (*), que se muestra por pantalla al ejecutar el programa.

4 Disco completo

Ahora consideraremos que la placa ocupa todo el disco ρ<6. La solución sólo dependerá de ρ y t, y la condición de frontera será [math] u(6,t)=0 [/math].

4.1 Problema de autovalores

Para obtener el problema de autovalores asociado (PA), partimos del problema (P), cuyas condiciones se imponen teniendo en cuenta la condición de frontera dada y que las soluciones tienen que ser difereciables en ρ=0. Debido a este último dato, sabemos que en ese punto tiene que existir tangente horizontal (pendiente nula) y obtenemos, así, la condición [math] u_ρ(0)=0 [/math].

El problema (P) inicial es:

[math]\left\{\begin{matrix} u_t -\frac{1}{ρ} u_ρ - u_{ρρ} = 0 \\ u_ρ(0)=0 ; u(6)=0\end{matrix}\right.[/math]

Aplicando el método de separación de variables, ya que nuestro problema es homogéneo en las condiciones de frontera, buscamos que [math] u(ρ,t) [/math] sea el producto de dos funciones que dependan respectivamente de una de las variables. Así: [math] u(ρ,t)=R(ρ)T(t) [/math].

Derivando, sustituyendo en (P) y dividiendo entre [math] R(ρ)T(t) [/math], llegamos al problema de autovalores asociado (PA):

[math]\left\{\begin{matrix} R''(ρ) +\frac{R'(ρ)}{ρ} = -λR(ρ) \\ R'(0)=0 ; R(6)=0\end{matrix}\right.[/math]

4.2 Relación con las funciones de Bessel

La función de Bessel de primera especie, [math] J_0(r) [/math], es la solución de la ecuación diferencial de Bessel:

[math]\left\{\begin{matrix} r^2 φ''(r) + r φ'(r) + r^2 φ(r) = 0 \\ φ(0)=1 ; φ'(0)=0\end{matrix}\right.[/math]

Para demostrar que las autofunciones del problema anterior son de la forma [math] φ_k = J_0 (ρλ^{1/2}) [/math], siendo λ el autovalor correspondiente, nos basamos en las similitudes entre a ecuación diferencial de Bessel y nuestro problema de autovalores. [math] R = J_0 (ρλ^{1/2}) [/math] tiene que satisfacer (PA). Esto se comprueba sin más que derivar aplicando la regla de la cadena y sustituir en el problema de autovalores, lo que nos proporciona una ecuación diferencial de Bessel de la cual R ya es solución.

Si [math] φ_k = J_0 (ρλ^{1/2}) [/math] son las autofunciones del problema que nos ocupa, siendo λ los autovalores, los ceros de la función de Bessel se relacionan con los autovalores de la siguiente forma: [math] φ_k (6)=0 \\ J_0 (6·λ^{1/2})=0 [/math]

Esto sucede si [math] 6·λ^{1/2} [/math] es un cero de la función de Bessel. Por tanto, si llamamos [math] r_0 [/math] a las raíces de la función de Bessel, obtenemos la relación: [math] r_0 = 6·λ^{1/2} [/math]

4.3 Autovalores y autofunciones

Autovalores

Una vez demostrado que las autofunciones de (P) son de la forma [math] φ_k = J_0 (ρλ^{1/2}) [/math], para calcular los cinco primeros autovalores utilizaremos la relación anterior e implementaremos un pequeño cálculo vectorial en Matlab.

Tomamos como condición [math] φ_k (6)=0 [/math], obtenemos la relación [math] r_0 = 6·λ^{1/2} [/math] y utilizamos los valores de los cinco primeros ceros de la función de Bessel dados. Así, tenemos:

%Cálculo vectorial para la obtención de los 5 primeros autovalores de (P)
r0=[2.4048,5.5201,8.6537,11.7915,14.9309]; %Vector de ceros de la función de Bessel
A=r0.^2/36 %Vector de autovalores

Los resultados que obtenemos son: [math] λ_1 = 0.1606, λ_2=0.8464, λ_3=2.0802, λ_4=3.8622, λ_5=6.1925 [/math]

Autofunciones

Para obtener la gráfica con las primeras cinco autofunciones, utilizaremos la función besselj de Matlab. J=besselj(nu,Z) computa la función de Bessel de orden nu (en nuestro caso, nu=0) para cada elemento del vector Z (en nuestro caso, el vector de [math] ρλ^{1/2} [/math]). Programando en Matlab, obtenemos los siguientes resultados:

%Datos del problema
clear all
a=0; b=6;
r0=[2.4048,5.5201,8.6537,11.7915,14.9309]; %Vector de ceros de la función de Bessel
A=r0.^2/36; %Vector de autovalores
K=5; %Número de autovalores y autofunciones buscadas
%Datos de la discretización
h=0.1;
N=(b-a)/h;
%Vector de radios
p=a:h:b;
%Matriz de ceros de la función de Bessel y matriz de autofunciones
J=zeros(N,K);
z=zeros(N,K);
%Inicialización de las matrices
j0=besselj(0,0);
J(:,1)=j0;
z0=0;
z(:,1)=z0;
for k=1:K
    for i=1:N
        z(i,k)=(p(1,i))'*sqrt(A(k));
        J(i,k)=besselj(0,z(i,k));
    end
end
%Representación gráfica de las autofunciones J
plot (J)

4.4 Resolución por el método de Fourier

Por último, aproximaremos la solución de la ecuación del calor utilizando cinco términos del desarrollo de Fourier cuando el dato inicial es [math] u(ρ,0)=-log \frac{|ρ+0.1|}{6.1} [/math]. Ensayamos con soluciones de la forma:: [math] u(ρ,t)= T_k(t)φ_k(ρ)\\ k=1,2,3,4,5 [/math]

Debe satisfacer la condición inicial del problema:: [math] u(ρ,0)= T_k(0)J_0(ρλ_k^{1/2}) =-log \frac{|ρ+0.1|}{6.1} [/math]

De la serie de Fourier de [math] h(ρ)=-log \frac{|ρ+0.1|}{6.1} [/math], obtenemos la expresión de los coeficientes de Fourier que se utilizará para programar en Matlab. En ella, tendremos en cuenta que, al tratarse de funciones de Bessel, las autofuciones son ortogonales respecto al producto escalar.

El programa de Matlab (*) para obtener la solución de la ecuación del calor por el método de Fourier es:

a=0; b=6;
r0=[2.4048,5.5201,8.6537,11.7915,14.9309]; %Vector de ceros de la función de Bessel
A=r0.^2/36; %Vector de autovalores
K=5; %Número de autovalores y autofunciones buscadas
%Datos de la discretización
h=0.1;
N=(b-a)/h;
%Vector de radios
p=a:h:b;
%Matriz de ceros de la función de Bessel y matriz de autofunciones
J=zeros(N+1,K);
z=zeros(N+1,K);
%Inicialización de las matrices
j0=besselj(0,0);
J(:,1)=j0;
z0=0*p;
z(:,1)=z0';
%Matriz de soluciones
U=0;
%Mallado en el tiempo
t0=0; tf=20;
m=0.1; %Longitud de paso
t=t0:m:tf;
[Mp,Mt]=meshgrid(p,t);
%Bucle
for k=1:K
   z(:,k)=(p)'*sqrt(A(k));
   J(:,k)=besselj(0,z(:,k));
   h=-log((abs(p+0.1))/6.1); %Vector de condiciones iniciales
   ck=trapz(p,p'.*h'.*J(:,k))/trapz(p,p'.*J(:,k).*J(:,k)); %Coeficientes de Fourier
   U=U+ck*exp(-A(k)*t)'*J(:,k)';
end
%Representación gráfica de la solución obtenida
mesh(Mp,Mt,U)

La gráfica obtenida es la siguiente:

(*) Ha sido modificado para la exposición oral.