Modelo Térmico en interior de edificio(Grupo 16-C)

De MateWiki
Revisión del 17:37 10 mar 2015 de Luismorenolopez94 (Discusión | contribuciones) (→‎Primera situación del Modelo)

Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Modelo Térmico en interior de edificio. Grupo 16-C
Asignatura Ecuaciones Diferenciales
Curso Curso 2014-15
Autores Laura Ramos Sangrós

Luis Moreno López

Alejandro Alcocer Jiménez

Álvaro Martínez de Andrés

Gonzalo Olmos de la Cruz

Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


Utilizando ecuaciones diferenciales se puede crear un modelo matemático que permita determinar la temperatura de un edificio en cualquier tiempo dado T(t), considerando que el edificio es una sola habitación. La razón de cambio de la temperatura está dada por todos los factores que generan o disipan el calor.

Es necesario tomar en cuenta tres factores: llamaremos al calor generado por las personas, computadoras, luces y demás artefactos H(t); también hay que considerar el calentamiento o enfriamiento provocado por la calefacción o aire acondicionado. A dicha función la llamaremos U(t) (las podríamos llamar razones respecto al tiempo).

El tercer factor es el efecto de la temperatura exterior sobre el edificio; este factor se puede modelar mediante la ya conocida ley de enfriamiento de Newton, que establece que hay una razón de cambio de la temperatura T(t) que es proporcional a la diferencia entre la temperatura exterior M(t) menos la interior T(t), es decir: 

                                [math]\frac{\partial T}{\partial t}=K[M(t)-T(t)] [/math]

La constante K depende de las propiedades físicas del edificio, es decir, de la cantidad de puertas, ventanas, etc. Pero no depende de M,T o t. Por lo tanto vemos que cuando M es mayor que T, la temperatura del edificio aumenta, y si M es menor que T, la temperatura del edificio disminuye. Resumiendo vemos que: 

                           [math]\frac{\partial T}{\partial t}=K[M(t)-T(t)] + H(t) + U(t) [/math]

1 Primera situación del Modelo

Teniendo el problema introducido,la primera situación que se nos plantea es averiguar cuánto tarda en cambiar considerablemente la temperatura del edificio. Para contestar a esto se nos supone que al final del día (to = 0), la temperatura exterior permanece a Mo grados C,la razón de calentamiento adicional H(t) en el interior del edificio es 0, y la razón de calefacción o aire acondicionado U(t) es también 0. Suponiendo por último que la temperatura en el interior del edificio a medianoche es T(0) = To grados C,la temperatura T(t) en función del tiempo t satisface: 

                           [math]\ T(t)-Mo = (To-Mo)e^{-kt} [/math]

                    [Esta solución la utilizaremos para obtener la curva de la solución exacta y compararla con la del método numérico]

La primera cuestión,que es la anterior citada, la abordamos suponiendo que la temperatura en el interior del edificio a media noche(t=0) es To=14 grados,que la temperatura exterior es constante e igual a 8 grados(Mo = 8) y que la constante del edificio es de 3 horas(k=1/3). Así, utilizando la ley de enfriamiento de Newton [math]\ T'(t)=K[Mo-T(t)] [/math] y con los datos dados,se procede a determinar como evoluciona la temperatura a lo largo del día que empieza mediante el método de Euler Implícito a partir del PVI:: 

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

Codificamos en Matlab para obtener la solución aproximada:: 

%Variacion de la temperatura del edificio con Euler implícito
clear all,clf

t0 = 0;                  %Tiempo incicial(media noche)
tn = 24; 
T0 = 14;                  %Tiempo final(24h del dia)
h = 0.01;                %Longitud de paso
N =(tn-t0)/h;            %Numero de nodos
t = linspace(0,24,N+1);  %Definimos var.indep.(Tiempo)

%Solución exacta
%T(t)=Mo+(To-Mo)*exp(-k*t) 
 Te = (14-8).*exp(-1/3*t)+8;

T = zeros(size(t));      %Creamos un vector de ceros en T
T(1) = T0;               %El valor de T en el primer instante

for i = 1:N
  %Despejando T(i+1) de la ecuacion de Euler Implícito sale
  T(i+1) = 1/(1+h/3)*(T(i)+h*(8/3));
end

