Modelo térmico aplicado a un edificio. Grupo 8C

Trabajo realizado por estudiantes
Título	Modelo térmico aplicado a un edificio. Grupo 8-C
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2014-15
Autores	Miguel Arbaizar, Cristino Pérez , Javier Mellado , Alejandro López
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

A lo largo de este trabajo, formularemos un modelo matemático para estudiar el comportamiento de la temperatura en el interior de un edificio durante un periodo de 24 horas. La temperatura dependerá de las siguientes variables:

• Temperatura exterior [M(t)]
• Calor producido en el interior por las personas, luces, máquinas, etc. [H(t)]
• Calentamiento y enfriamiento por la calefacción y/o aire acondicionado [U(t)]

Utilizaremos un análisis compartimental según el cual consideraremos a nuestro edificio como un solo compartimento en el que su temperatura estará designada por T(t).

2 Problema de valor inicial (P.V.I)

La razón de cambio de la temperatura quedará determinada según el problema de valor inicial o de Cauchy:

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

Donde k es una constante real que depende de las condiciones del edificio tales como el número de ventanas, aislamiento, etc.

3 Estudio de la temperatura interior

3.1 En función exclusivamente de la temperatura exterior [M(t)]

En este caso vamos a considerar que al final del día (to=0), la temperatura exterior permanece constante (M=8), y tanto el calor producido en el interior como el aportado por la calefacción/aire acondicionado es nulo (H=U=0). La temperatura para to=0 (medianoche) será de 14 grados ºC.

La constante de tiempo del edificio nos da una medida de cuánto tarda en cambiar considerablemente la temperatura del edificio. Esta constante también es conocida como la constante de tiempo de transferencia de calor entre el edificio y el exterior. La expresión será [math]k=\frac{1}{t}[/math], que en esta primera hipótesis será [math]\frac{1}{3}[/math].

El problema de valor inicial (P.V.I) en este caso será:

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

3.1.1 Método de Euler Implícito

Con la ayuda de Matlab, vamos a resolver el P.V.I con un tamaño de paso h=0.01

centro

Como podemos observar, la temperatura del edificio parte de 14º, para después bajar rapidamente (en las 5 primeras horas del día aproximadamente) hasta estabilizarse en 8º. Es decir, nuestro edificio pierde a lo largo de la noche su calor interior hasta igualarse a la temperatura exterior.

Tambien es fácil comprobar que a mayor diferencia de temperatura entre el interior y exterior del edificio, más rápido cambiará la función T. Cuando las dos temperaturas comienzan a ser similares, alrededor de las 5 de la mañana, la temperatura interior tiende a 8º y la velocidad de enfriamiento se reduce a 0 lentamente, a lo largo del resto del día.

3.2 En función de la temperatura exterior [M(t)] y el calor producido en el interior [H(t)]

3.2.1 Solución teórico-analítica

3.2.1.1 Demostración de la solución del P.V.I.

