Una placa uniformemente iluminada

En la página titulada Se calienta una placa expuesta al Sol Consideramos una placa de área A y espesor e, que está a la temperatura ambiente Ts. La superficie de la placa está pintada de negro. En un instante dado, se expone al Sol que ilumina la placa con una intensidad constante de I W/m2. Determinamos la evolución de la temperatura T de la placa a medida que transcurre el tiempo.

Supongamos una placa de espesor L, a la temperatura ambiente Ts, que se expone al sol que ilumina uniformemente su cara x=0, con q W/m2

Para calcular la temperatura T(x,t) en un punto x de la placa 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

Las condiciones de contorno en el extremo izquierdo x=0, y en el extremo derecho x=L son:

T(x,t) x | x=0 = h 1 ( T(0,t) T s ) q K t>0 T(x,t) x | x=L = h 2 ( T(L,t) T s )t>0

El primer término, proporcional a la diferencia de temperaturas T(0,t)-Ts entre el extremo izquierdo x=0 de la placa y la temperatura ambiente, da cuenta de las pérdidas de calor que experimenta la placa a través de su cara anterior. El segundo término, q es el fujo de calor uniforme (energía por unidad de área y de tiempo) que incide sobre dicha superficie. El término proporcional a la diferencia de temperaturas T(L,t)-Ts entre el extremo derecho x=L de la placa y la temperatura ambiente, da cuenta de las pérdidas de calor que experimenta la placa a través de su cara posterior. Las constantes de proporcionalidad h1 y h2 podrían ser diferentes

La condición inicial, es la de una placa a la temperatura uniforme T(x,0)=Ts del ambiente

El estado estacionario

Al cabo de cierto tiempo, teóricamente infinito, se establece un estado estacionario en el que la temperatura de cada punto de la placa no varía con el tiempo.

2 T x 2 =0 T(x)=Ax+B

donde A y B se determinan a partir de las condiciones de contorno en x=0

A= h 1 ( B T s ) q K

y en x=L

A= h 2 ( AL+B T s )

Despejamos A y B del sistema de dos ecuaciones

T( x, )=( q K )( h 2 x+ h 2 L+1 h 1 + h 2 + h 1 h 2 L )+ T s

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 ·contorno u(x,t) x | x=0 = h 1 u(0,t) u(x,t) x | x=L = h 2 u(L,t) ·inicialu(x,0)=( q K )( h 2 x h 2 L1 h 1 + h 2 + h 1 h 2 L )

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)=Acos(ωx)+Bsin(ωx)

Condiciones de contorno

Eliminando A y B en el sistema homogéneo de dos ecuaciones y expresando la ecuación resultante en términos de k=ωL

( h 1 h 2 L 2 k 2 )sin(k)+( h 1 + h 2 )Lkcos(k)=0

Calculamos las raíces kn de esta ecuación transcendente