%Sacamos tabla de resultados,Temperatura en abcisas y Tiempo en ordenadas
[T',t',Te'];

%Gráfica
hold on
ylabel('Temperatura'); xlabel('Tiempo(h)');
plot(t,Te,'r','LineWidth', 7);
plot(t,T, 'b','LineWidth', 5);
legend('Solucion exacta','Euler', 'location','best')
hold off
CurvaAP1bis.jpg

APRECIACIONES DEL GRÁFICO

Desde el punto de vista numérico se aprecia la solución exacta está superpuesta con la de Euler,lo que quiere decir que Euler es bastante exacto con ese tamaño de paso.

Se observa un descenso de temperatura a lo largo del día muy brusco las primeras horas pero a medida que pasa el día se va suavizando hasta terminar casi constante en las últimas horas.Esto sucede a causa de que no existe ningún calentamiento o enfriamiento adicional que lo modifique,que más tarde veremos que ocurre si varían estas condiciones.

2 Segunda situación del Modelo

En segundo lugar se exige que el calentamiento adicional H(t) es igual a una constante(Ho), pero aún no hay calefacción o enfriamiento (U(t)=0). En este segundo caso, la temperatura exterior M(t) varía en forma de onda senoidal durante un periodo de 24 horas,su mínimo en t=0(medianoche) y su máximo en t=12(mediodía),dato importante que se reflejará en las posteriores gráficas. 

                                                                   [math]\ M(t)=Mo-Bcos(wt) [/math]

donde Mo y B con constantes y [math]\ w=\frac{pi}{12}\ [/math] en radianes/hora y T(0)=To Demostramos que la solución de esta situación es:(dada por el enunciado del trabajo):

[math] T(t)=Bo - BF(t) + Ce^{-kt} [/math] donde\[\left\{\begin{matrix}\ Bo = Mo + \frac{Ho}{k}\ , & \\ F(t)= \frac{\cos(wt)+\frac{w}{k}\sin(wt)}{1+\frac{w^2}{k^2}}\ \, \\ C= To - Bo + \frac{B}{1+\frac{w^2}{k^2}}\\ & \end{matrix}\right.\]

Entonces el problema de valor inicial pedido será el siguiente \[\left\{\begin{matrix}\ T' = k[(Mo-Bcos(wt)-T]+Ho , & \\ T(To) = To & \end{matrix}\right.\]

Que es una EDO lineal de 1º Orden en T cuya expresión ordenada es 

                                                                       [math] T'+KT = kMo+Ho-kBcos(wt) [/math]

Resolviendo la parte homogénea (Variables separables) :: 

                           [math] T'+kT=0 \rightarrow \frac{dT}{T}=-kdt \rightarrow T_{H}=ce^{-kt} [/math]

La solución particular se resolverá por el principio de superposición de soluciones ya que lo hace mas cómodo.Se busca una solución para f1(t) = KMo + Ho y otra para f2(t) = -kBcos(wt).

--En la primera particular se tantea con Tp1 = A,teniendo Tp1'=0; así Tp1'+ kTp1= kMo + Ho y sale:: 

                                          [math] A = Mo + \frac{Ho}{k}\ \rightarrow Tp1 = Mo \frac{Ho}{k}\ [/math]

por la tanto la anterior expresión corresponde a Bo (primera solución particular) y queda demostrado.

--En la segunda solución particular se tantea con Tp2 = R(cos(wt)) + S(sen(wt)) para Tp2' + kTp2 = -kBcos(wt),que derivando Tp2 se obtiene

-Rwsen(wt)+Swcos(wt)+kRcos(wt)+kSsen(wt) = -kBcos(wt), que resolviéndose este sencillo sistema por cramer para sacar R y S: 

                 \[\left\{\begin{matrix}\ Sw + kR = -kB \, & \\  kS - Rw = 0 \ & \end{matrix}\right.\]

Así sustituyendo en la ecuación de tanteo Tp2 con los valores de R y S sale:: 

         [math] Tp2(t)=\frac{-B}{1+{(\frac{w}{k})}^2}cos(wt)-\frac{\frac{w}{k}B}{1+{(\frac{w}{k})}^2}sin(wt)[/math]

Por lo tanto la solución completa seráT(t) = Th + Tp1 + Tp2 

        [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)]}+ce^{-kt}[/math]

Claramente lo que va dentro de los corchetes es F1(t) ya que va multiplicado por B como muestra la solución del enunciado:: 

                           [math] F(t)= \frac{\cos(wt)+\frac{w}{k}\sin(wt)}{1+\frac{w^2}{k^2}}\\ [/math]


Que imponiéndose las condiciones iniciales del valor inicial T(0) = To se obtiene como dice el enunciado:: 

                                     [math] C= To - Bo + \frac{B}{1+\frac{w^2}{k^2}}\\ [/math]

Con esto quedan todas las expresiones que querían que demostrásemos analíticamente.


Lo siguiente será demostrar que la temperatura exterior media Mo coincide con Bo que es aproximadamente la temperatura media diaria dentro del eficicio, si Ho=0. Para ello planteamos lo siguiente:: 

       [math] Tm(t) = T.media = \frac{\int_0^{24}{T(t)dt}}{24}\ = \frac{\int_0^{24}{Bo-BF(t)+Ce^{-kt})dt}}{24}\ ≈ \frac{\int_0^{24}{Bo dt}}{24}\ ≈ Bo [/math]  pues [math] \int_0^{24}{F(t)dt} = 0 [/math]   al estar F(t) integrado por funciones periódicas,y [math]  \int_0^{24}{Ce^{-kt}dt} ≈ 0 [/math]

En el tercer apartado analítico de esta parte se pide demostrar que la variación senoidal dentro del edificio tiene un retraso de α horas respecto de la variación exterior teniendo [math] F(t) = [1+\frac{w^2}{k^2}]^{-1/2} cos(wt-α) ,con, t\gt0 ,donde, tanα=\frac{w}{k}\ [/math]


Si [math] F(t) = [1+\frac{w^2}{k^2}]^{-1/2}(cos(wt) + \frac{w}{k}\sin(wt)) = [1+\frac{w^2}{k^2}]^{-1/2}cos(wt-α) [/math]

Desarrollando el coseno [math] cos(wt) + \frac{w}{k}\sin(wt) = cos(wt)cos(α) + sin(wt)sin(α) [/math]

\[\left\{\begin{matrix}\ 1=cos(α), & \frac{w}{k}=sin(α) & \end{matrix}\right.\] Así [math] Tg(α) = \frac{w}{k}\ [/math]

Así concluimos la explicación analítica de la solución de la ecuación en esta situación del modelo.


Por último en esta parte se pide interpretar la temperatura interior del edificio teniendo en cuenta la exterior,a diferencia del primer apartado, con unos datos que se indican en el enunciado para introducir en el problema de valor inicial y así obtener cómo varía la temperatura.

Así con los valores dadosHo=3;Mo=7;k=1/3;B=5;To=13 el PVI queda:: \[\left\{\begin{matrix}\ T'(t) = 1/3[(7-5cos\frac{Πt}{12}\ -T(t)] + 3 , t>0 , & \\ T(To) = 13 & \end{matrix}\right.\] 

Veámoslo en un gráfico para más tarde interpretar los resultados.

%Apartado 2
clear all,clf
t0 = 0;
tn = 24;
to = 13;
h1 =0.1; h2=h1^2; h3=h1^3;
N1=(tn-t0)/h1; N2=(tn-t0)/h2; N3=(tn-t0)/h3;
%Var.indep.
t = t0:h3:tn;
t1= t0:h1:tn;
t2= t0:h2:tn;
t3= t0:h3:tn;
%Factores de la Solución exacta
B = 5; Bo = 16; C = to-Bo+ 1/(1+(((pi/12)/(1/3))^2))*B; w = pi/12; wk = w*3;
F = 1/(1+wk^2).*(cos(w.*t)+wk.*sin(w.*t));
M=7-B*cos(w*t);
%Solución exacta
Te= Bo-B.*F+C.*exp(-1/(1/3)*t);
%Vectores vacíos donde se guardaran las soluciones aproximadas
y1 = zeros(size(t1)); y2 = zeros(size(t2)); y3 = zeros(size(t3)); %EULER
y1(1) = to;  y2(1) = to;  y3(1) = to;
r1 = zeros(size(t1)); r2 = zeros(size(t2)); r3 = zeros(size(t3)); %KUTTA
r1(1) = to;  r2(1) = to;  r3(1) = to;  
 
for i = 1:N1
  y1(i+1) = y1(i)+h1*(1/3*((7-5*cos(w*t1(i)))-y1(i))+3);    %EULER
  K11 = 1/3*((7-5*cos(w*t1(i)))-r1(i))+3;                   %KUTTA       
  K21 = 1/3*((7-5*cos(w*(t1(i)+h1/2)))-(r1(i)+K11*h1/2))+3;
  K31 = 1/3*((7-5*cos(w*(t1(i)+h1/2)))-(r1(i)+K21*h1/2))+3;
  K41 = 1/3*((7-5*cos(w*(t1(i)+h1)-(r1(i)+K31*h1))))+3;
  r1(i+1) = r1(i)+h1/6*(K11+2*K21+2*K31+K41);
end
for i = 1:N2
  y2(i+1) = y2(i)+h2*(1/3*((7-5*cos(w*t2(i)))-y2(i))+3);    %EULER
  K12 = 1/3*((7-5*cos(w*t2(i)))-r2(i))+3;                   %KUTTA
  K22 = 1/3*((7-5*cos(w*(t2(i)+h2/2)))-(r2(i)+K12*h2/2))+3;
  K32 = 1/3*((7-5*cos(w*(t2(i)+h2/2)))-(r2(i)+K22*h2/2))+3;
  K42 = 1/3*((7-5*cos(w*(t2(i)+h2)-(r2(i)+K32*h2))))+3;
  r2(i+1) = r2(i)+h2/6*(K12+2*K22+2*K32+K42);
end
for i = 1:N3
  y3(i+1) = y3(i)+h3*(1/3*((7-5*cos(w*t3(i)))-y3(i))+3);    %EULER
  K13 = 1/3*((7-5*cos(w*t3(i)))-r3(i))+3;                   %KUTTA
  K23 =1/3*((7-5*cos(w*(t3(i)+h3/2)))-(r3(i)+K13*h3/2))+3;
  K33 = 1/3*((7-5*cos(w*(t3(i)+h3/2)))-(r3(i)+K23*h3/2))+3;
  K43 = 1/3*((7-5*cos(w*(t3(i)+h3))-(r3(i)+K33*h3)))+3;
  r3(i+1) = r3(i)+h3/6*(K13+2*K23+2*K33+K43);
end 
 
%Sacamos tabla de resultados
[t1',y1']; [t2',y2']; [t3',y3'];
[t1',r1']; [t2',r2']; [t3',r3'];
%Gráfico
hold on
xlabel('Tiempo'); ylabel('Temperatura');
plot(t,M,'g','LineWidth',4);
plot(t,Te,'r','LineWidth',  4);     
plot(t1,y1,'g','LineWidth', 4);   
plot(t1,r1,'y','LineWidth', 4);
plot(t2,y2,'c','LineWidth', 4); 
plot(t2,r2,'m','LineWidth', 4);
plot(t3,y3,'k','LineWidth', 4);   
plot(t3,r3,'b','LineWidth', 4);
legend('Temp Exterior','Solución exacta','E 0.1','RK 0.1','E 0.01','RK 0.01','E 0.001','RK 0.001','location','best');
hold off
Def2.jpg

APRECIACIONES IMPORTANTES DEL GRÁFICO::

Se aprecia la curva en forma senoidal por el efecto de la temperatura exterior Mo antes citado.La temperatura crece a lo largo del día hasta casi los 20 grados para bajar unas horas después a la temperatura inicial To.

Una apreciación importante es la temperatura media,antes demostrado que coincide con Bo,es igual a 16 que según la gráfica esto concuerda al unísono,lo que quiere decir que la temperatura media en un día es 16 grados centígrados.

Otra es el retraso en la variación senoidal al comienzo de la curva,donde en las primeras horas baja la temperatura un grado(efecto de la madrugada) hasta que el efecto de la temperatura exterior hace que suba en las horas posteriores(calor del día)llegando a un punto máximo a las 15 horas aproximadamente para luego volver a bajar(efecto del frío de la noche).

Desde el punto de vista numérico en este caso el método de Euler resulta mas eficaz que el de Kutta, al verse que la curva de los pasos de Euler está casi superpuesta en la solución exacta en rojo.Para Kutta la temperatura sube más a lo largo del día.

3 Tercera situación del Modelo

En esta nueva situación se añade la condición de un termostato sencillo que se utiliza para comparar la temperatura interior del edificio con la deseada aportando calor o enfriando dependiendo de si la interior es menor o mayor que la deseada,designada por Td.

Esta cantidad de calor o enfriamiento suministrado se traduce en la siguiente expresión:: 

                                                [math] U(t) = Ku[Td-T(t)] [/math]

Esta expresión condicionará la gráfica antes expuesta,veremos cuál será el resultado si sumamos esta condición con todo lo anterior.

Primeramente nos dicen que la solución analítica del PVI entero con esta última condición es: 

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

con 

  \[\left\{\begin{matrix}\ K1 = k+Ku\ , & \\ B2 = \frac{KuTd+kMo+Ho}{K1}\ \, \\ B1 = \frac{kB}{K1}\ \,  \\ F1(t)= \frac{\cos(wt)+\frac{w}{k1}\sin(wt)}{1+\frac{w^2}{k^2}}\ \, \\ C= To - B2 + F1(0) \ & \end{matrix}\right.\]

Para demostrar esto se procede a resolver analíticamente la ecuación lineal en T siguiente: 

                                 [math] T' = k[(Mo-Bcos(wt)]-T]+Ho+Ku[Td-T] [/math]

que ordenada para ver más claramente que es lineal queda:: 

                              [math] T' + kT + KuT = kMo + Ho - kBcos(wt) + KuTd [/math]

La ecuación homogénea es:: [math] T'+KT + KuT = 0 [/math] que es resuelta por variables separables da:: [math] T = Ce^{-(k+Ku)} [/math] siendo k+kU = K1 ,quedando demostrada la primera variable.

La solución particular la dividimos en dos mediante el principio de superposición como en la situación 2, Tp1 y Tp2 que las distintas funciones serán: f1(t) = kMo+Ho+KuTd y f2(t) = -kBcos(wt) que resolveremos por el método de tanteo de soluciones.

La primera solución vendrá determinada por el siguiente tanteo::

[math] Tp1 = A \rightarrow Tp1' = 0 [/math] sustituyendo en la EDOL [math] Tp1'+kTp1+KuTp1 = kMo+Ho+KuTd \rightarrow 0+Ak+AKu = kMo+Ho+KuTd \rightarrow A = \frac{KuTd+kMo+Ho}{k + Ku}\ [/math] Así la primera solución particular será:: 

                                             [math] Tp1 = \frac{KuTd+kMo+Ho}{K1}\ [/math]

que, claramente, esta expresión coincide con B2 por lo que queda demostrado que B2 es la primera solución particular.

Vamos con la segunda solución particular que tanteo con Tp2 = R cos(wt) + S sen(wt) , R y S reales que derivando Tp2' = -Vsen(wt) + Sw cos(wt). Voy a la EDOL y sustituyo:: 

                     [math] -Rw sen(wt)+Sw cos(wt)+kR cos(wt)+kS sen(wt)+KuR cos(wt)+KuS sen(wt) = -kBcos(wt) [/math]

que para encontrar los valores de R y S se plantea un sistema sencillo:

\[\left\{\begin{matrix}\ Sw + kR + KuR = -kB \, & \\ -Rw + kS + KuS = 0 \ & \end{matrix}\right.\]

Por Cramer sacamos sus valores y se sustituyen en la primera expresión de Tp2 quedando:: 

[math] Tp2(t)=-\frac{\frac{kB}{K_{1}}}{1+{(\frac{w}{K_{1}})}^2}cos(wt)-\frac{\frac{KB}{K_{1}^2}w}{1+{(\frac{w}{K_{1}})}^2}sin(wt) [/math]

Así la solución general estará dada por: T(t) = Th + Tp1 + Tp2::

[math]T(t)=\frac{kMo+Ho+kUT_{D}}{K_{1}}-B1{[\frac{\frac{k}{K_{1}}}{1+{(\frac{w}{K_{1}})}^2}cos(wt)+\frac{\frac{k}{K_{1}^2}w}{1+{(\frac{w}{k_{1}})}^2}sin(wt)]}+ce^{-kt}[/math]

Quedando claramente demostrado que,identificando coeficientes, F1(t) es lo que hay entre los corchetes multiplicando a B1 como se ve en la solución dada por el enunciado.

Solo queda demostrar C que, imponiendo las condiciones iniciales del PVI T(0) = 0 , queda:: 

                                         [math] C = To - Bo + \frac{B}{1+\frac{w^2}{k^2}}\\ [/math]

Terminada la parte analítica de demostración se nos exige deducir que B2 representa aproximadamente la temperatura media del edificio, muy próxima a Td, teniendo en cuenta que cuando el 'calentador' se apaga transcurren 30 minutos para que el término exponencial de la solución T(t) desaparezca; por lo tanto lo despreciamos en los cálculos [math] \int_0^{24}{Ce^{-K1t}} = 0 [/math].

Con esto, procedemos a calcular la temperatura media. 

 [math] Tm(t) = T.media = \frac{\int_0^{24}{T(t)dt}}{24}\ = \frac{\int_0^{24}{B2-B1F1(t)dt}}{24}\ ≈ \frac{\int_0^{24}{B2 dt}}{24}\ ≈ B2 [/math] 

pues       [math] \int_0^{24}{F1(t)dt} = 0  [/math]      al estar F1(t) formada por funciones periódicas(cos y sen).


Así, queda visto que el término B2 representa la temperatura media de la curva y su valor es 19'5 Cº, calculado a partir de la solución exacta codificada en Matlab que más tarde se demostrará.


Por último con los datos proporcionados de las constantes(Ho=3,Mo=4,k=1/3,B=5,To=24,Ku=7/8,Td=22), el PVI que modeliza esta situación es:: [math] \left \{\begin{matrix} T'(t) = 1/3[4-5cos\frac{Πt}{12}\ -T(t)] + 3 + \frac{7(22-T(t)}{8}\ , t\gt0\\ T(0)=24\end{matrix}\right . [/math]

Procedemos a la parte numérica elaborando un código que me permita encontrar la curva resultante con esta nueva condición,calculando la solución exacta y la aproximada por Euler y Runge-Kutta con pasos h=0.1, h=0.01, h=0.001; resultando:: 

%Apartado 3
%Euler y Rungge-Kutta
clear all,clf
t0 = 0;
tn = 24;
to = 24;
h1 =0.1; h2=h1^2; h3=h1^3;
N1=(tn-t0)/h1; N2=(tn-t0)/h2; N3=(tn-t0)/h3;
%Var.indep.
t = t0:h3:tn;
t1= t0:h1:tn;
t2= t0:h2:tn;
t3= t0:h3:tn;

%Factores de la Solución exacta
Td=22; Mo=4; Ho=3; k=1/3; KU=7/8; k1=k+KU; 
B1 = (k*5)/k1; 
B2 = 1/k1*(KU*Td+k*Mo + Ho) 
w = pi/12; 
F = 1/(1+(w^2/k1^2)).*(cos(w.*t)+w/k1.*sin(w.*t));
C = to-B2+1/(1+(w^2/k1^2));
%Solución exacta
Te= B2-B1.*F+C.*exp(-k1*t);

%Vectores vacíos donde se guardaran las soluciones aproximadas
y1 = zeros(size(t1)); y2 = zeros(size(t2)); y3 = zeros(size(t3)); %EULER
y1(1) = to;  y2(1) = to;  y3(1) = to;
r1 = zeros(size(t1)); r2 = zeros(size(t2)); r3 = zeros(size(t3)); %KUTTA
r1(1) = to;  r2(1) = to;  r3(1) = to;  

for i = 1:N1
  y1(i+1) = y1(i)+h1*(1/3*((Mo-5*cos(w*t1(i)))-y1(i))+3+7/8*(22-y1(i)));    %EULER
  K11 = 1/3*(Mo-5*cos(w*t1(i))-r1(i))+3+7/8*(22-y1(i));                   %KUTTA       
  K21 = 1/3*(Mo-5*cos(w*(t1(i)+h1/2))-(r1(i)+K11*h1/2))+3+7/8*(22-(r1(i)+K11*h1/2));
  K31 = 1/3*(Mo-5*cos(w*(t1(i)+h1/2))-(r1(i)+K21*h1/2))+3+7/8*(22-(r1(i)+K21*h1/2));
  K41 = 1/3*(Mo-5*cos(w*(t1(i)+h1))-(r1(i)+K31*h1))+3+7/8*(22-(r1(i)+K31*h1));
  r1(i+1) = r1(i)+h1/6*(K11+2*K21+2*K31+K41);
end
for i = 1:N2
  y2(i+1) = y2(i)+h2*(1/3*((Mo-5*cos(w*t2(i)))-y2(i))+3+7/8*(22-y2(i)));     %EULER
  K12 = 1/3*(Mo-5*cos(w*t2(i))-r2(i))+3+7/8*(22-y2(i));                   %KUTTA
  K22 = 1/3*(Mo-5*cos(w*(t2(i)+h2/2))-(r2(i)+K12*h2/2))+3+7/8*(22-(r2(i)+K12*h2/2));
  K32 = 1/3*(Mo-5*cos(w*(t2(i)+h2/2))-(r2(i)+K22*h2/2))+3+7/8*(22-(r2(i)+K22*h2/2));
  K42 = 1/3*(Mo-5*cos(w*(t2(i)+h2))-(r2(i)+K32*h2))+3+7/8*(22-(r2(i)+K32*h2));
  r2(i+1) = r2(i)+h2/6*(K12+2*K22+2*K32+K42);
end
for i = 1:N3
  y3(i+1) = y3(i)+h3*(1/3*((Mo-5*cos(w*t3(i)))-y3(i))+3+7/8*(22-y3(i)));     %EULER
  K13 = 1/3*(Mo-5*cos(w*(t3(i)))-r3(i))+3+7/8*(22-y3(i));                   %KUTTA
  K23 = 1/3*(Mo-5*cos(w*(t3(i)+h3/2))-(r3(i)+K13*h3/2))+3+7/8*(22-(r3(i)+K13*h3/2));
  K33 = 1/3*(Mo-5*cos(w*(t3(i)+h3/2))-(r3(i)+K23*h3/2))+3+7/8*(22-(r3(i)+K23*h3/2));
  K43 = 1/3*(Mo-5*cos(w*(t3(i)+h3))-(r3(i)+K33*h3))+3+7/8*(22-(r3(i)+K33*h3));
  r3(i+1) = r3(i)+h3/6*(K13+2*K23+2*K33+K43);
end 

%Sacamos tabla de resultados
[t1',y1']; [t2',y2']; [t3',y3'];
[t1',r1']; [t2',r2']; [t3',r3'];
%Gráficasio
hold on
xlabel('Tiempo'); ylabel('Temperatura');
plot(t,Te,'r','LineWidth',  4);     
plot(t1,y1,'g','LineWidth', 4);   
plot(t1,r1,'y','LineWidth', 4);
plot(t2,y2,'c','LineWidth', 4); 
plot(t2,r2,'m','LineWidth', 4);
plot(t3,y3,'k','LineWidth', 4);   
plot(t3,r3,'b','LineWidth', 4);
legend('Solución exacta','E 0.1','RK 0.1','E 0.01','RK 0.01','E 0.001','RK 0.001','location','best');
hold off


CurvaAP324bisbis.jpg

APRECIACIONES IMPORTANTES DEL GRÁFICO:: A primera vista se aprecia una disminución de la temperatura inicial lo cual es totalmente lógico porque al estar a 24 grados(To) el interior del edificio, que es mayor a la temperatura deseada, el termostato empieza a funcionar y a enfriar la habitación a una temperatura estable y confortante alrededor de los 20 grados que por efecto de la temperatura exterior varía de forma senoidal como lo mostramos en la segunda situación.

Efectivamente la temperatura media coincide con B2 que su valor es exactamente 19'5 grados.La temperatura varía hacia arriba y abajo de este valor

Desde el punto de vista analítico todas las curvas de los métodos con sus respectivos tamaños de paso están muy juntas, casi superpuestas, lo que indica que los dos métodos son bastante exactos.

Como curiosidad mostramos en un período de dos días cómo la temperatura se estabiliza y se comporta de la misma forma gracias a la acción del termostato y la temperatura exterior en el gráfico de más abajo. Cambiamos tn = 24 horas por tn = 48 horas.

CuvaAP348bisbis.jpg

4 Cuarta situación del Modelo

En esta última parte del trabajo ya no se añade ningún factor más para ver el comportamiento de la temperatura en una habitación del edificio sino que suponemos que el edificio consta de dos habitaciones A y B que se transfieren calor una a otra.

Las razones de calentamiento adicional estarán dadas por Ha(t),Ua(t) y Hb(t),Ub(t) respectivamente. M(t) que sigue siendo la influencia de temperatura exterior es la misma que antes [math] M(t) = Mo-Bcos\frac{Πt}{12}\ [/math]

Las constantes de tiempo de transferencia de calor son de 4 horas entre A y el exterior,5 horas entre B y el exterior y 2 horas entre A y B,que siendo esta constante: t = 1/k ,despejando k será en cada caso 1/4, 1/5 y 1/2, respectivamente.

La primera cuestión a resolver es determinar el problema de Cauchy que satisface Ta(t) y Tb(t) si cada una designa la temperatura de su zona(A y B) cuando hayan pasado 24 horas,es decir, t ε [0,24] y Ta(0) = 20 y Tb(0) = 18.

Claramente el PVI será un sistema de ecuaciones diferenciales ordinarias de orden uno al estar las ecuaciones conectadas con el paso de temperatura de una a otra. Tomando el problema de Cauchy del comportamiento de la temperatura de cada zona teniendo en cuenta el calentamiento adicional y la transferencia de calor de una a otra, el PVI del SEDOL quedará::

\[\left\{\begin{matrix}\ Ta'(t) = 1/4[M(t) -Ta(t)] + Ha(t) + Ua(t) + \frac{Ta(t)-Tb(t)}{2}\, & \\ Tb'(t) = 1/5[M(t) -Tb(t)] + Hb(t) + Ub(t) - \frac{Ta(t)-Tb(t)}{2}\ \, \\ Ta(0) = 20 \, \\ Tb(0) = 18 & \end{matrix}\right.\]

Como apartado final,teniendo ya el PVI, se nos dan una serie de valores a las constantes para que apliquemos el método de Euler, Euler implícito y Runge-Kutta con los mismos pasos que la situación anterior (h=0.1;0.01;0.001) para que mostremos,con ayuda del Matlab, como se comportan las curvas de ambas temperaturas en un mismo gráfico.

Los valores dados son los siguientes:: \[\left\{\begin{matrix}\ M(t) = 2-7cos\frac{Πt}{12}\ , & \\ Ha(t) = \frac{10t}{1+t}\ , & \\ Ua(t) = \frac{5(22-Ta(t))}{3}\ , & \\ Hb(t) = 4cos\frac{Πt}{6}\ , & \\ Ub(t) = \frac{13(23-Tb(t))}{7}\ & \end{matrix}\right.\]

El PVI definitivo del SEDOL queda:: 

\[\left\{\begin{matrix}\ Ta'(t) = 1/4[2-7cos\frac{Πt}{12}\ -Ta(t)] + \frac{10t}{1+t}\ + \frac{5(22-Ta(t))}{3}\ + \frac{Ta(t)-Tb(t)}{2}\, & \\ Tb'(t) = 1/5[2-7cos\frac{Πt}{12}\ -Tb(t)] + 4cos\frac{Πt}{6}\ + \frac{13(23-Tb(t))}{7}\ - \frac{Ta(t)-Tb(t)}{2}\ \, \\ Ta(0) = 20 \, \\ Tb(0) = 18 & \end{matrix}\right.\]

Codificamos en Matlab para sacar la solución aproximada por estos métodos 

%Apartado 4
%Euler,Euler implícito y Rungge-Kutta
clear all,clf
t0 = 0;
tn = 24;
to = [20,18]';
h1 =0.1; h2=h1^2; h3=h1^3;
N1=round((tn-t0)/h1); N2=round((tn-t0)/h2); N3=round((tn-t0)/h3);
%Var.indep.
%t = t0:h3:tn;
t1= t0:h1:tn;
t2= t0:h2:tn;
t3= t0:h3:tn;

%Vectores vacíos donde se guardaran las soluciones aproximadas
y1 = zeros(2,size(t1)); y2 = zeros(2,size(t2)); y3 = zeros(2,size(t3)); %EULER
y1(:,1) = to;  y2(:,1) = to;  y3(:,1) = to;
m1 = zeros(2,size(t1)); m2 = zeros(2,size(t2)); m3 = zeros(2,size(t3)); %EULER IMPLICITO
m1(:,1) = to;  m2(:,1) = to;  m3(:,1) = to;
r1 = zeros(2,size(t1)); r2 = zeros(2,size(t2)); r3 = zeros(2,size(t3)); %KUTTA
r1(:,1) = to;  r2(:,1) = to;  r3(:,1) = to;  
ka=1/4 ;kb=1/5;
%La matriz A se compone de los coeficientes de Ta y Tb en Ta', situados la primera fila, y de Ta y Tb en Tb' ,que van en la segunda fila
A = [-17/12,-1/2;-1/2,-109/70];
for i =1:N1
  y1(:,i+1) = y1(:,i)+h1.*(A*y1(:,i)+[ka*(2-7*cos(pi*t1(i)/12))+10*t1(i)/(1+t1(i))+110/3;kb*(2-7*cos(pi*t1(i)/12))+4*cos(pi*t1(i)/6)+299/7]);              %EULER
  m1(:,i+1) = (eye(2)-h1*A)\(m1(:,i)+h1*[ka*(2-7*cos(pi*t1(i+1)/12))+10*t1(i+1)/(1+t1(i+1))+110/3;kb*(2-7*cos(pi*t1(i+1)/12))+4*cos(pi*t1(i+1)/6)+299/7]); %EULER IMPLICITO
  K11 = A*r1(:,i)+[ka*(2-7*cos(pi*t1(i)/12))+(10*t1(i)/(1+t1(i)))+110/3;kb*(2-7*cos(pi*t1(i)/12))+(4*cos(pi*t1(i)/6))+299/7];                              %KUTTA       
  K21 = A*(r1(:,i)+h1/2*K11)+[ka*(2-7*cos(pi*(t1(i)+h1/2)/12))+(10*(t1(i)+h1/2)/(1+(t1(i)+h1/2)))+110/3;kb*(2-7*cos(pi*(t1(i)+h1/2)/12))+(4*cos(pi*(t1(i)+h1/2)/6))+299/7];
  K31 = A*(r1(:,i)+h1/2*K21)+[ka*(2-7*cos(pi*(t1(i)+h1/2)/12))+(10*(t1(i)+h1/2)/(1+(t1(i)+h1/2)))+110/3;kb*(2-7*cos(pi*(t1(i)+h1/2)/12))+(4*cos(pi*(t1(i)+h1/2)/6))+299/7];
  K41 = A*(r1(:,i)+h1/2*K31)+[ka*(2-7*cos(pi*(t1(i)+h1)/12))+(10*(t1(i)+h1)/(1+(t1(i)+h1)))+110/3;kb*(2-7*cos(pi*(t1(i)+h1)/12))+(4*cos(pi*(t1(i)+h1)/6))+299/7];
  r1(:,i+1) = r1(:,i)+h1/6*(K11+2*K21+2*K31+K41);
end
for i =1:N2
  y2(:,i+1) = y2(:,i)+h2.*(A*y2(:,i)+[ka*(2-7*cos(pi*t2(i)/12))+(10*t2(i)/(1+t2(i)))+110/3,kb*(2-7*cos(pi*t2(i)/12))+(4*cos(pi*t2(i)/6))+299/7]');         %EULER
  m2(:,i+1) = (eye(2)-h2*A)\(m2(:,i)+h2*[ka*(2-7*cos(pi*t2(i+1)/12))+10*t2(i+1)/(1+t2(i+1))+110/3;kb*(2-7*cos(pi*t2(i+1)/12))+4*cos(pi*t2(i+1)/6)+299/7]); %EULER IMPLICITO
  K12 = A*r2(:,i)+[ka*(2-7*cos(pi*t2(i)/12))+(10*t2(i)/(1+t2(i)))+110/3;kb*(2-7*cos(pi*t2(i)/12))+(4*cos(pi*t2(i)/6))+299/7];                              %KUTTA       
  K22 = A*(r2(:,i)+h2/2*K12)+[ka*(2-7*cos(pi*(t2(i)+h2/2)/12))+(10*(t2(i)+h2/2)/(1+(t2(i)+h2/2)))+110/3;kb*(2-7*cos(pi*(t2(i)+h2/2)/12))+(4*cos(pi*(t2(i)+h2/2)/6))+299/7];
  K32 = A*(r2(:,i)+h2/2*K22)+[ka*(2-7*cos(pi*(t2(i)+h2/2)/12))+(10*(t2(i)+h2/2)/(1+(t2(i)+h2/2)))+110/3;kb*(2-7*cos(pi*(t2(i)+h2/2)/12))+(4*cos(pi*(t2(i)+h2/2)/6))+299/7];
  K42 = A*(r2(:,i)+h2/2*K32)+[ka*(2-7*cos(pi*(t2(i)+h2)/12))+(10*(t2(i)+h2)/(1+(t2(i)+h2)))+110/3;kb*(2-7*cos(pi*(t2(i)+h2)/12))+(4*cos(pi*(t2(i)+h2)/6))+299/7];
  r2(:,i+1) = r2(:,i)+h2/6*(K12+2*K22+2*K32+K42);
end
for i =1:N3
  y3(:,i+1) = y3(:,i)+h3.*(A*y3(:,i)+[ka*(2-7*cos(pi*t3(i)/12))+(10*t3(i)/(1+t3(i)))+110/3,kb*(2-7*cos(pi*t3(i)/12))+(4*cos(pi*t3(i)/6))+299/7]');            %EULER
  m3(:,i+1) = (eye(2)-h3*A)\(m3(:,i)+h3*[ka*(2-7*cos(pi*t3(i+1)/12))+10*t3(i+1)/(1+t3(i+1))+110/3;kb*(2-7*cos(pi*t3(i+1)/12))+4*cos(pi*t3(i+1)/6)+299/7]); EULER IMPLICITO
  K13 = A*r3(:,i)+[ka*(2-7*cos(pi*t3(i)/12))+(10*t3(i)/(1+t3(i)))+110/3;kb*(2-7*cos(pi*t3(i)/12))+(4*cos(pi*t3(i)/6))+299/7];                                 %KUTTA       
  K23 = A*(r3(:,i)+h3/2*K13)+[ka*(2-7*cos(pi*(t3(i)+h3/2)/12))+(10*(t3(i)+h3/2)/(1+(t3(i)+h3/2)))+110/3;kb*(2-7*cos(pi*(t3(i)+h3/2)/12))+(4*cos(pi*(t3(i)+h3/2)/6))+299/7];
  K33 = A*(r3(:,i)+h3/2*K23)+[ka*(2-7*cos(pi*(t3(i)+h3/2)/12))+(10*(t3(i)+h3/2)/(1+(t3(i)+h3/2)))+110/3;kb*(2-7*cos(pi*(t3(i)+h3/2)/12))+(4*cos(pi*(t3(i)+h3/2)/6))+299/7];
  K43 = A*(r3(:,i)+h3/2*K41)+[ka*(2-7*cos(pi*(t3(i)+h3)/12))+(10*(t3(i)+h3)/(1+(t3(i)+h3)))+110/3;kb*(2-7*cos(pi*(t3(i)+h3)/12))+(4*cos(pi*(t3(i)+h3)/6))+299/7];
  r3(:,i+1) = r3(:,i)+h3/6*(K13+2*K23+2*K33+K43);
end 

%Sacamos tabla de resultados
[t1',y1']; [t2',y2']; [t3',y3'];
[t1',r1']; [t2',r2']; [t3',r3'];
%Gráficas
hold on
xlabel('Tiempo'); ylabel('Temperatura');     
plot(t1,y1,'g','LineWidth', 4);  
plot(t1,m1,'LineWidth', 4);
plot(t1,r1,'y','LineWidth', 4);
plot(t2,y2,'c','LineWidth', 4); 
plot(t2,m2,'LineWidth', 4);
plot(t2,r2,'m','LineWidth', 4);
plot(t3,y3,'k','LineWidth', 4);
plot(t3,m3,'LineWidth', 4); 
plot(t3,r3,'b','LineWidth', 4);
legend('E 0.1','EIm 0.1','RK 0.1','E 0.01','EIm 0.01','RK 0.01','E 0.001','EIm 0.001','RK 0.001','location','best');
hold off


Ap4definitivo.jpg

APRECIACIONES IMPORTANTES DEL GRÁFICO:: Dada la diferencia de temperaturas iniciales en cada habitación (22 en A y 18 en B), se ve que al comienzo de la gráfica la curva de B aumenta rápidamente de temperatura porque se le cede calor desde A ,al plantearse el modelo con ese sentido, mientras que la temperatura de A desciende porque está cediendo calor a B. Esta cedencia de calor tiene un máximo en cada curva pasadas unas 2 horas, así las gráficas se invierten porque no han alcanzado una temperatura estable y esta vez B cede calor a A para regularse de nuevo al intervenir muchos factores en el calentamiento es difícil que se igualen. Pasadas unas 10 horas parece que la curva se estabiliza y A queda siempre con una temperatura mayor que B aunque se siguen cediendo calor en forma senoidal.

La idea principal es que cuando una cede calor otra lo absorbe y viceversa, así sucesivamente.

Desde el punto de vista numérico los métodos dan soluciones muy parecidas, si ampliamos el gráfico se ven pequeñas diferencias siendo el Runge-Kutta teóricamente el más efectivo. Aunque se aprecia una ligera desviación de Euler explícito con paso 0.1 al principio de la curva al ser un método menos preciso y más con ese paso.

Adjuntamos una segunda gráfica pero esta será en un periodo de 48 para observar como se estabiliza esta gráfica en un intervalo más amplio.

CurvaAP4definitiva48.jpg

Con esto damos fin al trabajo, esperamos que les guste, nos ha servido para reforzar conocimientos y aprender cosas interesantes. También que sirva como buen ejemplo para futuros alumnos e investigaciones sobre este tema.

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Modelo_Térmico_en_interior_de_edificio(Grupo_16-C)&oldid=29140»