En este caso vamos a resolver analíticamente el P.V.I. siguiente: [math] \left \{ \begin{matrix} T'=k(Mo-B*cos(wt)-T)+Ho, t\gt0 \\ T(to)=T_0 \end{matrix} \right . [/math]

Operando la ecuación diferencial obtenemos: [math]T'+k*T=k*Mo-K*B*cos(wt)+Ho[/math]Se trata de una ecuación lineal de primer orden en T. La solución de este problema de valor inicial será del tipo: [math]T(t)=T_{H}(t)+T_{P}(t)[/math] En un primer lugar calculamos la homogénea: [math]T'+k*T=0 \rightarrow \frac{dT}{T}=-k*dt \rightarrow T_{H}=c*e^{-k*t}[/math] En segundo lugar calculamos la solución particular, que se realizará aplicando el principio de superposición de las soluciones: [math]f_{1}(t)=k*Mo+Ho; f_{2}(t)=-k*B*cos(wt) \rightarrow T_{P}(t)=T_{P_1}(t)+T_{P_2}(t)[/math] La solución para la primera particular será, sustituyendo en la completa: [math]T_{P_1}=A; T'_{P_1}=0 \rightarrow T'_{P_1}+k*T_{P_1}=k*Mo+Ho \rightarrow k*A=k*Mo+Ho \rightarrow A=Mo+\frac{Ho}{k} \rightarrow T_{P_1}(t)=Mo+\frac{Ho}{k}[/math] La solución para la segunda particular se obtendrá tanteando con: [math]T_{P_2}(t)=a*cos(wt)+b*sen(wt); T'_{P_2}(t)=-a*w*sen(wt)+b*w*cos(wt)[/math]Sustituyendo en la completa: [math]-a*w*sen(wt)+b*w*cos(wt)+k*a*cos(wt)+k*b*sen(wt)=-k*B*cos(wt)[/math]Identificando coeficientes obtenemos el sistema: [math] \left \{ \begin{matrix} b*w+k*a=-k*B \\ k*b-a*w=0 \end{matrix} \right . [/math] Del sistema, aplicando Cramer, obtenemos: [math]b=\frac{\begin{vmatrix} -k*B & k \\ 0 & -w \end{vmatrix}}{\begin{vmatrix} w & k \\ k & -w \end{vmatrix}}=\frac{w*k*B}{-w^2-k^2}=\frac{\frac{w}{k}*B}{1+{(\frac{w}{k})}^2};a=\frac{\begin{vmatrix} w & -k*B \\ k & 0 \end{vmatrix}}{\begin{vmatrix} w & k \\ k & -w \end{vmatrix}}=\frac{k^2*B}{-w^2-k^2}=\frac{-B}{1+{(\frac{w}{k})}^2}[/math] La segunda solucion particular será: [math]T_{P_2}(t)=\frac{-B}{1+{(\frac{w}{k})}^2}*cos(wt)-\frac{\frac{w}{k}*B}{1+{(\frac{w}{k})}^2}*sin(wt)[/math] Finalmente la solución completa será: [math]T(t)=Mo+\frac{Ho}{k}-B*\frac{1}{1+{(\frac{w}{k})}^2}*cos(wt)+\frac{\frac{w}{k}}{1+{(\frac{w}{k})}^2}*sin(wt)+c*e^{-k*t}[/math] Identificando los coeficientes obtenidos con los de la solución propuesta obtenemos los valores que pretendíamos demostrar:[math]Bo=Mo+\frac{Ho}{k};F(t)=\frac{1}{1+{(\frac{w}{k})}^2}*cos(wt)+\frac{\frac{w}{k}}{1+{(\frac{w}{k})}^2}*sin(wt)[/math] Sustituyendo la condición inicial en la completa (T(0)=T_0) se obtiene el último coeficiente pedido:[math]T_0=Mo+\frac{Ho}{k}-B*\frac{1}{1+{(\frac{w}{k})}^2}*cos(w*0)+\frac{\frac{w}{k}}{1+{(\frac{w}{k})}^2}*sin(w*0)+c*e^{-k*0} \rightarrow c=T_0-Bo+\frac{B}{1+{(\frac{w}{k})}^2}[/math]

3.2.1.2 Demostración de Bo es aproximadamente igual a Mo

La temperatura media del edificio queda definida como: [math]Tmedia(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] Siendo:[math]\int_{0}^{24} F(t)\, dt=0;:\int_{0}^{24} c*e^{-k*t}\, dt=0[/math] Ya que la primera de ellas está formada por funciones periódicas y la segunda indica el enunciado que es nula.

3.2.1.3 Demostración de la variación senoidal del edificio tiene un retraso igual a [math]\boldsymbol{\alpha}[/math]

Si F(t)=(1+(\frac{w}{k})}^2)^{-\frac{1}{2}}*(cos(wt)+\frac{w}{k}*sin(wt))=(1+(\frac{w}{k})}^2)^{-\frac{1}{2}}*cos{(wt-\boldsymbol{\alpha})}</math>

Haciendo el desarrollo de [math]cos(wt)+\frac{w}{k}*sin(wt)=cos(wt)*cos(\boldsymbol{\alpha})+sin(wt)*sin(\boldsymbol{\alpha})[/math] llegamos a la igualdad:

[math] \left \{ \begin{matrix} cos(\boldsymbol{\alpha})=1 (1) \\ sen(\boldsymbol{\alpha})=\frac{w}{k} (2) \end{matrix} \right . [/math] De donde se sigue, si dividimos [math]\frac{(2)}{(1)}[/math] obtenemos que: [math]\frac{sen(\boldsymbol{\alpha})}{cos(\boldsymbol{\alpha})}=tan(\boldsymbol{\alpha})=\frac{w}{k}[/math]

3.2.2 Solución numérica

En este ocasión vamos a considerar que sí existe una razón de calentamiento adicional [H(t)], que supondremos constante a lo largo del día, a H(t)=Ho=3. De la misma manera que en el apartado anterior, el calor aportado por la calefacción/aire acondicionado es nulo (U=0), es decir, nos vamos a adentrar en el estudio de la temperatura del edificio durante las estaciones de primavera y otoño, cuando no se precisan calefactores o aparatos de aire acondicionado para alcanzar la temperatura deseada. La constante de tiempo será [math]k=\frac{1}{3}[/math] y la temperatura para to=0 (medianoche) será de 13 grados ºC.

A diferencia del caso anterior, la temperatura exterior [M(t)] no permanece constante, sino que varía en forma de onda senoidal durante un periodo de 24 horas, con su mínimo en t=0 (medianoche) y su máximo en t=12 (mediodía), de la siguiente manera: [math]M(t)=Mo-B*cos(wt)[/math]En el caso que nos ocupa, Mo recibirá el valor de 7ºC, B el valor de 5ºC y w el valor de [math]\frac{2*π}{24}=\frac{π}{12}[/math]radianes/hora. La expresión de la temperatura exterior quedará definida de la siguiente forma: [math]M(t)=7 - 5*cos(\frac{π}{12}*t)[/math]

A partir de la situación anteriormente descrita, planteamos el problema de valor inicial (P.V.I):

[math] \left \{ \begin{matrix} T'=\frac{1}{3}*(7-5*cos(\frac{π}{12}*t)-T)+3, t\gt0\\ T(0)=13 \end{matrix} \right . [/math]

3.2.2.1 Método de Euler

Con la ayuda de Matlab, vamos a resolver el P.V.I con un tamaño de paso h1=0.1, h2=0.01 y h3=0.001. Sin embargo, por la baja resolución de la gráfica adjunta y las pequeñas diferencias milimétricas entre las funciones obtenidas para cada tamaño de paso, apenas son apreciables las diferencias entre ellas.

centro

Mediante esta gráfica, podemos extraer varios aspectos interesantes de este estudio. Primeramente, podríamos decir que la temperatura exterior esta actuando sobre el edificio como modulador de su temperatura, aumentando y disminuyendo ésta de forma similar, con variaciones muy parecidas. Es importante señalar que la expresión del P.V.I. obtenida confiere a la función de la temperatura interior T un comportamiento más "equilibrado" que la temperatura exterior, pues ambas fluctúan alrededor de los 16ºC y los 7ºC respectivamente, pero abarcando un rango de temperaturas menor : Tε[12,20] y Mε[2,12].

Sin embargo, los periodos de variación de ambas funciones son iguales, de 12 horas exactamente. Tambien, debemos mencionar el retraso de α horas de la temperatura interior respecto a la exterior, que conllevan pequeños periodos durante los que la temperatura esta disminuyendo dentro del edificio y aumentando afuera. (Este curioso fenómeno producido se trata con mayor profundidad en la parte analítica de este trabajo).

3.2.2.2 Método de Runge-Kutta de orden 4

Con la ayuda de Matlab, vamos a resolver el P.V.I con un tamaño de paso h1=0.1, h2=0.01 y h3=0.001.

centro

3.3 En función de la temperatura exterior [M(t)], el calor producido en el interior [H(t)] y el efecto de la calefacción y/o aire acondicionado [U(t)]

3.3.2 Solución numérica

A partir de este momento el estudio de la variación de la temperatura en el interior del edificio lo realizaremos suponiendo que está afectado por todas las variables. Existe una razón de calentamiento adicional [H(t)], que seguimos suponiendo constante (H=3). La temperatura para t0=0 (medianoche) será de 24 grados ºC.

La temperatura exterior [M(t)] sigue variando de forma de senoidal durante un periodo de 24 horas, con su mínimo en t=0 (medianoche) y su máximo en t=12 (mediodía), de la siguiente manera: [math]M(t)=Mo-B*cos(wt)[/math]En esta situación, Mo recibirá el valor de 4ºC, B el valor de 5ºC y w el valor de [math]\frac{2*π}{24}=\frac{π}{12}[/math]radianes/hora. La expresión de la temperatura exterior quedará definida de la siguiente forma: [math]M(t)=4 - 5*cos(\frac{π}{12}*t)[/math]

Suponemos que en el edificio del apartado anterior está instalado un termostato sencillo que se utiliza para comparar la temperatura real interior del edificio con una temperatura deseada TD. Si la temperatura es menor que la deseada, el calefactor proporciona calentamiento; si es mayor se mantiene apagado. Si la temperatura es mayor que la deseada, el airea acondicionado proporciona enfriamiento; en caso contrario se encuentra apagado. Suponemos que la cantidad de calor o enfriamiento suministrado es proporcional a la diferencia de temperatura U(t) = Ku [TD − T(t)] , donde Ku es una constante positiva. En esta situación la constante Ku tiene un valor de [math]k=\frac{7}{8}[/math] y la temperatura deseada es de 22 grados ºC.

Por lo tanto el P.V.I. que modeliza el esta situación será:

[math] \left \{ \begin{matrix} T'=\frac{1}{3}*(4-5*cos(\frac{π}{12}*t)-T)+3+\frac{7}{8}*(22-T), t\gt0\\ T(0)=24 \end{matrix} \right . [/math]

3.3.2.1 Método de Euler

Con la ayuda de Matlab, vamos a resolver el P.V.I con un tamaño de paso h1=0.1, h2=0.01 y h3=0.001.

centro

3.3.2.2 Método de Runge-Kutta de orden 4

Con la ayuda de Matlab, vamos a resolver el P.V.I con un tamaño de paso h1=0.1, h2=0.01 y h3=0.001.

centro

4 Estudio de la temperatura interior del edificio separado en dos zonas, A y B

Ahora vamos a suponer que el edificio está dividido en dos zonas, zona A y zona B. La temperatura en la zona A estará designada por [Ta(t)] en tanto que la temperatura de la zona B estará designada por [Tb(t)]. La temperatura inicial en la zona A será de 20 grados ºC mientras que la temperatura inicial de la zona B será de 18 grados ºC

En la zona A la razón de calentamiento adicional será [Ha(t)] y el efecto de la calefacción y/o aire acondicionado [Ua(t)]. En la zona B la razón de calentamiento adicional será [Hb(t)] y el efecto de la calefacción y/o aire acondicionado [Ub(t)]. En este caso, para las zonas A y B respectivamente: [math] \left \{ \begin{matrix} Ha(t)=10*\frac{t}{1+t} \\ Ua(t)=\frac{5}{3}*(22-Ta(t)) \end{matrix} \right . [/math]

[math] \left \{ \begin{matrix} Hb(t)=4*cos(\frac{π}{6}*t)\\ Ub(t)=\frac{13}{7}*(23-Tb(t)) \end{matrix} \right . [/math]

La temperatura exterior variará en forma de onda senoidal durante un periodo de 24 horas, con su mínimo en t=0 (medianoche) y su máximo en t=12 (mediodía), de la siguiente manera tanto en la zona A como B: [math]M(t)=Mo-B*cos(wt)[/math]En esta situación será igual tanto para la zona A como para la zona B y su valor será: [math]M(t)=2 - 7*cos(\frac{π}{12}*t)[/math]

La constante de tiempo de cada zona nos da una medida de cuánto tarda en cambiar considerablemente la temperatura de dicha zona. Esta constante también es conocida como la constante de tiempo de transferencia de calor entre el edificio y el exterior. Esta constante tendrá un valor de [math]\frac{1}{4}[/math] entre la zona A y el exterior, un valor de [math]\frac{1}{5}[/math] entre la zona B y el exterior y un valor de [math]\frac{1}{2}[/math] entre la zona A y la zona B.

Para la resolución de este caso, tendremos que resolver un sistema de ecuaciones diferenciales, que será el siguiente:

[math] \left \{ \begin{matrix} Ta'= \frac{1}{4}*(2-7*cos(\frac{π}{12}*t)-Ta) + 10*\frac{t}{1+t} + \frac{5}{3}*(22-Ta) + \frac{1}{2}*(Tb-Ta) \\ Tb'= \frac{1}{5}*(2-7*cos(\frac{π}{12}*t)-Tb) + 4*cos(\frac{π}{6}*t) + \frac{13}{7}*(23-Tb) - \frac{1}{2}*(Tb-Ta) \end{matrix} \right . [/math]

4.1 Método de Euler

Con la ayuda de Matlab, vamos a resolver el P.V.I con un tamaño de paso h1=0.1, h2=0.01 y h3=0.001.

centro

4.2 Método de Runge-Kutta de orden 4

Con la ayuda de Matlab, vamos a resolver el P.V.I con un tamaño de paso h1=0.1, h2=0.01 y h3=0.001.

centro

4.3 Método de Euler Implícito

Con la ayuda de Matlab, vamos a resolver el P.V.I con un tamaño de paso h1=0.1, h2=0.01 y h3=0.001.

centro