La solución completa u(x,t) es la superposición de los productos de las funciones Fn(xGn(x) para cada kn

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

Ya que la condición de contorno en el extremo izquierdo x=0, requiere que Bω=h1A o bien, An=Bnkn/(h1L)

Condición inicial

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

u(x,0)= n=1 B n X n X n (x)= k n h 1 L cos( k n x L )+sin( k n x L )

Multiplicando ambos miembros por Xm(x) e integrando entre 0 y L

0 L u(x,0) X m (x)dx = n=1 B n 0 L X m (x) X n (x)dx

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

0 L X m (x) X n (x)dx =0mn

Obtenemos los coeficientes Bm

B m = 0 L ( q K )( h 2 x h 2 L1 h 1 + h 2 + h 1 h 2 L ) X m (x)dx 0 L X m 2 (x)dx

Cambiando m por n e integrando

B n =( q/K h 1 + h 2 + h 1 h 2 L ) ( 1 h 1 + h 2 L 2 k n 2 )sin( k n )+ ( h 1 + h 2 )L h 1 k n cos( k n ) L( h 1 + h 2 + h 1 h 2 L ) h 1 k n ( 1+ k n 2 ) 2 h 1 + 1 4 ( k n h 1 L k n )sin(2 k n ) 1 2 h 1 cos(2 k n )

Solución completa

La solución completa es

u(x,t)= n=1 B n ( k n h 1 L cos( k n x L )+sin( k n x L ) )exp( α k n 2 L 2 t ) T(x,t)=u(x,t)+T(x,)= T s +( q K )( h 2 x+ h 2 L+1 ( h 1 + h 2 + h 1 h 2 L ) ){ n=1 B n ( k n h 1 L cos( k n x L )+sin( k n x L ) )exp( α k n 2 L 2 t ) }

Ejemplo

Supongamos que la anchura de la placa es L=0.5 y los coeficientes h1=10 y h2=10. Vamos a representar gráficamente la función y=(h1h2L2-x2)sin(x)+(h1+h2)Lxcos(x), para estimar aproximadamente donde se encuentran las raíces de la ecuación transcendente.

> L=0.5; % anchura placa
>> h1=10; %pérdidas por radiación
>> h2=10; %pérdidas por radiación
>> f=@(x) (h1*h2*L^2-x^2)*sin(x)+(h1+h2)*L*x*cos(x);
>> fplot(f,[1,50])
>> 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_10 para calcular la distribución de temperaturas a lo largo de la dirección perpendicular a la placa en el instante t

function [x,T]=temperatura_10(Ts,L,h1,h2,q_K,k,a2,t)
   x=linspace(0,L,100);
   if(t==0)
       T=Ts*ones(1,length(x));
       return;
   end
      
   cte=q_K/(h1+h2+h1*h2*L);
   T=Ts+cte*(-h2*x+h2*L+1);
    
    for n=1:length(k)
       num=(-1/h1+h2*L^2/k(n)^2)*sin(k(n))+(h1+h2)*L*cos(k(n))
/(h1*k(n))-L*(h1+h2+h1*h2*L)/(h1*k(n));
       den=(1+k(n)^2)/(2*h1)+(k(n)/h1-L/k(n))*sin(2*k(n))/4-cos(2*k(n))/(2*h1);
       T=T+(cte*num/den)*exp(-k(n)^2*t/(a2*L^2))*(k(n)*cos(k(n)*x/L)
/(h1*L)+sin(k(n)*x/L));
   end
end

Creamos un script en el que establecemos la temperatura ambiente Ts, la anchura L de la placa, el material del cual está hecho la placa, (parámetro α), los coeficientes de pérdida de calor por radiación en ambas caras h1 y h2 y el cociente q/K entre el flujo de calor q que incide uniformemente sobre la cara izquierda y la conductividad térmica del material K. El script calcula un número elevado de raíces kn de la ecuación transcendente y representa la distribución de temperaturas a lo largo de la dirección perpendicular a la placa, en varios instantes. Se sugiere al lector probar con materiales de baja conductividad.

Ts=20; %temperatura ambiente
alfa=11352; %Coeficiente, 1/alfa del aluminio
L=0.5; % longitud placa
h1=10; %pérdidas por radiación
h2=10; %pérdidas por radiación
q_K=10000/209.3; %conductividad del aluminio

%raíces
x=linspace(0,30,10);
f=@(x) (h1*h2*L^2-x.^2).*sin(x)+(h1+h2)*L*x.*cos(x);
k=raices(f,x);

hold on
for t=[5 50 100 200]
    [x,T]=temperatura_10(Ts,L,h1,h2,q_K,k,alfa,t);
    plot(x,T,'displayName',num2str(t));
end
title('Evolución de la temperatura de la placa')
xlabel('x')
ylabel('T')
legend('-DynamicLegend','location','northeast')
grid on
hold off  

Creamos un script para mostrar la evolución de la temperatura de la superficie iluminada x=0, que tiende hacia un valor límite constante

Ts=20; %temperatura ambiente
alfa=11352; %Coeficiente, 1/alfa del aluminio
L=0.5; % longitud placa
h1=10; %pérdidas por radiación
h2=10; %pérdidas por radiación
q_K=10000/209.3; %conductividad del aluminio

%raíces
x=linspace(0,30,10);
f=@(x) (h1*h2*L^2-x.^2).*sin(x)+(h1+h2)*L*x.*cos(x);
k=raices(f,x);
t=5:5:500; %tiempo
T=zeros(1,length(t));
cte=q_K/(h1+h2+h1*h2*L);
for i=1:length(t)    
    T(i)=Ts+cte*(h2*L+1);
    for n=1:length(k)
        num=(-1/h1+h2*L^2/k(n)^2)*sin(k(n))+(h1+h2)*L*cos(k(n))/
(h1*k(n))-L*(h1+h2+h1*h2*L)/(h1*k(n));
        den=(1+k(n)^2)/(2*h1)+(k(n)/h1-L/k(n))*sin(2*k(n))/4-cos(2*k(n))
/(2*h1);
        T(i)=T(i)+(cte*num/den)*exp(-k(n)^2*t(i)/(alfa*L^2))*(k(n)/(h1*L));
    end
end
plot(t,T)

title('Temperatura de la superficie iluminada')
xlabel('t')
ylabel('T')
grid on

Una vez obtenidas las raíces kn de la ecuación transcendente, comprobamos las relaciones de ortogonalidad entre las funciones Xn(x) y Xm(x) con m≠n

>> f1=@(x) cos(k(1)*x/L)+h1*L*sin(k(1)*x/L)/k(1);
>> f2=@(x) cos(k(2)*x/L)+h1*L*sin(k(2)*x/L)/k(2);
>> f3=@(x) f1(x).*f2(x);
>> integral(f3,0,L)
ans =   9.3675e-17

>> f1=@(x) cos(k(7)*x/L)+h1*L*sin(k(7)*x/L)/k(7);
>> f2=@(x) cos(k(5)*x/L)+h1*L*sin(k(5)*x/L)/k(5);
>> f3=@(x) f1(x).*f2(x);
>> integral(f3,0,L)
ans =   1.0408e-16

Nota

La primera raíz de la ecuación transcendente es k=0, esta raíz no está incluida por que produce divisiones entre cero o una indeterminación 0/0 al calcular los coeficientes Bn

Comprobamos con MATLAB que el coeficiente B es cero cuando k→0

>> syms k h1 h2 L;
>> num=(-1/h1+h2*L^2/k^2)*sin(k)+(h1+h2)*L*cos(k)/(h1*k)-
L*(h1+h2+h1*h2*L)/(h1*k);
>> den=(1+k^2)/(2*h1)+(k/h1-L/k)*sin(2*k)/4-cos(2*k)/(2*h1);
>> limit(num/den,k,0)
ans =0

Referencias

E Marín, A Lara-Bernal, A Calderón and O Delgado-Vasallo On the heat transfer through a solid slab heated uniformly and continuously on one of its surfaces Eur. J. Phys. 32 (2011) 783-791