Una varilla cuyo extremo radia calor

Supongamos una varilla metálica de longitud L, que se ha calentado a una temperatura T0, uno de sus extremos x=0, se mantiene a temperatura ambiente Ts y el otro extremo x=L radia calor al ambiente.

Para calcular la temperatura T(x,t) en un punto x de la varilla y en un instante t, resolvemos la ecuación de la conducción del calor

1 α T t = 2 T x 2 0<x<Lt>0

La condición de contorno en el extremo izquierdo x=0, es T(0,t)=Ts

La condición de contorno en el extremo derecho x=L, es

T(x,t) x | x=L =h( T(L,t) T s )

La condición inicial, es la de una varilla a temperatura uniforme T(x,0)=T0

El estado estacionario

Al cabo de cierto tiempo, teóricamente infinito, que en la práctica depende del tipo de material que empleamos, se establece un estado estacionario en el que la temperatura de cada punto de la varilla no varía con el tiempo. La temperatura de la varilla será constante e igual a la temperatura ambiente T(x,∞)=Ts

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

Definimos la función u(x,t)=T(x,t)-T(x,∞), 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 x=0 y en x=L se escriben

1 α u t = 2 u x 2 0<x<Lt>0 ·contornou(0,t)=0 u(x,t) x | x=L +hu(L,t)=0 ·inicialu(x,0)= T 0 T s

Buscamos una solución de la forma u(x, t)=F(xG(t), variables separadas

1 α 1 G(t) dG(t) dt = 1 F(x) d 2 F(x) d x 2 = ω 2

Integramos la primera 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(x) d x 2 + ω 2 F(x)=0

Es una ecuación diferencial similar a la de un Movimiento Armónico Simple, cuya solución es F(x)=Asin(ωx)+Bcos(ωx)

Condiciones de contorno

La solución completa u(x,t) es la superposición de los productos de las funciones Fn(xGn(x) para cada kn que calcularemos resolviendo la ecuación transcendente anterior en k=ωL

u(x,t)= n=1 A n sin( k n x L )exp( α k n 2 L 2 t )

Condición inicial

Solamente, queda por determinar los coeficientes An, identificando la solución para t=0 con la condición inicial u(x,0)

u(x,0) = n=1 A n sin( k n x L ) u(x,0) = T 0 T s

Teniendo en cuenta la relación de ortogonalidad, véase más abajo

0 L sin( k n x L ) sin( k m x L )dx=0 n m

Obtenemos los coeficientes An

A n = 0 L ( T 0 T s )sin( k n x/L)dx 0 L sin 2 ( k n x/L)dx = 8( T 0 T s ) sin 2 ( k n /2 ) 2 k n sin( 2 k n ) = 4( T 0 T s ) k n sin 2 ( k n /2 ) k n 2 +hL sin 2 ( k n )

Para obtener esta última expresión, hemos tenido en cuenta que kn son las raíces de la ecuación transcendente kcos(k)+hLsin(k)=0

Solución completa

La solución completa es

u(x,t)= n=1 A n sin( k n x L )exp( α k n 2 L 2 t ) T(x,t)=u(x,t)+T(x,)= T s +4( T 0 T s ) n=1 k n sin 2 ( k n /2 ) k n 2 +hL sin 2 ( k n ) sin( k n x L )exp( α k n 2 L 2 t )

Ejemplo

Supongamos que hL=1, vamos a representar gráficamente la función y=xcos(x)+hLsin(x), para estimar aproximadamente donde se encuentran las raíces de la ecuación transcendente. Cambiamos el valor del parámetro hL para ver el efecto sobre las raíces

>> hL=1;
>> f=@(x) x*cos(x)+hL*sin(x);
>> fplot(f,[1,100])
>> grid on
>> xlabel('x')
>> ylabel('y')

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_4 para calcular la distribución de temperaturas a lo largo de la varilla en el instante t

function [x,T]=temperatura_4(T0,Ts,L,k,a2,t)
   x=linspace(0,L,100);
   if(t==0)
       T=T0*ones(1,length(x));
       return;
   end
      
   T=Ts*ones(1,length(x));
   for n=1:length(k)
        an=8*(T0-Ts)*sin(k(n)/2)^2/(2*k(n)-sin(2*k(n)));
        T=T+an*exp(-k(n)^2*t/(a2*L^2))*sin(k(n)*x/L);
   end
end

Creamos un script en el que establecemos la temperatura inicial T0, la temperatura ambiente Ts, la longitud de la varilla L, el material del cual está hecho la varilla, (parámetro α), el coeficiente de pérdida de calor por radiación hL. El script calcula un número elevado de raíces kn de la ecuación transcendente y representa la distribución de temperaturas de la varilla en varios instantes. Se sugiere al lector cambiar el valor del parámetro hL y observar el efecto sobre la distribución de temperaturas en la varilla

T0=100; %temperatura inicial
Ts=0; %temperatura ambiente (final)
hL=1; %pérdida de calor por radiación
L=1; % longitud de la varilla
alfa=11352; %Coeficiente, 1/alfa del aluminio

f=@(x) x.*cos(x)+hL*sin(x);
x=linspace(0,400,400);
k=raices(f,x);

%instantes
t=[10 100 500 2000, 5000];
hold on
axis([0 0.5 -5 100]);
for i=1:length(t)
    [x,T]=temperatura_4(T0,Ts,L,k,alfa,t(i));
    plot(x,T,'displayName',num2str(t(i)));
end
title('Evolución de la temperatura de una varilla')
xlabel('Longitud')
ylabel('Temperatura')
legend('-DynamicLegend','location','northeast')
grid on
hold off  

Una vez obtenidas las raíces kn de la ecuación transcendente, comprobamos las relaciones de ortogonalidad entre las funciones sin(knx/L) y sin(kmx/L) con m≠n

>> f=@(x) sin(k(1)*x/L).*sin(k(2)*x/L);
>> integral(f,0,L)
ans =  -6.9389e-17
>> f=@(x) sin(k(5)*x/L).*sin(k(9)*x/L);
>> integral(f,0,L)
ans =     0