Equilibrio térmico (II)

Estudiaremos una situación en la que los extremos de las varillas se mantienen aislados.

T 1 (x,t) x | x=0 =0, T 2 (x,t) x | x=L =0

Sus temperaturas iniciales son distintas T1(x,0)=T01 y T2(x, 0)=T02.

En un instante t, la distribución de temperaturas en la parte izquierda es T1(x, t) y en la parte derecha, T2(x, t). Soluciones de las ecuaciones diferencial en derivadas parciales

T 1 (x,t) t = K 1 ρ 1 c 1 2 T 1 (x,t) x 2 ,0x<d T 2 (x,t) t = K 2 ρ 2 c 2 2 T 2 (x,t) x 2 ,dxL

Con las condiciones iniciales y de controno especificadas. En la posición d de contacto de los dos tramos se ha de cumplir

T 1 (d,t)= T 2 (d,t) K 1 T 1 (x,t) x | x=d = K 2 T 2 (x,t) x | x=d

La función T(x,t) es continua y el flujo de calor a través de la sección de la varilla en la posición de contacto d, a uno y otro lado, es el mismo

Una situación similar la encontramos en la propagación de ondas en una cuerda no homogénea

Estado estacionario

Consideremos dos cuerpos, el primero fabricado con un material de calor específico c1 y densidad ρ1 y el segundo, con un material de calor específico c2 y densidad ρ2

Inicialmente, cada cuerpo está a una temperatura T01 y T02, en recintos adiabáticos separados. En el instante t=0, se elimina la pared de separación y se ponen en contacto. Después de un tiempo, teóricamente infinito, los cuerpos adquieren una temperatura común de equilibrio Teq, ya que el calor cedido por el cuerpo caliente Q es absorbido por el cuerpo más frío (se supone que no hay pérdidas de calor).

Q= m 1 c 1 ( T eq T 01 )Q= m 2 c 2 ( T eq T 02 ) T eq = m 1 c 1 T 01 + m 2 c 2 T 02 m 1 c 1 + m 2 c 2 = ( ρ 1 d ) c 1 T 01 +( ρ 2 (Ld) ) c 2 T 02 ( ρ 1 d ) c 1 +( ρ 2 (Ld) ) c 2

En el estado estacionario, T1(x,∞)=Teq, T2(x,∞)=Teq

Supondremos que los dos cuerpos (varillas) tienen la misma sección S. La masa de cada varilla es proporcional a su longitud

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

Definimos las funciones: u1(x,t)=T1(x,t)-T1(x,∞) y u2(x,t)=T2(x,t)-T2(x,∞). En términos de estas nuevas funciones, las ecuaciones de la conducción del calor se esciben,

u 1 (x,t) t = K 1 ρ 1 c 1 2 u 1 (x,t) x 2 ,0x<d u 2 (x,t) t = K 2 ρ 2 c 2 2 u 2 (x,t) x 2 ,dxL

Variables separadas

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

u(x,t) t =α 2 u(x,t) x 2 F(x) dG(t) dt =αG(t) d 2 F(x) d x 2 1 G(t) dG(t) dt = α 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 d 2 F(x) d x 2 + μ 2 ω 2 F(x)=0, μ 2 = 1 α F(x)=Asin( μωt )+Bcos( μωt )

Las soluciones de la segunda, para cada una de las tramos de la varilla, son

F 1 (x)= A 1 sin( μ 1 ωx)+ B 1 cos( μ 1 ωx),0x<d F 2 (x)= A 2 sin( μ 2 ωx)+ B 2 cos( μ 2 ωx),dxL

Tenemos un sistema homogéneo de tres ecuaciones con tres incógnitas B1 ,A2 y B2. Despejamos A2 y B2 en función de B1

B 2 = B 1 cos( μ 1 ωd)·cos( μ 2 ωL) cos( μ 2 ω(Ld)) A 2 = B 1 cos( μ 1 ωd)·sin( μ 2 ωL) cos( μ 2 ω(Ld))

Las funciones F1(x) y F2(x) se expresan

