Diferencia entre revisiones de «Modelo Térmico. (grupo 6C)»

Revisión actual del 21:35 12 mar 2015

Este artículo está en versión beta. El autor de este artículo no lo ha terminado todavía, por favor no lo edites hasta que elimine este mensaje.

1 INTRODUCCIÓN

El trabajo que a continuación desarrollaremos, se basa en la formulación de un modelo matemático que describa el comportamiento de la temperatura en el interior de un edificio a lo largo de 24 horas en función de:

• Temperatura exterior. (M(t))
• Calor que se genera dentro del edificio. (H(t))
• Enfriamiento y calentamiento producido por el sistema de aire acondicionado y calefacción respectivamente. (U(t))

Para formular el modelo hemos utilizado un análisis compartimental ya que consideramos el edificio como una sola habitación en el que la temperatura está definida como T(t)

2 PROBLEMA DE CAUCHY

La ley de enfriamiento de Newton establece que:

"Cuando la diferencia entre un cuerpo y su medio externo no es demasiado grande, el calor transferido en la unidad de tiempo hacia el cuerpo o desde el cuerpo por conducción, convección y radiación es aproximadamente proporcional a la diferencia de temperatura entre el cuerpo y el medio externo"

Teniendo esto en cuenta, esta ley establece que hay una razón de cambio de temperatura que es proporcional a la temperatura interior y exterior del edificio designadas en nuestro caso como T(t) y M(t) respectivamente. Dicha razón queda determinada según el P.V.I siguiente:

