Conducción del calor en un cilindro muy largo (I)

El cilindro se calienta a una temperatura uniforme T0 y se deja enfriar sumergiéndolo en un baño térmico a temperatura Ts.

En el estado estacionario, después de un tiempo t→∞, la temperatura final del cilindro 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 = 2 T r 2 + 1 r T r 0r<R,t>0

r es la distancia radial desde el eje del cilindro.

La condición de contorno es T(R,t)=Ts

La distribución inicial de temperaturas es T(r,0)=T0

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 + 1 r u r 0<r<Rt>0 ·contornou(R,t)=0 ·inicialu(x,0)= T 0 T s

Buscamos la solución mediante el procedimiento de separación de variables, de la forma u(r,t)=F(rG(t)

F α dG dt =G( d 2 F d 2 r + 1 r dF dr ) 1 F ( d 2 F d 2 r + 1 r dF dr )= 1 α 1 G dG dt = k 2 R 2

Resolvemos la ecuación diferencial de la derecha

dG dt + k 2 α R 2 G=0 G(t)=G(0)exp( k 2 α R 2 t )

Resolvemos la ecuación diferencial de la izquierda

1 F ( d 2 F d 2 r + 1 r dF dr )= k 2 R 2

F(r)=J0(kr/R), J0 es la función de Bessel de primera especie y de orden cero.

Condiciones de contorno

La condición de contorno en r=R es T(R,t)=Ts, o bien, u(R,t)=F(RG(t)=0. Las raíces kn, son los valores de k que cumplen la ecuación transcendente J0(k)=0

La solución completa u(r,t) es la superposición de los productos Fn(rGn(t)

u(r,t)= n=1 A n J 0 ( k n r R ) exp( k n 2 α R 2 t )

Condición inicial

Los coeficientes An se determinan a partir de la distribución inicial de temperaturas u(r,0)

u(r,0)= T 0 T s = n=1 A n J 0 ( k n r R )

Despejamos An utilizando las relaciones de ortogonalidad de las funciones de Bessel, más abajo

( T 0 T s ) 0 R r J 0 ( k m r R )dr= n=1 A n 0 R r J 0 ( k m r R ) J 0 ( k n r R )dr ( T 0 T s ) 0 R r J 0 ( k m r R )dr= A m 0 R r J 0 2 ( k m r R )dr

Haciendo el cambio de variable x=r/R y teniendo en cuenta los resultados de las integrales

( T 0 T s ) 0 1 x J 0 ( k m x )dx= A m 0 1 x J 0 2 ( k m x )dx ( T 0 T s ) J 1 ( k m ) k m = 1 2 A m ( J 0 2 ( k m )+ J 1 2 ( k m ) ) A m = 2( T 0 T s ) J 1 ( k m ) ( J 0 2 ( k m )+ J 1 2 ( k m ) ) k m

>> syms x k;
>> int('x*besselj(0,k*x)^2',x,0,1)
ans= besselj(0, k)^2/2 + besselj(1, k)^2/2
>> int('x*besselj(0,k*x)',x,0,1)
ans =besselj(1, k)/k

Cambiando m por n y teniendo en cuenta que kn son las raíces de la ecuación transcendente J0(k)=0. Los coeficientes An valen

A n = 2( T 0 T s ) k n J 1 ( k n )

Solución completa

u(r,t)=2( T 0 T s ) n=1 1 k n J 1 ( k n ) J 0 ( k n r R ) exp( k n 2 α R 2 t ) T(r,t)=u(r,t)+T(r,)= T s +2( T 0 T s ) n=1 1 k n J 1 ( k n ) J 0 ( k n r R ) exp( k n 2 α R 2 t )

Ejemplo

Representamos gráficamente la función J0(x), para estimar aproximadamente donde se encuentran las raíces de la ecuación transcendente J0(x)=0.

>> f=@(x) besselj(0,x);
>> fplot(f,[0,100])
>> xlabel('x')
>> ylabel('J_0(x)')
>> title('Función J_0(x)')
>> grid on

La función raices calcula las raíces múltiples de la función f(x) buscando los intervalos en los que dicha función cambia de signo y utilizando la función MATLAB fzero para encontrarlas

function r = raices(f, x)
    y=f(x);
    indices=find(y(1:end-1).*y(2:end)<0);
    r=zeros(1,length(indices));
    for k=1:length(indices)
        r(k)=fzero(f, [x(indices(k)), x(indices(k)+1)]);
    end
end

Definimos la función temperatura_5 para calcular la distribución de temperaturas a lo largo de de la dirección radial del cilindro en el instante t

function [r,T]=temperatura_5(T0,Ts,R,k,a2,t)
   r=linspace(0,R,100);
   if(t==0)
       T=T0*ones(1,length(r));
       return;
   end
      
   T=Ts*ones(1,length(r));
   for n=1:length(k)
        an=2*(T0-Ts)/(k(n)*besselj(1,k(n)));
        T=T+an*exp(-k(n)^2*t/(a2*R^2))*besselj(0,(k(n)*r/R));
   end
end

Creamos un script en el que establecemos la temperatura inicial T0, la temperatura ambiente Ts, el radio R de un cilindro muy largo, el material del cual está hecho el cilindro, (parámetro α). El script calcula un número elevado de raíces kn de la ecuación transcendente y representa la distribución radial de temperaturas del cilindro en varios instantes.

T0=100; %temperatura inicial
Ts=0; %temperatura ambiente
alfa=11352; %Coeficiente, 1/alfa del aluminio
R=1; % radio del cilindro

%raíces
x=linspace(0,100,100);
f=@(x) besselj(0,x);
k=raices(f,x);

hold on
axis([0 1 -5 100]);
for t=[10 500 1000 2000]
    [x,T]=temperatura_5(T0,Ts,R,k,alfa,t);
    plot(x,T,'displayName',num2str(t));
end
title('Evolución de la temperatura de un cilindro')
xlabel('r')
ylabel('Temperatura')
legend('-DynamicLegend','location','southwest')
grid on
hold off  

Una vez obtenidas las raíces kn de la ecuación transcendente, J0(k)=0, comprobamos las relaciones de ortogonalidad que hemos utilizado para calcular los coeficientes An

0 1 x J 0 ( k m x) J 0 ( k n x)dx=0mn

>> f=@(x) x.*(besselj(0,k(1)*x).*besselj(0,k(2)*x));
>> integral(f,0,1)  
ans =   8.9311e-18