{ F 1 (x)= B 1 cos( μ 1 ωx) F 2 (x)= B 1 cos( μ 1 ωd)·sin( μ 2 ωL) cos( μ 2 ω(Ld)) sin( μ 2 ωx)+ B 1 cos( μ 1 ωd)·cos( μ 2 ωL) cos( μ 2 ω(Ld)) cos( μ 2 ωx) { F 1 (x)= B 1 cos( μ 1 ωx),0x<d F 2 (x)= B 1 cos( μ 1 ωd) cos( μ 2 ω(Ld)) cos( μ 2 ω(Lx)),dxL

Además, obtenemos una ecuación transcendente (el determinante de los coeficientes del sistema homogéneo ha de ser nulo)

K 1 μ 1 sin( μ 1 ωd)cos( μ 2 ω(Ld) )+ K 2 μ 2 sin( μ 2 ω(Ld) )cos( μ 1 ωd)=0

La primera raíz es ω0=0

Utilizamos la función raices que llama al procedimiento fzero de MATLAB para calcular las raíces múltiples de la ecuación transcendente.

Ejemplo

Consideremos dos varillas que difieren en su conductividad térmica

k1=209.3; %aluminio
c1=880;
r1=2700;
k2=45; %acero
c2=460;
r2=7800;
a1=sqrt(c1*r1/k1);
a2=sqrt(c2*r2/k2);
L=1; %longitud varilla
d=L/3; %posición contacto

%raíces wn
w=linspace(0,0.2,50);
f=@(w) sin(a1*w*d).*cos(a2*w*(L-d))+(k2*a2/(k1*a1))*sin(a2*w*(L-d)).*cos(a1*w*d);
raiz=raices(f,w); 
fplot(f,[0,0.2])
xlabel('\omega')
ylabel('f(\omega)')
grid on
title('Raíces de la ecuación transcendente')

raiz =
    0.0129    0.0270    0.0420    0.0570    0.0714    0.0847    0.0974    0.1111    
0.1259    0.1409    0.1557    0.1693    0.1819    0.1953

La representación gráfica nos ayuda a establecer el paso de exploración de la función y el intervalo o el número de raíces que deseamos obtener. Es conveniente, comprobar que las raíces obtenidas numéricamente se corresponden con los ceros de la función f(ω)

Solución completa

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

u 1 (x,t)= n=0 B n cos( μ 1 ω n x)exp( ω n 2 t) ,0x<d u 2 (x,t)= n=0 B n cos( μ 1 ω n d) cos( μ 2 ω n (Ld)) cos( μ 2 ω n (Lx) )exp( ω n 2 t) ,dxL

Donde ω0=0

De forma abreviada expresamos la superposición

u 1 (x,t)= n=0 B n X 1n (x)exp( ω n 2 t) ,0x<d u 2 (x,t)= n=0 B n X 2n (x)exp( ω n 2 t) ,dxL

Relación de ortogonalidad

Se puede demostar que

c 1 ρ 1 0 d X 1n (x) X 1m (x)dx+ c 2 ρ 2 d L X 2n (x) X 2m (x)dx=0

Evaluamos esta relación, utilizando procedimientos numéricos

k1=209.3; %aluminio
c1=880;
r1=2700;
k2=45; %acero
c2=460;
r2=7800;
a1=sqrt(c1*r1/k1);
a2=sqrt(c2*r2/k2);
L=1; %longitud varilla
d=L/3; %posición contacto

%raíces wn
w=linspace(0,0.2,50);
f=@(w) sin(a1*w*d).*cos(a2*w*(L-d))+(k2*a2/(k1*a1))*sin(a2*w*(L-d)).*cos(a1*w*d);
raiz=raices(f,w); 
w1=raiz(3); 
w2=raiz(4);
X1_1=@(x) cos(a1*w1*x);
X1_2=@(x) cos(a1*w2*x);
X2_1=@(x) cos(a1*w1*d)*cos(a2*w1*(L-x))/cos(a2*w1*(L-d));
X2_2=@(x) cos(a1*w2*d)*cos(a2*w2*(L-x))/cos(a2*w2*(L-d));
Z1=@(x) X1_1(x).*X1_2(x);
Z2=@(x) X2_1(x).*X2_2(x);
res=c1*r1*integral(Z1,0,d)+c2*r2*integral(Z2,d,L);
disp(res)
   7.0315e-08

Coeficientes Bn

Los coeficientes Bn se calculan a partir de las distribución inicial de temperatura u(x,0)

u(x,0){ u 1 (x,0)= n=0 B n X 1n ,0x<d u 2 (x,0)= n=0 B n X 2n ,dxL

Multiplicamos la primera ecuación por X1m y la segunda por X2m e integramos entre 0 y d, la primera y entre d y L, la segunda

{ 0 d u 1 (x,0) X 1m dx = 0 d ( n=0 B n X 1n X 1m ) dx d L u 2 (x,0) X 2m dx = d L ( n=0 B n X 2n X 2m ) dx { 0 d u 1 (x,0) X 1m dx = n=0 B n · 0 d ( X 1n X 1m dx ) dx d L u 2 (x,0) X 2m dx = n=0 B n · d L ( X 2n X 2m dx ) dx

Sumamos

ρ 1 c 1 0 d u 1 (x,0) X 1m dx + ρ 2 c 2 d L u 2 (x,0) X 2m dx = n=0 B n ·{ ρ 1 c 1 0 d ( X 1n X 1m dx ) dx+ ρ 2 c 2 d L ( X 2n X 2m dx ) dx } ρ 1 c 1 0 d u 1 (x,0) X 1m dx + ρ 2 c 2 d L u 2 (x,0) X 2m dx = B m { ρ 1 c 1 0 d X 1m 2 dx+ ρ 2 c 2 d L X 2m 2 dx }

El primer coeficiente B0, correspondiente a ω0=0, merece un cálculo aparte. X10=1, X20=1

ρ 1 c 1 0 d u 1 (x,0) X 10 dx + ρ 2 c 2 d L u 2 (x,0) X 20 dx = B 0 { ρ 1 c 1 0 d X 10 2 dx+ ρ 2 c 2 d L X 20 2 dx } ρ 1 c 1 0 d ( T 01 T eq )dx + ρ 2 c 2 d L ( T 02 T eq )dx = B 0 { ρ 1 c 1 0 d dx + ρ 2 c 2 d L dx } ρ 1 c 1 ( T 01 T eq )d+ ρ 2 c 2 ( T 02 T eq )(Ld)= B 0 ( ρ 1 c 1 d+ ρ 2 c 2 (Ld) )

Introduciendo la expresión de la temperatura final de equilibrio, Teq. El resultado es B0=0

Los otros coeficientes Bn son

ρ 1 c 1 0 d ( T 01 T eq )cos( μ 1 ω m x)dx + ρ 2 c 2 d L ( T 02 T eq ) cos( μ 1 ω m d) cos( μ 2 ω m (Ld) ) cos( μ 2 ω m (Lx) )dx = B m { ρ 1 c 1 0 d X 1m 2 dx+ ρ 2 c 2 d L X 2m 2 dx } ρ 1 c 1 μ 1 ω m ( T 01 T eq )sin( μ 1 ω m d)+ ρ 2 c 2 μ 2 ω m ( T 02 T eq ) cos( μ 1 ω m d) cos( μ 2 ω m (Ld) ) sin( μ 2 ω m (Ld) )= B m { ρ 1 c 1 0 d X 1m 2 dx+ ρ 2 c 2 d L X 2m 2 dx }

Utilizaremos procedimientos numéricos para calcular la expresión entre corchetes del segundo miembro de la igualdad.

Una vez calculados los coeficientes Bn, se ha completado la función u(x, t). Ahora queda la temperatura de la varilla, T(x, t)=u(x, t)+T(x, ∞)=u(x, t)+Teq

T 1 (x,t)= n=1 B n cos( μ 1 ω n x)exp( ω n 2 t) + T eq ,0x<d T 2 (x,t)= n=1 B n cos( μ 1 ω n d) cos( μ 2 ω n (Ld) ) cos( μ 2 ω n (Lx) )exp( ω n 2 t) + T eq ,dxL

Dado que la exponencial disminuye rápidamente con el tiempo y por otra parte, los valores de las raíces ωn se hacen grandes, basta sumar unos pocos términos a la serie para obtener buenos resultados.

Completamos el programa para representar la distribución de temperaturas en los instantes t: 50, 500, 2000, 10000 s. La línea horizontal a trazos corresponden al estado estacionario. La parte derecha de la varilla, está hecha de acero y la parte izquierda, de aluminio de conductividad térmica K mucho mayor

k1=209.3; %aluminio
c1=880;
r1=2700;
k2=45; %acero
c2=460;
r2=7800;
a1=sqrt(c1*r1/k1);
a2=sqrt(c2*r2/k2);

T01=100; %temperaturas iniciales
T02=0;
L=1; %longitud varilla
d=L/3; %posición contacto

%estado estacionario, temperatura de equilibrio
Teq=(r1*d*c1*T01+r2*(L-d)*c2*T02)/(r1*d*c1+r2*(L-d)*c2);

%raíces wn
w=linspace(0,0.2,50);
f=@(w) sin(a1*w*d).*cos(a2*w*(L-d))+(k2*a2/(k1*a1))*sin(a2*w*(L-d)).*cos(a1*w*d);
raiz=raices(f,w); 

%amplitudes An
A=zeros(1,length(raiz));
for n=1:length(raiz)
    w=raiz(n);
     X1=@(x) cos(a1*w*x);
     X2=@(x) cos(a1*w*d)*cos(a2*w*(L-x))/cos(a2*w*(L-d));
    Z1=@(x) X1(x).^2;
    Z2=@(x) X2(x).^2;
    den=integral(Z1,0,d)+(r2*c2/(r1*c1))*integral(Z2,d,L);
    num=(T01-Teq)*sin(a1*w*d)/(a1*w)+(r2*c2/(r1*c1))*(T02-Teq)*
cos(a1*w*d)*sin(a2*w*(L-d))/(a2*w*cos(a2*w*(L-d)));
    A(n)=num/den; 
end
hold on
line([0,L],[Teq,Teq],'lineStyle','--','color','k')
line([d,d],[-20,100],'lineStyle','--','color','k')
x1=linspace(0,d,50);
x2=linspace(d,L,100);
 for t=[50,500,2000, 10000]
    u1=zeros(1,length(x1));
    u2=zeros(1,length(x2));
    for n=1:length(raiz)
        w=raiz(n);
        u1=u1+A(n)*cos(a1*w*x1)*exp(-w^2*t);
        u2=u2+A(n)*cos(a1*w*d)*cos(a2*w*(L-x2))*exp(-w^2*t)/cos(a2*w*(L-d));
    end
    T=[u1+Teq,u2+Teq];
    x=[x1,x2];
    plot(x,T)
end
hold off
xlabel('x')
ylabel('T')
grid on
title('Conducción del calor, no homogéneo')