Ver el código de «Modelo Térmico de un Edificio.(Grupo 13A)»


{{ TrabajoED | Modelo Térmico de un Edificio.(Grupo 13A) | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | María Aguilera, Paula Martínez, Miguel Sánchez, Laura García, Isabel Roselló, Sarah Boufounas }}
[[Categoría:Ecuaciones Diferenciales]]
[[Categoría:ED14/15]]
[[Categoría:Trabajos 2014-15]]

==Introducción==
Este trabajo consiste en estudiar, mediante un modelo matématico, el comportamiento de la temperatura interior de un edificio en un periodo de 24 horas. Esta temperatura dependerá del calor generado dentro, de la temperatura exterior y del sistema de calefacción o aire acondicionado. Designaremos estos factores como:
:* H(t) → Calor generado por las personas, luces, máquinas...
:* M(t) → Temperatura exterior
:* U(t) → Calentamiento producido por la calefacción o enfriamiento debido al aire acondicionado

==Problema de Cauchy==
Designando la temperatura del edificio como T(t) el problema de valor inicial queda como:
[[Archivo:Jeje2.PNG|centro]]
Donde k es una constante del edificio que depende del número de puertas, ventanas, aislamiento, etc.
===Problema de Cauchy particularizado===
Denominamos t'=1/k al tiempo en que el edificio sufre un cambio de temperatura considerable, este tiempo es denominado como la constante de tiempo del edificio o constante de tiempo de transferencia de calor entre el edificio y el exterior.
Respondiendo a la primera pregunta el tiempo que tarda en cambiar considerablemente la temperatura para un edificio cerrado esta entre 2 y 4 horas.
En este caso particular tomaremos un tiempo t'=3 por lo que la constante k es igual a ⅓, el tiempo a medianoche es t=0 con una To=14, la temperatura exterior M(t) es constante e igual a 8 ºC y el calor generado por las personas H(t) y por la calefaccion U(t) es igual a 0, y nos queda el problema siguiente:
[[Archivo:Jeje3.PNG|centro]]
===Euler Implícito para el problema de valor inicial===
Tomamos una longitud de paso h=0.01 y comparamos el resultado con la ecuación exacta que resolviendo el problema analíticamente da T=6*e^(-1/3*t)+8 y podemos comprobar que es prácticamente idéntica.
[[Archivo:Jojo2.JPG|centro]]
[[Archivo:Jojo1.JPG|centro]]