[math] \left \{ \begin{matrix} T'=k(M(t)-T)+H(t)+U(t), t\gt0 \\ T(to)=T_0 \end{matrix} \right . [/math]

k: es una constante real que depende del número de ventanas, aislamiento del edificio..

3 ANÁLISIS DE LA TEMPERATURA EN EL INTERIOR DEL EDIFICIO

3.1 EDIFICIO "CERRADO"

A la hora de analizar la temperatura del edificio vamos a ir suponiendo distintas situaciones,en está primera situación vamos a tratar de saber cuanto tiempo tarda en cambiar considerablemente la temperatura del edificio cuando esta cerrado, esto quiere decir que en el interior del edificio no hay actividad que afecte al calentamiento o enfriamiento del mismo ni tampoco hay calefacción ni aire acondicionado, por tanto las variables U(t) y H(t) son iguales a cero; por sólo vamos a tener en cuenta la temperatura exterior M(t).

Suponemos que la temperatura en el interior del edificio y a media noche es 14ºC, y la temperatura exterior es constante e igual 8ºC, (M(t)=8ºC), por otra parte, la constante de tiempo del edificio está definida como t'= [math]\dfrac{1}{k} [/math], nosotros suponemos que t'=3, por lo tanto la constante k=[math]\dfrac{1}{3} [/math].

Con las hipótesis expuestas en el párrafo anterior, obtenemos el siguiente P.V.I:

[math] \left \{ \begin{matrix} T'=\dfrac{1}{3}(8-T), t\gt0 \\ T(to)=14 \end{matrix} \right . [/math]

Utilizamos el método de Euler implícito para el PVI anterior, con una longitud de paso h=0,001 y compararemos este resultado con el resultado de la ecuación exacta obtenida analíticamente. Implementado el código que a continuación se muestra en Matlab, se ha llegado al gráfico que responde a nuestra cuestión inicial.

%Variacion de la temperatura del edificio cerrado con Euler implícito
t0 = 0;                  
tn = 24; 
T0 = 14;                 
h = 0.01;                
N =(tn-t0)/h;            
t = linspace(0,24,N+1);  
%solución exacta 
 Te = (14-8).*exp(-1/3*t)+8;

T = zeros(size(t));     
T(1) = T0;               
for i = 1:N
  T(i+1) = 1/(1+h/3)*(T(i)+h*(8/3));
end
[T1,t1,Te1];
%Gráfica
hold on
ylabel('Temperatura'); xlabel('Tiempo(h)');
plot(t1,Te1,'g');
plot(t1,T1, 'r',);

hold off

centro

Ante todo cabe destacar que la solución exacta y la obtenida mediante el método numérico son prácticamente idénticas. Además en el gráfico se observa un brusco descenso de la temperatura durante las primeras horas del día, para luego estabilizarse cuando ya han pasado aproximadamente 10 horas, es decir el edificio trata de alcanzar el equilibrio térmico entre la temperatura exterior e interior perdiendo calor durante la noche.

3.2 EDIFICIO EN UN DÍA HÁBIL DURANTE LA PRIMAVERA U OTOÑO

En esta segunda situación vamos a analizar la variación de la temperatura en un día hábil, esto quiere decir que a diferencia del aparado anterior en el edificio se esta desarrollando una actividad que genera un calentamiento como por ejemplo las luces, máquinas... Dicho calentamiento está definido mediante la variable H(t). Al suponer estar en una estación del año intermedia no se precisa calefacción ni aire acondicionado por lo tanto, U(t)=0.Por otra parte la temperatura exterior M(t) variará de forma senoidal en un periodo de 24 horas con un mínimo en t=0h y un máximo en t=12h. Con Mo y B constantes >0 y [math]\ ω=\frac{π}{12}\ [/math] radianes/hora y un [math]\ T(0)=To [/math]. Para demostrar que: [math] T(t)=Bo - BF(t) + Ce^{-kt} [/math](1) siendo \[\left\{\begin{matrix}\ Bo = Mo + \frac{Ho}{k}\ , & \\ F(t)= \frac{\cos(ωt)+\frac{ω}{k}\sin(ωt)}{1+\frac{ω^2}{k^2}}\ \, \\ C= To - Bo + \frac{B}{1+\frac{ω^2}{k^2}}\\ & \end{matrix}\right.\] Hay que resolver el problema de valor inicial: \[\left\{\begin{matrix}\ T'=k[Mo-Bcos(ωt)-T]+Ho, \\ T(0)=To & \end{matrix}\right.\] Se trata de una ecuación lineal de primer orden en T. La solución será del tipo: [math]T(t)=T_{hom}(t)+T_{par}(t)[/math] La solución homogénea es: [math]T_{hom}=ce^{-kt}[/math] Para resolver la solución particular aplicamos el principio de la superposición y obtenemos dos soluciones particulares una para [math] y_{1}(t)=kMo+Ho [/math] y otra para [math] y_{2}(t)=-kBcos(ωt) [/math] Tanteamos una solucion [math] T_{p1}=Α \ y \ T'_{p1}=0 [/math] y resolviendo nos da [math] Α=Mo+\frac{Ho}{k}\ [/math] y [math] T_{p1}=Mo+\frac{Ho}{k}\ [/math] Y seguidamente para la segunda solucion particular tanteamos [math] T_{p2}=fcos(ωt)+gsin(ωt) [/math] que metiendola en la ecuación [math] T'_{p2}+kT_{p2}=-kBcos(ωt)[/math] y resolviendo por Cramer obtenemos [math] f=\frac{ωkB}{-{ω}^2-{k}^2}\ [/math] [math] g=\frac{{k}^2B}{-{ω}^2-{k}^2}\ [/math] y por tanto [math] T_{p2}=\frac{-B}{1+{(\frac{ω}{k})}^2}cos(ωt)-\frac{\frac{ω}{k}B}{1+{(\frac{ω}{k})}^2}sin(ωt)[/math]

La solución [math] T(t)=T_{h}+T_{p1}+T_{p2} [/math] nos queda [math] T(t)=Mo+\frac{Ho}{k}\ -B[\frac{1}{1+{(\frac{ω}{k})}^2}cos(ωt)+\frac{\frac{ω}{k}}{1+{(\frac{ω}{k})}^2}sin(ωt)] +Ce^{-kt} [/math] e imponiendo la condición inicial [math] T(0)=To [/math] se obtiene un [math] C=To -Bo+\frac{B}{1+{(\frac{ω}{k})}^2}\ [/math] exactamente como dice el enunciado.

DEMOSTRACIÓN DE QUE [math] B_0[/math] ES APROXIMADAMENTE IGUAL A [math] M_0[/math]

[math]TM(t)=\frac{\int_{0}^{24} T(t)\, dt}{24}=\frac{\int_{0}^{24} (Bo-B*F(t)+c*e^{-k*t})\, dt}{24}\simeq\frac{\int_{0}^{24} Bo\, dt}{24}\simeq Bo[/math]

donde [math]\int_{0}^{24} F(t)\, dt=0;:\int_{0}^{24} c*e^{-k*t}\, dt\simeq 0[/math].

Por último vamos a demostrar que la variación senoidal del edificio tiene un retraso de α horas respecto a la variación exterior con un [math] F(t) = (1+\frac{ω^2}{k^2})^{-1/2} cos(ωt-α), t\gt0 [/math] y que [math] tan(α)=\frac{ω}{k}\ [/math]

De modo que nos queda [math] F(t) = (1+\frac{ω^2}{k^2})^{-1/2}(cos(ωt) + \frac{ω}{k}\sin(ωt)) = (1+\frac{ω^2}{k^2})^{-1/2}cos(ωt-α) [/math] → [math] cos(ωt) + \frac{ω}{k}\sin(ωt) = cos(ωt)cos(α) + sin(ωt)sin(α) [/math]

\[\left\{\begin{matrix}\ cos(ωt)=cos(ωt)cos(α) & \frac{w}{k}\sin(ωt)=sin(ωt)sin(α) & \end{matrix}\right.\] \[\left\{\begin{matrix}\ 1=cos(α) & \frac{w}{k}=sin(α) & \end{matrix}\right.\]

Y por lo tanto [math] tan(α)=\frac{ω}{k}\ [/math]

Ahora con la ayuda de MATLAB,mediante el método de Euler y Runge-Kutta 4 resolvemos el P.V.I, que se nos presenta para los datos: [math]H_0=3, M_0=7, K=1/3, B=5, T_0=13[/math]

t0=0;
tN=24;
y0=13;
H0=3;
M0=7;
k=1/3;
B=5;
w=pi/12;
B0=M0+(H0/k);
C=y0-B0+(B/(1+(w/k)^2));
h=input('Introduce tamaño de paso');
N=(tN-t0)/h;

t=t0:h:tN;

F=(1/(1+(w/k)^2))*(cos(w*t)+(sin(w*t)*(w/k)));
%solucion exacta
x=B0-B*F+C*(exp(-k*t));

y=zeros(1,N+1); %euler
y(1)=y0;

u=zeros(1,N+1); %runge-kutta
u(1)=y0;


for i=1:N
 y(i+1)=y(i)+h*((k*((M0-(B*cos(w*t(i))))-y(i)))+3); %euler

 U1=(k*((M0-(B*cos(w*t(i))))-y(i)))+3;
 U2=(k*((M0-(B*cos(w*(t(i)+(0.5*h))))-(y(i)+(0.5*U1*h))))+3);
 U3=k*((M0-(B*cos(w*(t(i)+(0.5*h))))-(y(i)+(0.5*U2*h))))+3;
 U4=k*((M0-(B*cos(w*(t(i)+(h))))-(y(i)+(U3*h))))+3;
 u(i+1)=u(i)+h/6*(U1+2*U2+2*U3+U4); %runge-kutta4
end
M=M0-(B*cos((w)*t));
figure(1)
hold on
plot(t,y,'b')
plot(t,u,'g')
plot(t,M,'r')
plot(t,x,'y')
legend('Temperatura interior, Euler','Temperatura interior, Runge-kutta','Temperatura exterior')
xlabel('Tiempo en horas');
ylabel('Temperatura en ºC');
hold off

ANÁLISIS DEL GRÁFICO:

En el gráfico se observa que la temperatura media en un día es aproximadamente 16ºC ya que anteriormente hemos demostrado que Bo=Tm. Otro aspecto importante es el retraso de 3 horas de la temperatura interior respecto a la temperatura exterior, que da lugar a periodos en los que la temperatura esta disminuyendo dentro del edificio aumentando la de fuera.

3.3 EDIFICIO EN UN DÍA HÁBIL DE VERANO O INVIERNO

En esta situación vamos a añadir una variable más, aparte de la temperatura exterior y el calor generado por la actividad en el edificio, se trata de un termostato sencillo que regula la temperatura del interior del edificio en función de la temperatura deseada, designada por [math]T_D[/math]. El funcionamiento del termostato es el siguiente: el sistema mide la temperatura real en el interior del edificio y aporta frío o calor dependiendo de la estación, si la temperatura real y la deseada no coinciden, en caso de coincidir el termostato permanece apagado.

El calor o el frío suministrado viene dado por la función:

[math]U(t)=K_u[T_D-T(t)][/math]

Del enunciado sabemos que la solución analítica del P.V.I es:

[math] T(t) = B2-B1F1(t)+Ce^{-k1t} [/math]

Donde: \[\left\{\begin{matrix}\ K1 = k+Ku\ , & \\ B2 = \frac{KuT_{D}+kMo+Ho}{K1}\ \, \\ B1 = \frac{kB}{K1}\ \, \\ F1(t)= \frac{\cos(ωt)+\frac{ω}{k1}\sin(ωt)}{1+\frac{ω^2}{k^2}}\ \, \\ C= To - B2 + F1(0) \ & \end{matrix}\right.\] Para demostrar que lo anterior es cierto, resolvemos el siguiente P.V.I :

[math] T' = k[(Mo-Bcos(ωt)]-T]+Ho+Ku[T_{D}-T] [/math]

Como se puede observar, se trata de una ecuación de primer orden cuya solución será del tipo:

[math]T_G=T_H+T_P[/math]

* Solución de la ecuación homogénea:

[math]T'+KT + KuT = 0[/math] Aplicando el método de variables separadas tenemos que:

[math] T = Ce^{-(k+Ku)} [/math] 

y así se llega a la conclusión de que [math]K1= k+K_u[/math]

• Sólución particular:

Para hallarla hemos aplicado el principio de superposición y obtenemos dos soluciones particulares una para cada una de las siguientes ecuaciones:

• [math] y_{1}(t)=kMo+Ho+KuT_{D} (1) [/math]
• [math] y_{2}(t)=-kBcos(ωt) (2) [/math]

En (1) aplicamos el método de resolución por tanteo, y obtenemos lo siguiente:

[math] T_{P1}=Α \ y \ T'_{P1}=0 [/math] obtenemos [math] A = \frac{KuT_{D}+kMo+Ho}{k + Ku}\ [/math] y [math] T_{p1}= \frac{KuT_{D}+kMo+Ho}{K1}\ [/math].

Vemos que [math]T_{P1}=B2[/math] (q.c.d)

Para calcular la segunda solución particular: [math] T_{p2}=fcos(ωt)+gsin(ωt) [/math] . Derivando esta expresión [math] T_{P2}=-fsen(ωt)+gωcos (ωt) [/math] y sustituyéndola en la EDOL llegamos a un sistema de ecuaciones que se resuelve por Cramer, obteniendo los valores de f y g y sustituyéndolos en [math] T_{P2}[/math]

[math] T_{P2} (t)=-\frac{\frac{kB}{K_{1}}}{1+{(\frac{ω}{K_{1}})}^2}cos(ωt)-\frac{\frac{KB}{K_{1}^2}ω}{1+{(\frac{ω}{K_{1}})}^2}sin(ωt) [/math]

La solución es [math] T(t)=T_{H}+T_{P1}+T_{P2} [/math] donde: [math]T(t)=\frac{kMo+Ho+kUT_{D}}{K_{1}}-B1{[\frac{\frac{k}{K_{1}}}{1+{(\frac{ω}{K_{1}})}^2}cos(ωt)+\frac{\frac{k}{K_{1}^2}ω}{1+{(\frac{ω}{k_{1}})}^2}sin(ωt)]}+Ce^{-kt}[/math] Idenfificando coeficientes observamos que el coeficiente que multiplica a B1 es F1(t). (c.q.d)

Por último debemos demostrar que [math] C = To – B2 + +F1(0) [/math], para lograrlo basta con imponer la condición: [math]T(0)=T_0[/math].

Por otro lado vamos a demostrar que B2 representa aproximadamente la temperatura media del edificio, muy próxima a la temperatura deseada. Primero suponemos suponemos un valor medio de 24 horas. La temperatura media esta definida como

[math]T_{media}(t)=\frac{\int_{0}^{24} T(t)\, dt}{24}=\frac{\int_{0}^{24} (B_{2}-B_{1}*F_{1}(t)+c*e^{-K_{1}*t})\, dt}{24}\simeq\frac{\int_{0}^{24} B_{2}\, dt}{24}\simeq B_{2}[/math]

Siendo:

• [math]\int_{0}^{24} B_{2}*F_{1}(t)\, dt=0[/math] Esto es así porque las funciones son periodicas

• [math]\int_{0}^{24} c*e^{-K_{1}*t}\, dt\simeq 0[/math] Debido a que el término exponecial nos dice el enunciado que desaparece a los 30 minutos.

[math]B2=\dfrac{K_uT_D+KM_O+H_O}{K1}=\dfrac{K_uT_D+KM_O+H_O}{K+K_u}=\dfrac{K_uT_D}{K+K_u}+\dfrac{KM_o+H_o}{K+K_u}=T_D[/math]

ya que [math]\dfrac{K_u}{K+K_u}[/math] es aprox 1 y [math]\dfrac{KM_O+H_O}{K+K_u}[/math] es aproximadamente 0.

Con esto hemos conseguido demostrar que B2 y Tm es aproximadamente la temperatura deseada.