Conducción del calor en una esfera homogénea (I)

Consideremos una esfera de radio R en la cual la distribución inicial de temperaturas y las condiciones de contorno tienen simetría esférica. Las superficies isotérmicas son superficies esféricas concéntricas y la temperatura es una función únicamente de la distancia radial r y del tiempo t.

La esfera se calienta hasta una temperatura uniforme T0 (distribución inicial de temperaturas T(r, 0)=T0. En el instante t=0, se sumerge en un recipiente grande de agua a temperatura Ts que es continuamente agitada. La condición de contorno es por tanto, T(R, t)=Ts.

En el estado estacionario, después de un tiempo t→∞, la temperatura final de la esfera será T(r,∞)=Ts, la temperatura del baño térmico.

La ecuación de la conducción del calor apropiada para resolver este problema es

1 α T t = 1 r 2 r ( r 2 T r ) 1 α T t = 2 r T r + 2 T r 2

Solución de la ecuación de la conducción del calor

Definimos la función u(r,t)=T(r,t)-T(r,∞), en términos de esta nueva función, la ecuación de la conducción del calor, las condición inicial en el instante t=0, y las condiciones de contorno en r=R se escriben

1 α u t = 2 u r 2 + 2 r u r 0<r<Rt>0 ·contornou(R,t)=0 ·inicialu(x,0)= T 0 T s

Si hacemos la sustitución v(r,t)=r·u(r,t) nos queda la educación diferencial en derivadas parciales

v t =α 2 v 2 r

Buscamos soluciones de la forma v(r, t)=F(rG(t) (variables separadas)

v(r,t) t =α 2 v(r,t) 2 r 1 α 1 G(t) dG(t) dt = 1 F(r) d 2 F(r) d 2 r = ω 2

Integramos la primara ecuación diferencial

dG(t) dt +α ω 2 G(t)=0 G(t)=G(0)·exp(α ω 2 t)

Integramos la segunda ecuación diferencial

d 2 F(r) d 2 r + ω 2 F(r)=0

Es una ecuación diferencial similar a la de un Movimiento Armónico Simple. La solución es F(r)=A·sin(ωr)+B·cos(ωr)

En primer lugar, F(0)=0 ya que u(r,t) tiene que ser finito cuando r=0, recuérdese que u(r,t)=v(r,t)/r y v(r, t)=F(rG(t). Esto implica que B=0

Condiciones de contorno

u(R,t)=0, por lo que F(R)=A·sin(ωR)=0, es decir, ωn=nπ/R

La solución v(r,t) es la superposición de los productos Fn(rGn(t). La solución de la ecuación de la conducción del calor r·u(r,t)=v(r,t)

r·u(r,t)= n=1 A n sin( ω n r)·exp(α ω n 2 t)

Condición inicial

La condición inicial es u(r,0)=T0-Ts, que nos permite determinar los coeficientes An

r·u(r,0)= n A n sin( nπ R r )

Multiplicamos ambos miembros por sin(mπr/R) e integramos entre 0 y R, para ello, hacemos el cambio de variable z=πr/R, dz=πdr/R

0 R r·u(r,0)sin( mπ r R )dr = n=1 A n 0 R sin( nπ r R ) sin( mπ r R )dr

El resultado de la integral del segundo miembro cuando m≠n y cuando m=n es

0 R sin( mπ r R )sin( nπ r R )dr = R π 0 π sin( mz )sin( nz )dz ( 1 n 2 m 2 ) 0 π sin( mz )sin( nz )dz = 1 m sin( nz )cos( mz )+ n m 2 cos( nz )sin( nz ) | 0 π =0 0 π sin 2 ( nz )dz = 1 2 0 π ( 1cos( 2nz ) )dz = 1 2 ( z 1 2n sin( 2nz ) ) | 0 π = π 2

Despejamos An

A n = 2 R 0 R r·u(r,0)sin( nπ r R ) dr

Teniendo en cuenta que u(r,0)=T0-Ts y efectuando el cambio de variable zr/R

A n = 2R π 2 ( T 0 T s ) 0 π z·sin(nz)dz = 2R π ( T 0 T s ) ( 1 ) n+1 n

Solución completa

La temperatura T(r,t) en cualquier punto a una distancia r del centro de la esfera, en un instante t, se compone de la suma de la temperatura en el estado estacionario Ts y de una serie rápidamente convergente que describe el estado transitorio.

r·u(r,t)= 2R π ( T 0 T s ) n=1 (1) n+1 n sin( nπ r R )exp( α n 2 π 2 R 2 t ) T(r,t)=u(r,t)+T(r,)= T(r,t)= T s + 2R πr ( T 0 T s ) n=1 (1) n+1 n sin( nπ r R )exp( α n 2 π 2 R 2 t )

Ejemplo

Creamos la función temperatura_1, que admite como parámetros: las temperatura del baño Ts, la temperatura inicial, T0, el radio de la esfera R, el parámetro 1/α y el tiempo t en el que queremos calcular la distribución de temperaturas en la esfera. Devuelve un vector x de puntos a lo largo de la dirección radial y sus correspondientes temperturas en un vector T.

function [x,T]=temperatura_1(Ts,T0,R,a2,t)
   x=linspace(eps,R,100);
    if(t==0)
       T=T0*ones(1,length(x));
       return;
   end  
   T=Ts*ones(1,length(x));
   n=0;
   v=1;
   cte=pi*pi/(a2*R*R);
   k=2*R*(T0-Ts)./(pi*x);
    while(v>0.01)
        n=n+1;
        v=exp(-cte*n*n*t);
        T=T+v*(-1)^(n+1)*k.*sin(pi*n*x/R)/n;
   end
end

Una serie infinita tenemos que interrumpirla en algún término de acuerdo con algún criterio. Un posible criterio sería que la exponencial sea menor que un valor prefijado, por ejemplo 0.01. Nos surge otro problema en el origen x=0, se produce una indeterminación 0/0. Calculamos la temperatura en el centro de la esfera mediante la expresión simplificada, o bien, evitar el punto x=0, tomando un punto muy próximo.

Creamos un script que establece

Define el vector de los instantes t, en segundos, en el que queremos representar la distribución de temperaturas a lo largo de la dirección radial.

Llama a la función temperatura_1 y representa mediante el comando plot cada una de las graficas T(x) en la misma ventana.

Ts=100; %Temperatura, baño
T0=0; %Temperatura inicial 
R=0.5; %Radio esfera
alfa=11352; %Coeficiente, 1/alfa

hold on
axis([0 R -5 100]);
for t=[2 100 200 300 400]
    [x,T]=temperatura_1(Ts,T0,R,alfa,t);
    plot(x,T,'displayName',num2str(t));
end
title('Evolución de la temperatura de una esfera')
xlabel('r')
ylabel('T')
legend('-DynamicLegend','location','northwest')
grid on
hold off 

Temperatura en el centro de la esfera

Temperatura en el centro de la esfera cuando r→0

lim r0 sin( nr π R ) nr = π R

en función del tiempo t

T(0,t)= T s +2( T 0 T s ) n=1 (1) n+1 exp( α n 2 π 2 R 2 t )

Ts=100; %Temperatura, baño
T0=0; %Temperatura inicial 
R=0.5; %Radio esfera
alfa=11352; %Aluminio, coeficiente, 1/alfa

t=4:4:400;
T=Ts*ones(1,length(t));
cte=pi*pi/(alfa*R*R);
for i=1:length(t)
    v=1;
    n=0;
    while(v>0.01)
        n=n+1;
        v=exp(-cte*n*n*t(i));
        T(i)=T(i)+v*2*(T0-Ts)*(-1)^(n+1);
    end
end

plot(t,T);
title('Temperatura en el centro de una esfera')
xlabel('t')
ylabel('T')
grid on  

Referencias

Unsworth J. Duarte F. J. Heat diffusion in a solid sphere and Fourier theory: An elementary practical example. Am. J. Phys. 47 (11) November 1979, pp. 981-983