==Modelo con una variación senoidal==
En este apartado consideramos la razon de calentamiento  <math>\ H(t)=Ho=cte </math>, no hay calefaccion o enfriamiento U(t)=0, y la temperatura exterior M(t) varia de forma senoidal en un perido de 24 horas con un mínimo en t=0 horas y un máximo en t=12 horas. <math>\ M(t)=Mo - Bcos(ωt) </math> 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> 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.\] 
Tenemos que resolver el problema de valor inicial: \[\left\{\begin{matrix}\ T'=k[Mo-Bcos(ωt)-T]+Ho, \\ T(0)=To & \end{matrix}\right.\]
Es una ecuación lineal de primer orden y la resolvemos como <math> Tgeneral=Thomogenea+Tparticular </math>.
:Para la homogénea <math> T'+kT=0 </math> la solución es <math> T_{h}=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 solucion <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.

==Apartado3b==
A continuación vamos a demostrar que la temperatura media es igual a B2

Suponemos un valor medio de 24 horas:
[[Archivo:Susu1.PNG|centro]]
[[Archivo:Susu3.PNG|centro]]

Seguidamente el ejercio nos pide que demostremos que la temperatura media está muy próxima a la temperatura deseada.
[[Archivo:Susu4.PNG|centro]]

==Estudio de la temperatura del edificio dividido en dos zonas==

El edificio consta ahora de dos zonas A y B. Designaremos a la temperatura en la zona A por Ta(t), y la temperatura en la zona B por Tb(t). Como dato del enunciado sabemos que la temperatura inicial en la zona A es de 20 grados y  la de la zona B es de 18 grados.
Ha(t) designara a la razón de calentamiento adicional en la zona A y  Ua(t) el efecto de la calefacción y/o aire acondicionado.
Asi mismo, Hb(t) designara a la razón de calentamiento adicional en la zona B  y  Ub(t) el efecto de la calefacción y/o aire acondicionado.

Sabemos que las constantes de tiempo de transferencia de calor son de 4 h entre A y el exterior, 5 h entre B y el exterior y 2 h entre A y B. 
Despejando de t=1/k obtenemos : 1/4 , 1/5 y 1/2 respectivamente.

Recordamos que el problema de valor inicial (PVI) o problema de Cauchy es el siguiente sistema de dos ecuaciones:

[[Archivo:Mariaecus2.PNG]]
 
La temperatura exterior  variara en forma de onda senoidal durante un t[0,24] de este modo :
M(t) - Mo - B * cos ( wt ) y nos queda que : M(t) - 2- 7 * cos (( pi/12)*t)

[[Archivo:Mariaecus3.PNG]]

En nuestro caso, el problema de Cauchy que satisface las temperaturas Ta(t) y Tb(t) es:

[[Archivo:Mariaecus1.PNG]]

===Solución numérica===

Para resolver el sistema de ecuaciones resultante nos ayudaremos de matlab. Se ha creado un programa que nos dará las soluciones a este sistema y nos dibujará las gráficas de las temperaturas (temperatura interior en zona A, temperatura interior en zona B y temperatura exterior) en función del tiempo durante 24 horas. Para ello deberemos dar un número de paso (h). En nuestro caso se ha utilizado h=0.1;h=0.01;h=0.001. Cuanto mayor sea el número de paso, mayor será la exactitud de la solución pues más se acercará esta a la solución exacta, la cual se lograría integrando. Se utilizan tres métodos distintos: Euler, Euler implícito y Runge-Kutta de orden 4. Cada uno de ellos nos dará una solución distinta, pues todos son métodos aproximados, aunque todos ellos nos darán soluciones muy similares, que serán también similares a la solución exacta.
Por último cabe señalar que al ser un sistema en lugar de una sola ecuación diferencial estos programas deberán de operan con matrices.

{{matlab|codigo=
%Datos del problema
t0=0; tN=24; h=0.001; y0=[20;18]; %Variaremos h dándole como valores 0.1;0.01;0.001
%Definimos la variable independiente
t=t0:h:tN ;
%Preparamos la matriz de la solución
x=zeros(2,length(t)); %Primera fila Temperatura A, segunda fila Temperatura B, Euler
%Inicializamos
x(:,1)=y0;
y=zeros(2,length(t)); %Primera fila Temperatura A, segunda fila Temperatura B,  Euler modificado
%Inicializamos
y(:,1)=y0;
z=zeros(2,length(t)); %Primera fila Temperatura A, segunda fila Temperatura B, Runge Kutta
%Inicializamos
z(:,1)=y0;
%Aplicamos el método
for i =1:length(t)-1
%Comenzacmos con Euler
 x(:,i+1)=x(:,i)+h*[(1/4)*(2-7*cos((pi/12)*t(i))-x(1,i))+((10*t(i))/(1+t(i)))+(5/3)*(22-x(1,i))+1/2*(x(2,i)-x(1,i)) ; 1/5*(2-7*cos(pi/12*t(i))-x(2,i))+4*cos(pi/6*t(i))+13/7*(22-x(2,i))-1/2*(x(2,i)-x(1,i))];
%Ahora Euler implícito
K=[(1/4)*(2-7*cos((pi/12)*t(i))-y(1,i))+((10*t(i))/(1+t(i)))+(5/3)*(22-y(1,i))+1/2*(y(2,i)-y(1,i)) ; 1/5*(2-7*cos(pi/12*t(i))-y(2,i))+4*cos(pi/6*t(i))+13/7*(22-y(2,i))-1/2*(y(2,i)-y(1,i))];
y(:,i+1)=y(:,i)+h*[(1/4)*(2-7*cos((pi/12)*t(i))-(y(1,i))+1/2*h*K(1,1))+((10*t(i))/(1+t(i)))+(5/3)*(22-(y(1,i)+1/2*h*K(1,1)))+1/2*((y(2,i)+1/2*h*K(2,1))-(y(1,i)+1/2*h*K(1,1))) ; 1/5*(2-7*cos(pi/12*t(i))-(y(2,i)+1/2*h*K(2,1)))+4*cos(pi/6*t(i))+13/7*(22-(y(2,i)+1/2*h*K(2,1)))-1/2*((y(2,i)+1/2*h*K(2,1))-(y(1,i)+1/2*h*K(1,1)))];
%Ahora Runge-Kutta
K1=[1/4*(2-7*cos((pi/12)*t(i))-z(1,i))+(10*t(i)/(1+t(i)))+5/3*(22-z(1,i))+1/2*(z(2,i)-z(1,i)) ; 1/5*(2-7*cos(pi/12*t(i))-z(2,i))+4*cos(pi/6*t(i))+13/7*(22-z(2,i))-1/2*(z(2,i)-z(1,i))];
K2=[1/4*(2-7*cos((pi/12)*t(i))-(z(1,i)+1/2*h*K1(1)))+(10*t(i)/(1+t(i)))+5/3*(22-(z(1,i)+1/2*h*K1(1)))+1/2*((z(2,i)+1/2*h*K1(2))-(z(1,i)+1/2*h*K1(1))) ; 1/5*(2-7*cos(pi/12*t(i))-(z(2,i)+1/2*h*K1(2)))+4*cos(pi/6*t(i))+13/7*(22-(z(2,i)+1/2*h*K1(2)))-1/2*((z(2,i)+1/2*h*K1(2))-(z(1,i)+1/2*h*K1(1)))];
K3=[1/4*(2-7*cos((pi/12)*t(i))-(z(1,i)+1/2*h*K2(1)))+(10*t(i)/(1+t(i)))+5/3*(22-(z(1,i)+1/2*h*K2(1)))+1/2*((z(2,i)+1/2*h*K2(2))-(z(1,i)+1/2*h*K2(1))) ; 1/5*(2-7*cos(pi/12*t(i))-(z(2,i)+1/2*h*K2(2)))+4*cos(pi/6*t(i))+13/7*(22-(z(2,i)+1/2*h*K2(2)))-1/2*((z(2,i)+1/2*h*K2(2))-(z(1,i)+1/2*h*K2(1)))];
K4=[1/4*(2-7*cos((pi/12)*t(i))-(z(1,i)+h*K3(1)))+(10*t(i)/(1+t(i)))+5/3*(22-(z(1,i)+h*K3(1)))+1/2*((z(2,i)+h*K3(2))-(z(1,i)+h*K3(1))) ; 1/5*(2-7*cos(pi/12*t(i))-(z(2,i)+h*K3(2)))+4*cos(pi/6*t(i))+13/7*(22-(z(2,i)+h*K3(2)))-1/2*((z(2,i)+h*K3(2))-(z(1,i)+h*K3(1)))];
z(:,i+1)=z(:,i)+h/6*(K1+2*K2+2*K3+K4);
end
%Temperatura exterior
M=2-7*cos(pi*t/12);
%Dibujamos las gráficas
figure (1)
hold on
plot(t,M,'c-*')
plot(t,x(1,:),'y')
plot(t,x(2,:),'m')
plot(t,y(1,:),'r')
plot(t,y(2,:),'g')
plot(t,z(1,:),'b')
plot(t,z(2,:),'k')
title('Paso h=0.001')
legend('Tª Exterior','TªA Euler','TªB Euler','TªA Euler implícito','TªB Euler implícito','TªA Runge-Kutta','TªB Runge-Kutta','Location','best')
xlabel('Tiempo en horas')
ylabel('Temperatura A y Temperatura B')
hold off
}}
====Euler, Euler implícito y Runge-Kutta de orden 4 con h=0.1 ====
En la siguiente gráfica se muestran las temperaturas interiores en cada zona calculadas con h=0.1 por los tres métodos citados, así como la temperatura del exterior. Se puede observar que con este tamaño de paso, que es relativamente pequeño las gráficas de cada temperatura interior obtenidas con los distintos métodos difieren entre ellas. Aunque la diferencia no es muy grande si es apreciable por lo que no sería un tamaño de paso demasiado conveniente. 
[[Archivo:H01.jpg|marco|centro]]
La siguiente imagen es un zoom de anterior gráfica para que se pueda apreciar con mayor detalle la diferencia de las soluciones halladas con cada método.
[[Archivo:H01zoom.jpg|marco|centro]]

==== Euler, Euler implícito y Runge-Kutta de orden 4 con h=0.01  ====
En la siguiente gráfica se muestran las temperaturas interiores en cada zona calculadas con h=0.01 por los tres métodos citados, así como la temperatura del exterior. Se puede observar que con este tamaño de paso, que ya comienza a ser más considerable las gráficas de cada temperatura interior obtenidas con los distintos métodos prácticamente no difieren entre ellas. A primera vista y sin realizar un zoom se podría decir que las curvas están superpuestas. Este tamaño ya puede resultarnos lo suficientemente aproximado.
[[Archivo:H001.jpg|marco|centro]]
En la siguiente imagen se puede observar que para apreciar la separación de las curvas el zoom que se debe realizar es mucho mayor que el anterior.
[[Archivo:H001zoom.jpg|marco|centro]]

==== Euler, Euler implícito y Runge-Kutta de orden 4 con h=0.001 ====
En la siguiente gráfica se muestran las temperaturas interiores en cada zona calculadas con h=0.001 por los tres métodos citados, así como la temperatura del exterior. Se puede observar que con este tamaño de paso las curvas se acercan aún más que en el caso anterior, por lo cual este tamaño de paso nos daría un salución aún más aproximada a la realidad.
[[Archivo:H0001.jpg|marco|centro]]
En la imagen se observa que para lograr ver que las curvas no están realmente superpuestas, el zoom que hay que realizar es enorme y por supuesto incluso mayor que el anterior.
[[Archivo:H0001zoom.jpg|marco|centro]]