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

El cilindro se calienta a una temperatura uniforme T0 y se deja enfriar en un ambiente cuya temperatura constante es Ts.

En el estado estacionario, después de un tiempo t→∞, la temperatura final del cilindro será T(r,∞)=Ts, la temperatura ambiente.

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) r | r=R =h( T(R,t) T s )

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<Lt>0 ·contorno u(r,t) r | r=R +hu(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

Dado que F(r)=J0(kr/R), la condición de contorno da lugar a una ecuación transcendente, cuyas raíces kn son los valores de k que cumplen

d J 0 (kr/R) dr | r=R +h J 0 ( k )=0 k R J 1 (k)+h J 0 ( k )=0

>> syms x k;
>> diff(besselj(0,k*x))
ans =-k*besselj(1, k*x)

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 k·J1(k)=hR·J0(k). Los coeficientes An valen

A n = 2( T 0 T s ) ( k n 2 + h 2 R 2 ) J 0 ( k n )

Solución completa

u(r,t)=2( T 0 T s ) n=1 1 ( k n 2 + h 2 R 2 ) J 0 ( 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 2 + h 2 R 2 ) J 0 ( k n ) J 0 ( k n r R ) exp( k n 2 α R 2 t )

Ejemplo

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

>> hR=10;
>> f=@(x) x*besselj(1,x)-hR*besselj(0,x);
>> fplot(f,[0,100])
>> xlabel('x')
>> ylabel('f(x)')
>> title('Función f(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_7 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_7(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)*besselj(1,k(n))/(k(n)*
(besselj(1,k(n))^2+besselj(0,k(n))^2));
        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 α) y el parámetro hR, cuanto mayor sea su valor, mayor son las pérdidas por radiación de la superficie del cilindro. 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
hR=10; %pérdidas por radiación

%raíces
x=linspace(0,100,100);
f=@(x) x.*besselj(1,x)-hR*besselj(0,x);
k=raices(f,x);

hold on
axis([0 1 -5 100]);
for t=[50 500 1000 2000]
    [x,T]=temperatura_7(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, k·J1(k)-hR·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 =  -2.1034e-17