Ver el código de «Modelo Térmico. (grupo 6C)»

{{Beta}}
{{ TrabajoED | Modelo térmico. Grupo 22C | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] |Antonio Diaz Margarit, Pablo Rodríguez Sandoval, Marta Serrano Grande, David Sánchez Carretero, Ricardo García Dengra}}
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]

==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)

==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>0 \\
         T(to)=T_0
      \end{matrix}
   \right . </math>

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

==ANÁLISIS DE LA TEMPERATURA EN EL INTERIOR DEL EDIFICIO==
===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>0 \\
         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.

{{matlab|codigo=
%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, 'b',);
legend('Solucion exacta','Euler implicito',)
hold off

}}



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.

===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>0 </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>


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.

=== 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>