Revisión actual del 04:01 12 may 2016

El objetivo del trabajo es formular un modelo matemático que describa el comportamiento de la temperatura en el interior de un edificio en el lapso de 24h en función de la temperatura exterior, del calor que se genera dentro del edificio y del sistema de calefacción y aire acondicionado. Se trata de responder a estas tres preguntas: ¿Cuánto tarda en cambiar considerablemente la temperatura del edificio? ¿Cómo varía la temperatura del edificio durante la primavera y el otoño, cuando no se emplea el aire acondicionado o calefacción? ¿Cómo varia la temperatura del edificio durante el verano, cuando se utiliza aire acondicionado, o en invierno, cuando se emplea calefacción?

Para modelizar este suceso utilizaremos la teoría de ecuaciones diferenciales y modelos numéricos resueltos mediante el programa MatLab.

1 Introducción y modelización

Dentro de un edificio existen multitud de parámetros que regulan la temperatura interior, por ejemplo, dependiendo del tipo de paredes, cubierta y forjados que lo compongan, podrá tener una mayor transmitancia térmica y quedar más expuesto a la temperatura exterior. Además hay que tener en cuenta que el calor se puede transmitir por conducción, convección y radiación.

El número, tamaño y calidad de ventanas también influye, ya que representan los puentes térmicos que hacen más o menos vulnerable al edificio.

Planta de un Edificio

Muros tipo de un edificio

Detalle de cubierta ajardinada

Detalle de cubierta no transitable

Las personas y las máquinas también influyen en la temperatura interior, ya que generan calor que transmiten al interior del edificio, de tal manera que según el manual de instalaciones de calefacción por agua caliente de Franco Martín Sánchez alrededor de tres personas generan el mismo calor que un radiador. Evidentemente, los calefactores, radiadores y máquinas de aire acondicionado varían la temperatura interior y son usados para mantener una temperatura ideal que permita a las personas hacer las actividades para las que está diseñado el edificio.

De tal manera que para simplificar todas estas variables se tomará el edificio como un bloque donde todos sus paramentos tendrán la misma capacidad de transmisión de calor (k), el calor producido por las personas, máquinas y calefacción no será definido por zonas puntuales de fuentes o sumideros, sino que serán funciones que tendrán la misma densidad en todos los puntos del edificio; quedando un modelo así:

centro

Funciones y parámetros que participan en la temperatura:

• [math]T(t)=[/math] Temperatura interior del edificio en un instante [math]t[/math]. [Medido en ºC].
• [math]H(t)=[/math] Calor producido por personas y elementos del edificio. "(Expresados en términos de energía por unidad de tiempo que multiplicado por la capacidad calorífica del edificio da temperaturas por unidad de tiempo)"
• [math]U(t)=[/math] Calentamiento producido por la calefacción. "(Expresados en términos de energía por unidad de tiempo que multiplicado por la capacidad calorífica del edificio da temperaturas por unidad de tiempo)"
• [math]M(t)=[/math] Temperatura exterior. "(Medido en ºC)"
• [math]k=[/math] Constante que depende de las propiedades del edificio. "(Expresado en unidades de grados de cambio de temperatura por energía calorífica)"
• [math]t_0=\frac{1}{k}=[/math] Constante del tiempo del edificio, correspondiente con el instante inicial. "(Expresado en horas)"

En nuestro caso particular estudiaremos la evolución de la temperatura en un lapso de 24h, para ello definiremos la variable [math]T(t)[/math] que representa la temperatura en el interior del edificio en un instante [math]t[/math]. Como se puede observar [math]T[/math] será una función dependiente únicamente del tiempo ya que consideramos el interior del edificio como un único habitáculo en el cuál la temperatura es igual para todo el espacio.

Para deducir la ecuación diferencial que gobierna el comportamiento de la temperatura, aceptaremos como válidas dos hipótesis:

• Aceptaremos como válida y extrapolable a nuestro problema la ley de transmisión del calor de Fourier, que afirma que la conducción de calor entre dos cuerpos es proporcional a la diferencia de temperatura de ambos para una cierta constante [math]k[/math] de proporcionalidad.

[math]\frac{dT_1}{dt}=k(T_1-T_2)[/math]

• Consideraremos la temperatura [math]T(t)[/math] igual para todo el interior del edificio.

Con estas hipótesis como premisas deduciremos la ecuación diferencial que rige la temperatura [math]T(t)[/math] en base a que la variación de grados centígrados es igual a los grados centígrados que aumentan menos los que disminuyen. El problema de Cauchy que modela el fenómeno es, por tanto:

\begin{equation*} \begin{cases} T'(t)=\Delta(ºC)=k·(M(t)-T(t))+H(t)±U(t);\hspace{1cm} t≥0 \\ T(0)=T_0 \end{cases} \end{equation*}

2 Tiempo de variación significativa de [math]T(t)[/math]

La primera pregunta que se nos plantea es cuánto tarda en variar considerablemente la temperatura del edificio. Para ello, vamos a resolver numéricamente la ecuación diferencial planteada anteriormente utilizando el método de Euler implícito con un paso h=0.01. La dificultad de usar un método de este tipo reside en que la variable que queremos calcular no aparece de manera explícita en la ecuación, por lo que hay que despejarla manualmente para introducirla en el algoritmo.

Para despejarla, tendremos en cuenta que al final del día (que coincide con el inicio del día, al ser períodos de 24h) la temperatura exterior permanece a [math]M_0[/math] grados; que la razón de calentamiento adicional [math](H_0)[/math] es 0 grados, al igual que la razón de calefacción [math](U_0)[/math]; y que la temperatura en el interior del edificio a medianoche es [math]t_0[/math] grados.

Es decir, la ecuación diferencial quedaría como:

\begin{equation*} \begin{cases} T'(t)=\Delta(ºC)=k(M(t)-T(t))+0+0 \\ T(0)=T_0 \end{cases} \end{equation*}

Ecuación diferencial definitiva:

\begin{equation*} \begin{cases} T'(t)=\Delta(ºC)=k(M(t)-T(t))\\ T(0)=T_0 \end{cases} \end{equation*}

Que tiene como solución de la ecuación la función:

[math]T(t)=(T_0-M_0)·e^{-k·t}+M_0[/math]

En este apartado hemos tomado una temperatura exterior [math]M(t)[/math] constante e igual a 8°C, una constante térmica del edificio igual a 3 y una temperatura inicial a las 00:00h de 14°C. Usando Euler implícito como indica el enunciado, quedaría un bucle:

\begin{equation*} \begin{cases} T_{n+1}=T_n+h·f(t_{n+1},T{n+1})\\ T(0)=T_0 \end{cases} \end{equation*}

[math]T_{n+1}=T_n+h·f(t_{n+1},T{n+1}) = T_n+h·f·[k·(M_0-T_1)] \Longrightarrow[/math]

[math]T_{n+1}=T_n+h·[\frac{1}{3}·(8-T_{n+1})]=T_n+h·[\frac{8}{3}-\frac{1}{3} T_{n+1}]=T_n+\frac{8}{3}h-\frac{1}{3}h·T_{n+1}\Longrightarrow[/math]

[math]T_{n+1}=T_n+\frac{8}{3}h-\frac{1}{3}h·T_{n+1} \Longrightarrow T_{n+1}+\frac{1}{3}h·T_{n+1}=T_n+\frac{8}{3}h\Longrightarrow[/math]

[math]T_{n+1}·(1+\frac{1}{3}h)=T_n+\frac{8}{3}h⟹\boxed{(T_{n+1}=\frac{(T_n+\frac{8}{3}h)}{(1+\frac{1}{3}h)})}[/math]

A continuación, resolvemos el problema numérico en Matlab con los datos anteriormente expuestos:

Solución Exacta comparada con la obtenida por el método de Euler implícito

El error máximo obtenido con el método numérico es [math]0.0036736935ºC[/math]

% CAMBIO CONSIDERABLE DE LA TEMPERATURA DEL EDIFICIO
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all, clc, close all
% Ecuación diferencial: T'(t)=1/3*(M0-T(t));
% Solución exacta T(t)=(T0-M0)*exp(-1/3*t)+M0;

% DATOS DE PARTIDA:
t0=0; tN=24; % Intervalo de tiempo.
k=1/3;       % Cte de tiempo del edificio: t=1/k=3
M0=input('Valor inicial de la temperatura exterior: ');
T0=input('Valor inicial de la temperatura interior: ');

% DISCRETIZACIÓN DE LA VARIABLE INDEPENDIENTE "t":
h=input('Introduzca el tamaño de paso h: '); 
N=round((tN-t0)/h);    % Número de subintervalos.
t=linspace(t0,tN,N+1); % Vector de tiempo.
%--------------------------------------------------------------------------
% SOLUCIÓN EXACTA:
Tex=(T0-M0)*exp(-1/3*t)+M0; % Función solución exacta (solución analítica).

% CÁLCULO NUMÉRICO: Método de Euler Implícito (Backward Euler Method).
T=zeros(size(t)); % Matriz "vacía" de la solución numérica.
T(1)=T0;
for i=1:N         % Inicio del bucle.
  T(i+1)=(T(i)+8*h/3)/(1+h/3);  
end
% ¡Cuidado! Si no hay valor de "T" ó "t" en la EDO, escribir "0*T" ó "0*t".
%--------------------------------------------------------------------------
% GRÁFICAS:
hold on
subplot(1,3,1)  % Primera gráfica.
plot(t,Tex,'r') % Gráfica de la solución exacta (solución analítica).
xlabel('Time'); ylabel('Temperature') 
title('Temperature variation')        
legend('Exact Ecuation')              
subplot(1,3,2)  % Segunda gráfica.
plot(t,T,'b')   % Gráfica de la solución aproximada (solución númerica).
xlabel('Time'), ylabel('Temperature') 
title('Temperature variation')         
legend('Backward Euler Method')       
subplot(1,3,3)          % Tercera gráfica.
plot(t,Tex,'r',t,T,'b') % Gráfica de comparación de métodos.
xlabel('Time'); ylabel('Temperature') 
title('Comparation between methods')        
legend('Exact Ecuation','Backward Euler Method','Location','best')
hold off
%--------------------------------------------------------------------------
% ESTUDIO DE LOS ERRORES:  
format long;       % Todos los decimales.
err=abs(Tex-T);    % Cálculo del error en valor absoluto.
maxi= max(err);    % Máximo error obtenido.
fprintf('El error máximo es %7.10f\n',maxi)
pos=1:1:length(t); % Vector posición (mismo tamaño que el vector "t").
% Tabla recopilación de soluciones y errores:
disp('No=0, Sí=1')
quieres=input('Quieres la tabla? ')
if quieres==1
  Tabl=[pos', t',T',Tex',err'] 
else
  Tabl=[pos', t',T',Tex',err'];
end

3 Variación de la temperatura interior [math]T(t)[/math] sin calefacción ni aire acondicionado [math](U(t)=0)[/math]

En este apartado, el calor producido por personas y elementos del edificio [math]H(t)[/math] se mantiene constante siendo igual a [math]H_0[/math]. No hay calefacción, por lo que [math]U(t)[/math] es igual a cero y la temperatura exterior varía de forma senoidal con un período de 24 horas de la manera: [math]M(t)=M_0-B·cos⁡(ϖ·t)[/math]

Es decir:

1. [math]T(t)=[/math] Temperatura interior del edificio.
2. [math]H(t)=[/math] Calor producido por personas y elementos del edificio[math]=H_0[/math]
3. [math]U(t)=[/math] Calentamiento producido por la calefacción[math]=0[/math]
4. [math]M(t)=[/math] Temperatura exterior[math]=M_0-B·cos⁡(ϖ·t)=M_0-B·cos⁡(\frac{π}{12·t})[/math]
5. [math]k=[/math] Constante que depende de las propiedades del edificio.
6. [math]t_0=\frac{1}{k}=[/math] Constante del tiempo del edificio.

Es decir, la ecuación diferencial quedaría como:

\begin{equation*} \begin{cases} T'(t)=k(M(t)-T(t))+H(t)+U(t)=k(M_0-B·cos⁡(\frac{π}{12·t})-T(t))+H_0+0\\ T(0)=T_0 \end{cases} \end{equation*}

Ecuación diferencial definitiva:

\begin{equation*} \begin{cases} T'(t)=k(M_0-B·cos⁡(\frac{π}{12·t})-T(t))+H_0\\ T(0)=T_0 \end{cases} \end{equation*}

Que tiene como solución:

\begin{equation*} T(t)=B_0-B·F(t)+C·e^{-k·t} \begin{cases} B_0=M_0+\frac{H_0}{k}\\ F(t)=\frac{cos(⁡(ϖ·t)+\frac{ϖ}{k·sen(ϖ·t)}}{(1+(\frac{ϖ}{k})^2)}\\ C=T_0-B_0+\frac{B}{(1+(\frac{ϖ}{k}^2)} \end{cases} \end{equation*}

3.1 Demostración del modelo

Para resolver esta ecuación diferencial usamos el principio de superposición, para poder usarlo, la ecuación diferencial debe ser lineal y los términos independientes deben satisfacer la misma ecuación diferencial homogénea.

[math]T'(t)=k(M_0-B·cos⁡(ϖ·t)-T(t))+H_0=k·M_0-k·B·cos⁡(ϖ·t)-k·T(t)+H_0⟹[/math]

[math]\boxed{T'(t)+k·T(t)=k·M_0-k·B·cos⁡(ϖ·t)+H_0}[/math]

Dividimos la ecuación en dos, una con las constantes y otra con la función coseno:

\begin{equation*} \begin{cases} T'(t)+k·T(t)=k·M_0+H_0\\ T'(t)+k·T(t)=-k·B·cos⁡(ϖ·t) \end{cases} \end{equation*}

1. Cálculo de la solución general de la homogénea:

[math]T'(t)+k·T(t)=0[/math]

Usamos el polinomio característico:

[math]m+k=0 ⇔ m=-k ⇒ S.G.H=c_1·e^{-k·t}[/math]

2. Cálculo de la solución particular de la ecuación con las constantes:: [math]T'(t)+k·T(t)=k·M_0+H_0[/math] Conjunto asociado a la S.G.H:: [math]S.G.H=c_1·e^{-k·t}[/math]

[math]S.G.H≡\{e^{-k·t}\}[/math]

Conjunto asociado al término independiente:: [math]p_1(t)=k·M_0+H_0[/math]

[math]p_1(t)≡\{1\} \nsubseteq S.G.H[/math]

Comparando ambos conjuntos, podemos decir que la solución particular de la completa es del tipo:

[math]T_p=A= \begin{equation*} \begin{cases} T_p=A\\ T'_p=0 \end{cases} \end{equation*} [/math]

Sustituimos en nuestra ecuación diferencial:

[math]T'(t)+k·T(t)=k·M_0+H_0⟺0+k·A=k·M_0+H_0⇒[/math]

[math]k·A=k·M_0+H_0⇒A=M_0+\frac{H_0}{k}[/math]

Es decir, nos quedaría una solución:

[math]T(t)=c_1·e^{-k·t}+M_0+\frac{H_0}{k}[/math]

3. Cálculo de la solución particular de la ecuación con el coseno:

[math]T'(t)+k·T(t)=-k·B·cos⁡(ϖ·t)[/math]

Conjunto asociado a la S.G.H: :[math]S.G.H=c_1·e^{-k·t}[/math]

[math]S.G.H≡\{e^{-k·t}\}[/math]

Conjunto asociado al término independiente: :[math]p_2(t)=-k·B·cos⁡(ϖ·t)[/math]

[math]p_2 (t)≡\{cos⁡(ϖ·t),sen(ϖ·t)\} \nsubseteq S.G.H[/math]

Comparando ambos conjuntos, podemos decir que la solución particular de la completa es del tipo:

[math]T_p=C·cos⁡(ϖ·t)+D·sen(ϖ·t)⇒ \begin{equation*} \begin{cases} T_p=C·cos⁡(ϖ·t)+D·sen(ϖ·t)\\ T'_p=-C·ϖ·sen(ϖ·t)+D·ϖ·cos⁡(ϖ·t) \end{cases} \end{equation*} [/math]

Sustituimos en nuestra ecuación diferencial:

[math]T'(t)+k·T(t)=-k·B·cos⁡(ϖ·t)⇒[/math]

[math]-C·ϖ·sen(ϖ·t)+D·ϖ·cos⁡(ϖ·t)+k·[C·cos⁡(ϖ·t)+D·sen(ϖ·t)]=-k·B·cos⁡(ϖ·t)⇒[/math]

[math][-ϖ·C+D·k]sen(ϖ·t)+[ϖ·D+C·k] cos⁡(ϖ·t)=-k·B·cos⁡(ϖ·t)[/math]

Subdividimos según las funciones:

• [math]sen(ϖ·t): -ϖ·C+B·k=0[/math]

  [math]cos⁡(ϖ·t): ϖ·B+A·k=-k·B[/math]

Y nos queda el sistema de ecuaciones siguiente:

[math] \begin{equation*} \begin{cases} -ϖ·C+D·k=0\\ ϖ·D+C·k=-k·B \end{cases} ⇒ \begin{pmatrix} -ϖ & k \\ k & ϖ \\ \end{pmatrix} · \begin{pmatrix} C \\ D \\ \end{pmatrix} = \begin{pmatrix} 0\\ -k·B\\ \end{pmatrix} \end{equation*} [/math]

Que resolvemos por el método de Cramer:

[math]C=\frac{det⁡(f_3,f_2)}{det⁡(f_1,f_2)}=\frac{k^2·B}{-(ϖ^2+k^2 )}=-\frac{k^2·B}{ϖ^2+k^2}[/math]

[math]D=\frac{det⁡(f_1,f_3)}{det⁡(f_1,f_2)}=\frac{k·ϖ·B}{-(ϖ^2+k^2 )}=-\frac{k·ϖ·B}{ϖ^2+k^2}[/math]

Es decir, nos quedaría una solución:

[math]T(t)=c_1·e^{-k·t}-\frac{k^2·B}{ϖ^2+k^2}·cos⁡(ϖ·t)-\frac{k·ϖ·B}{ϖ^2+k^2}·sen(ϖ·t)=[/math]

[math]c_1·e^{-k·t}-\frac{k^2·B·cos⁡(ϖ·t)+k·ϖ·B·sen(ϖ·t)}{ϖ^2+k^2}=[/math]

[math]c_1·e^{-k·t}-B·\frac{cos⁡(ϖ·t)+\frac{ϖ}{k}·sen(ϖ·t)}{\frac{ϖ}{k})^2+1}[/math]

Como se ha comentado antes, hemos aplicado el principio de superposición, por lo que la solución completa de nuestra ecuación queda como:

[math]T(t)=M_0+\frac{H_0}{k}-B·\frac{cos⁡(ϖ·t)+\frac{ϖ}{k}·sen(ϖ·t)}{(\frac{ϖ}{k})^2+1}+c_1·e^{-k·t}[/math]

La condición inicial nos dice que T(0)=T_0, por lo que:

[math]T(0)=M_0+\frac{H_0}{k}-B·\frac{(cos⁡(ϖ·0)+\frac{ϖ}{k}·sen(ϖ·0)}{(\frac{ϖ}{k})^2+1)}+c_1·e^{-k·0}=T_0⟹[/math]

[math]M_0+\frac{H_0}{k}·\frac{1+\frac{ϖ}{k}·0}{(\frac{ϖ}{k})^2+1}+c_1=T_0⟹\boxed{c_1=T_0+\frac{B}{(\frac{ϖ}{k})^2+1}-(M_0+\frac{H_0}{k})}[/math]

Por lo que finalmente la Función solución queda como:

[math]\boxed{T(t)=M_0+\frac{H_0}{k}-B·\frac{cos⁡(ϖ·t)+\frac{ϖ}{k}·sen(ϖ·t)}{(\frac{ϖ}{k})^2+1}+(T_0+\frac{B}{(\frac{ϖ}{k})^2+1}-(M_0+\frac{H_0}{k}))·e^{-k·t}}[/math]

Y comparándola con el enunciado:

\begin{equation*} T(t)=B_0-B·F(t)+C·e^{-k·t} \begin{cases} B_0=M_0+\frac{H_0}{k}\\ F(t)=\frac{cos(⁡(ϖ·t)+\frac{ϖ}{k·sen(ϖ·t)}}{(1+(\frac{ϖ}{k})^2)}\\ C=T_0-B_0+\frac{B}{(1+(\frac{ϖ}{k}^2)} \end{cases} \end{equation*}

Nos queda de la misma forma.

3.2 Temperatura media diaria exterior e interior. Influencia de [math]H_0[/math]

Siendo [math]B_0[/math] aproximadamente la temperatura media diaria dentro del edificio, se pide comprobar que si [math]H_0=0[/math], entonces [math]B_0[/math] coincide con la temperatura media exterior.

[math]T_{med}=\frac{1}{nº horas}·\sum_{i=1}^nT=\frac{1}{24}·\int_0^{24}T(t)·dt=\frac{1}{24}·\int_0^{24}B_0-B·F(t)+C·e^{-k·t}·dt[/math]

Para que sea más fácil de calcular y seguir, calcularemos la integral separando cada una de los términos:

[math]\frac{1}{24}·\int_0^{24}B_0-B·F(t)+C·e^{-k·t}·dt=\frac{1}{24}·(\int_0^{24}B_0·dt+\int_0^{24}-B·F(t)·dt+\int_0^{24}C·e^{-k·t}·dt)[/math]

• Término [math]B_0[/math]:

[math]\int_0^{24}B_0·dt=B_0\int_0^{24}dt=B_0[t]_0^{24}=24·B_0=24·[M_0+\frac{H_0}{k}][/math]

[math]\boxed{\int_0^{24}B_0·dt=24·[M_0+\frac{H_0}{k}]}[/math]

• Función [math]F(t)[/math]:

[math]\int_0^{24}-B·F(t)·dt=-\frac{B}{1+(\frac{ω}{k})^2}\int_0^{24}cos(ω·t)+\frac{ω}{k}sen(ω·t)·dt[/math]

[math]={k=-\frac{B}{(1+(\frac{ω}{k})^2}}[/math]

[math]=k\biggl[\int_0^{24}cos(ω·t)+\int_0^{24}\frac{ω}{k}sen(ω·t)·dt\biggr][/math]

[math]=k\biggl[\frac{1}{ϖ}\int_0^{24}ω·cos(ω·t)-\frac{1}{k}\int_0^{24}\frac{ω}{k}-ω·sen(ω·t)·dt\biggr][/math]

[math]=k·(\frac{1}{ϖ} [sen(ϖ·t)]_0^{24}+\frac{1}{k} [cos⁡(ϖ·t) ]_{0^24})=\{ϖ=\frac{π}{12}\}=[/math]

[math]=k·(\frac{1}{ϖ} (sen(2·π)-sen(0))+\frac{1}{k} (cos(2·π)-cos(0)))[/math]

[math]=k·(\frac{1}{ϖ} (0-0)+\frac{1}{k} (1-1))=0[/math]

[math]\boxed{\int_0^{24}-B·F(t)·dt=0}[/math]

• Función [math]C·e^{-k·t}[/math]:

[math]\int_0^{24}C·e^{-k·t}·dt=C·\int_0^{24}e^{-k·t}·dt=C·(-\frac{1}{k})\int_0^{24}-k·e^{-k·t}·dt=-\frac{C}{k}·[e^{-k·t}]_0^{24}[/math]

[math]-\frac{C}{k}·(e^{-24·k}-e^0)=-\frac{C}{k}·(e^{-24·k}-1)[/math]

Nos queda una integral completa:

[math]T_{med}=\frac{1}{nº horas}·\sum_{i=1}^nT=\frac{1}{24}·\int_0^{24}T(t)·dt[/math]

[math]=\frac{1}/{24}·[24·[M_0+\frac{H_0}{k}]--\frac{C}{k}·(e^{-24·t}-1)][/math]

[math]=1/24·[24·[M_0+\frac{H_0}{k}]--\frac{1}{k}·(T_0-B_0+\frac{B}{(1+(\frac{ϖ}{k}^2)}·(e^{-24·k}-1)][/math]

Que operando:

[math]=[M_0+\frac{H_0}{k}]+\frac{1}{24k}·\biggl(T_0-B_0+\frac{B}{(1+(\frac{ϖ}{k})^2}\biggr)-\frac{1}{24·k·e^{24·k}}\biggl(T_0-B_0+\frac{B}{(1+(\frac{ϖ}{k}^2)}\biggr)[/math]

El enunciado nos indica que [math]t_0=1/k ∈(2,4)[/math], es decir, que [math]k∈(1/2,1/4)[/math]

• Teniendo en cuenta que [math]k=1/2[/math]

[math][M_0+\frac{H_0}{k}]+\frac{1}{24k}·\biggl(T_0-B_0+\frac{B}{1+(\frac{ϖ}{k})^2}-\frac{1}{24·k·e^{24·k}}·\biggl(T_0-B_0+\frac{B}{1+(\frac{ϖ}{k})^2}\biggr)\biggr)[/math]

[math]=[M_0+\frac{H_0}{k}]+\frac{1}{12}·\biggl(T_0-B_0+\frac{B}{1+(\frac{π}{6})^2}-\frac{1}{12·e^{12}}·\biggl(\frac{B}{1+(\frac{π}{6})^2}\biggr)\biggr)[/math]

[math]=[M_0+\frac{H_0}{k}]+0.0833(T_0-B_0+0.784·B)-5.12·10^{-7}·(T_0-B_0+0.784·B)[/math]

Los términos [math]0.0833(T_0-B_0+0.784·B)-5.12·10^{-7}·(T_0-B_0+0.784·B)[/math] resultan despreciables comparados con el primer término.

• Teniendo en cuenta que [math]k=1/4[/math]

[math][M_0+\frac{H_0}{k}]+\frac{1}{24k}·\biggl(T_0-B_0+\frac{B}{1+(\frac{ϖ}{k})^2}-\frac{1}{24·k·e^{24·k}}·\biggl(T_0-B_0+\frac{B}{1+(\frac{ϖ}{k})^2}\biggr)\biggr)[/math]

[math]=[M_0+\frac{H_0}{k}]+\frac{1}{6}·\biggl(T_0-B_0+\frac{B}{1+(\frac{π}{3})^2}-\frac{1}{6·e^{6}}·\biggl(\frac{B}{1+(\frac{π}{3})^2}\biggr)\biggr)[/math]

[math]=[M_0+\frac{H_0}{k}]+0.166(T_0-B_0+0.476·B)-4.13·10^{-4}·(T_0-B_0+0.476·B)[/math]

Los términos [math]0.166(T_0-B_0+0.476·B)-4.13·10^{-4}·(T_0-B_0+0.476·B)[/math] resultan despreciables comparados con el primer término, por lo que:

[math]T_{med}≅M_0+\frac{H_0}{k}[/math]

Y si [math]H_0=0[/math], entonces:

[math]T_{med}≅M_0[/math]

3.3 Desfase de la temperatura interior

Se pide comparar la función [math]F(t)[/math] con otra que nos demuestra que la variación senoidal de la temperatura del edificio tiene un retraso de [math]\alpha[/math] horas respecto a la variación exterior:

[math]F(t)=[1+(\frac{ϖ}{k})^2]^{-1}·cos⁡(ϖ·t-α)=\frac{cos⁡(ϖ·t)+\frac{ϖ}{k}·sen(ϖ·t)}{1+(\frac{ϖ}{k})^2}[/math]

Usamos las identidades trigonométricas:

[math]cos⁡(ϖ·t-α)=cos⁡(ϖ·t)·cos⁡(α)+sen(ϖ·t)·sen(α)[/math]

De tal manera que:

[math]\frac{(cos⁡(ϖ·t)·cos⁡(α)+sen(ϖ·t)·sen(α)}{(1+(\frac{ϖ}{k})^2}=\frac{cos⁡(ϖ·t)⁡+\frac{ϖ}{k}·sen(ϖ·t)}{1+(\frac{ϖ}{k})^2}[/math]

• [math]cos⁡(ϖ·t)·cos⁡(α)= cos(ϖ·t)·1→cos⁡(α)=1[/math]
• [math]sen(ϖ·t)·sen(α)=\frac{ϖ}{k}·sen(ϖ·t)→sen(α)=\frac{ϖ}{k}[/math]

Quedando:

[math] \frac{sen(α)}{cos⁡(α)}=\frac{\frac{ϖ}{k}}{1}⇒tg(α)=\frac{ϖ}{k}[/math]

3.4 Representación por método numérico e interpretación

Vamos a resolver numéricamente el probelma de Cauchy definido anteriormente, para ello tomaremos una temperatura interior inicial igual a 13°C y supondremos que la temperatura exterior varía con el tiempo adoptando la forma de una función senoidal [math]M(t)=7-5·cos(\frac{π}{12}·t).[/math] La representaremos mediante los métodos de Euler y de Runge-Kutta de orden 4 con longitudes de paso h=0.1, 0.01 y 0.001 para un lapso de tiempo de 24h.

Variación de la temperatura con el tiempo para paso h=0.1

Variación de la temperatura con el tiempo para paso h=0.01

Variación de la temperatura con el tiempo para paso h=0.001

% VARIACIÓN DE TEMPERATURA EN EL EDIFICIO SIN CALEFACCIÓN NI AIRE
% ACONDICIONADO
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc; clf; clear all;
% DATOS DE PARTIDA:
t0=0; tN=24; % Intervalo de tiempo en horas.
T0=input('Introduzca el valor de la temperatura interna inicial: ');       
h=input('Introduzca vector de longitudes de paso: ');
for i=1:length(h)
    h2=h(i);
% DISCRETIZACIÓN DE LA VARIABLE INDEPENDIENTE "t":
t=t0:h2:tN;           % Vector de tiempos.
f=inline('1/3*(7-5*cos((pi/12)*t)-T)+3','t','T');
figure(i)
T=zeros(1,length(t)); % Matriz "vacía" de la solución numérica.
T(1)=T0;
% FUNCIÓN TEMPERATURA EXTERIOR
M=7-5*cos((pi/12)*t);
% CÁLCULO NUMÉRICO Y GRÁFICAS:
%-------------------------- Método de Euler -------------------------------
TE=T;
for j=1:length(t)-1;
    TE(j+1)=TE(j)+h2*f(t(j),TE(j));
end
subplot(2,1,1) % Primera gráfica: Método de Euler.
hold on
plot(t,M,'b')  % Gráfica de temperaturas exteriores.
plot(t,TE,'r') % Gráfica de temperaturas interiores.
k=h(i);
legend('Exterior temperature','Indoor temperature','Location','best')
title(['Euler Method with  ' num2str(k) ''])
xlabel('Time (hours)'); ylabel('Temperature (ºC)')
hold off
%------------------------- Método de Runge-Kutta --------------------------
TRK=T;
for j=1:length(t)-1
  K1=f(t(j),TRK(j));
  K2=f(t(j)+h2/2,TRK(j)+(h2/2)*K1);
  K3=f(t(j)+h2/2,TRK(j)+(h2/2)+K2);
  K4=f(t(j)+h2,TRK(j)+h2*K3);
  TRK(j+1)=TRK(j)+(h2/6)*(K1+2*K2+2*K3+K4);
end
subplot(2,1,2)  % Segunda gráfica: Método de Runge-kutta de orden 4
hold on
plot(t,M,'b')   % Gráfica de temperaturas exteriores.
plot(t,TRK,'r') % Gráfica de temperaturas interiores.
legend('Exterior temperature','Indoor temperature','Location','best')
title(['Runge Kutta Method with ' num2str(k) ''])
xlabel('Time (hours)'); ylabel('Temperature (ºC)')
hold off
end

4 Efecto en la temperatura interior [math]T(t)[/math] del uso de calefacción y aire acondicionado [math](U(t)≠0)[/math]

En este apartado, el edificio tiene todos los componentes del apartado anterior pero además tiene calor producido por calefacción de la forma

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

donde [math]U(t)=T_D[/math] es la temperatura deseada y [math]K_u[/math] es la constante de transmisividad de la calefacción.

Es decir:

• [math]T(t)=[/math]Temperatura interior del edificio.
• [math]H(t)=[/math]Calor producido por personas y elementos del edificio[math]=H_0[/math]
• [math]U(t)=[/math]Calentamiento producido por la calefacción [math]=0[/math]
• [math]M(t)=[/math]Temperatura exterior [math]=M_0-B·cos⁡(ϖ·t)=M_0-B·cos⁡(π/12·t)[/math]
• [math]k= [/math]Constante que depende de las propiedades del edificio.
• [math]t_0=\frac{1}{k}=[/math]Constante del tiempo del edificio.
• [math]T_D=[/math]Temperatura deseada.
• [math]K_u=[/math]Constante de la calefacción
• [math]\frac{1}{k+K_D}=[/math]Constante de tiempo del edificio con calefacción

Es decir, la ecuación diferencial quedaría como:

\begin{equation*} \begin{cases} T'(t)=k(M(t)-T(t))+H(t)+U(t)=k(M_0-B·cos⁡(\frac{π}{12}·t)-T(t))+H_0+K_u·(T_D-T(t)) \\ T(0)=T_0 \end{cases} \end{equation*}

Que tiene como solución:

[math] \begin{equation*} T(t)=B_2-B_1·F_1(t)+C·e^{-k·t} \begin{cases} K_1=k+K_u \\ B_2=\frac{K_u·T_D+k·M_0+H_0}{K_1} \\ B_1=\frac{k·B}{K_1} \\ F_1 (t)=\frac{cos⁡(ϖ·t)+\frac{ϖ}{K_1}·sen(ϖ·t)}{1+(\frac{ϖ}{K_1})^2} \\ C=T_0-B_2+F_1 (0) \end{cases} \end{equation*} [/math]

4.1 Demostración del modelo

Para resolver esta ecuación diferencial usamos el principio de superposición, para poder usarlo, la ecuación diferencial debe ser lineal y los términos independientes deben satisfacer la misma ecuación diferencial homogénea.

[math]T'(t)=k(M_0-B·cos⁡(ϖ·t)-T(t))+H_0+K_u·(T_D-T(t))=k·M_0-k·B·cos⁡(ϖ·t)-k·T(t)+H_0+K_u·T_D-K_u·T(t)⟹[/math]

[math]\boxed{T'(t)+(k+K_u)·T(t)=k·M_0-k·B·cos⁡(ϖ·t)+H_0+K_u·T_D}[/math]

Dividimos la ecuación en dos, una con las constantes y otra con la función coseno:

[math] \begin{equation*} \begin{cases} T'(t)+(k+K_u)·T(t)=k·M_0+H_0+K_u·T_D \\ T'(t)+k·T(t)=-k·B·cos⁡(ϖ·t) \end{cases} \end{equation*} [/math]

1. Cálculo de la solución general de la homogénea:

Usamos el polinomio característico: [math]m+(k+K_u)=0⇔m=-(k+K_u)⇒S.G.H=c_1·e^{-(k+K_u)·t}[/math]

2. Cálculo de la solución particular de la ecuación con las constantes:

[math]T'(t)+(k+K_u)·T(t)=k·M_0+H_0+K_u·T_D[/math]

Conjunto asociado a la S.G.H: [math]S.G.H=c_1·e^{-(k+K_u)·t}[/math]

[math]S.G.H≡{e^{-(k+K_u)·t} }[/math]

Conjunto asociado al término independiente: [math]p_1 (t)=k·M_0+H_0+K_u·T_D[/math]

[math]p_1 (t)≡{1}\nsubseteq S.G.H[/math]

Comparando ambos conjuntos, podemos decir que la solución particular de la completa es del tipo:

\begin{equation*} T_p=A⇒ \begin{cases} T_p=A \\ T_p'=0 \end{cases} \end{equation*}

Sustituimos en nuestra ecuación diferencial:

[math]T'(t)+(k+K_u)·T(t)=k·M_0+H_0+K_u·T_D⟺0+k·A=k·M_0+H_0+K_u·T_D⇒[/math]

[math](k+K_u)·A=k·M_0+H_0+K_u·T_D⇒A=\frac{k·M_0+H_0+K_u·T_D}{k+K_u}[/math]

Es decir, nos quedaría una solución:

[math]T(t)=c_1·e^{-(k+K_u)·t}+\frac{k·M_0+H_0+K_u·T_D}{k+K_u}[/math]

3.Cálculo de la solución particular de la ecuación con el coseno:

[math]T'(t)+(k+K_u)·T(t)=-k·B·cos⁡(ϖ·t)[/math]

Conjunto asociado a la S.G.H::[math]S.G.H=c_1·e^{-(k+K_u)·t}[/math]

[math]S.G.H≡\{e^{-(k+K_u)·t} \}[/math]

Conjunto asociado al término independiente::[math]p_2 (t)=-k·B·cos⁡(ϖ·t)[/math]

[math]p_2 (t)≡\{cos⁡(ϖ·t),sen(ϖ·t)\}\nsubseteq S.G.H[/math]

Comparando ambos conjuntos, podemos decir que la solución particular de la completa es del tipo:

\begin{equation*} T_p=C·cos⁡(ϖ·t)+D·sen(ϖ·t)⇒ \begin{cases} T_p=C·cos⁡(ϖ·t)+D·sen(ϖ·t)\\ T_p'=-C·ϖ·sen(ϖ·t)+D·ϖ·cos⁡(ϖ·t) \end{cases} \end{equation*}

Sustituimos en nuestra ecuación diferencial:

[math]T'(t)+(k+K_u)·T(t)=-k·B·cos⁡(ϖ·t)⇒[/math]

[math]-C·ϖ·sen(ϖ·t)+D·ϖ·cos⁡(ϖ·t)+(k+K_u)·[C·cos⁡(ϖ·t)+D·sen(ϖ·t)]=-k·B·cos⁡(ϖ·t)⇒[/math]

[math][-ϖ·C+D·(k+K_u)]sen(ϖ·t)+[ϖ·D+C·(k+K_u)]cos⁡(ϖ·t)=-k·B·cos⁡(ϖ·t)[/math]

Subdividimos según las funciones:

• [math]sen(ϖ·t): -ϖ·C+B·(k+K_u)=0[/math]
• [math]cos⁡(ϖ·t): ϖ·B+A·(k+K_u)=-k·B[/math]

Y nos queda el sistema de ecuaciones siguiente:

[math] \begin{equation*} \begin{cases} -ϖ·C+D·(k+K_u)=0\\ ϖ·D+C·(k+K_u)=-k·B \end{cases} ⇒ \begin{pmatrix} -ϖ & k+K_u \\ k+K_u & ϖ \\ \end{pmatrix} · \begin{pmatrix} C \\ D \\ \end{pmatrix} = \begin{pmatrix} 0\\ -k·B\\ \end{pmatrix} \end{equation*} [/math]

Que resolvemos por el método de Cramer:

[math]C=\frac{det⁡(f_3,f_2)}{det⁡(f_1,f_2)}=\frac{k·(k+K_u)·B}{-(ϖ^2+(k+K_u)^2)}=-\frac{k·B·(k+K_u)}{ϖ^2+(k+K_u)^2}[/math]

[math]D=\frac{det⁡(f_1,f_3)}{det⁡(f_1,f_2)}=\frac{k·ϖ·B}{-(ϖ^2+(k+K_u)^2)}=-\frac{k·ϖ·B}{ϖ^2+(k+K_u)^2}[/math]

Es decir, nos quedaría una solución:

[math]T(t)=c_1·e^{-k·t}-\frac{k·B·(k+K_u)}{ϖ^2+(k+K_u)^2}·cos⁡(ϖ·t)-\frac{k·ϖ·B}{ϖ^2+(k+K_u)^2}·sen(ϖ·t)=[/math]

[math]c_1·e^{-k·t}-\frac{k·(k+K_u)·B·cos⁡(ϖ·t)+k·ϖ·B·sen(ϖ·t)}{ϖ^2+(k+K_u)^2}=[/math]

[math]c_1·e^{-k·t}-\frac{k·B}{k+K_u}·\frac{(cos⁡(ϖ·t)+\frac{ϖ}{k+K_u}·sen(ϖ·t)}{(\frac{ϖ}{k+K_u})^2+1}[/math]

Como se ha comentado antes, hemos aplicado el principio de superposición, por lo que la solución completa de nuestra ecuación queda como:

[math]\boxed{T(t)=\frac{k·M_0+H_0+K_u·T_D}{k+K_u}-\frac{k·B}{k+K_u}·\frac{(cos⁡(ϖ·t)+\frac{ϖ}{k+K_u}·sen(ϖ·t)}{(\frac{ϖ}{k+K_u})^2+1+c_1·e^{-(k+K_u)·t}}}[/math]

La condición inicial nos dice que T(0)=T_0, por lo que:

[math]T(0)=\frac{k·M_0+H_0+K_u·T_D}{k+K_u}-\frac{k·B}{k+K_u}·\frac{cos⁡(0)+\frac{ϖ}{k+K_u}·sen(0)}{(\frac{ϖ}{k+K_u})^2+1}+c_1·e^{-(k+K_u)·0} =T_0⟹[/math]

[math]\frac{k·M_0+H_0+K_u·T_D}{k+K_u}-\frac{k·B}{k+K_u}·\frac{1}{(\frac{ϖ}{k+K_u})^2+1}+c_1=T_0⟹[/math]

[math]\boxed{c_1=T_0+\frac{k·B}{k+K_u}·\frac{1}{(\frac{ϖ}{k+K_u})^2+1}-\frac{k·M_0+H_0+K_u·T_D}{k+K_u}}[/math]

Por lo que finalmente la Función solución queda como:

[math]\boxed{T(t)=\frac{k·M_0+H_0+K_u·T_D}{k+K_u}-\frac{k·B}{k+K_u}·\frac{cos⁡(ϖ·t)+\frac{ϖ}{k+K_u}·sen(ϖ·t)}{(\frac{ϖ}{k+K_u})^2+1}+\biggl(\frac{T_0+k·B}{k+K_u}·\frac{1}{(\frac{ϖ}{k+K_u})^2+1}-\frac{k·M_0+H_0+K_u·T_D}{k+K_u}\biggr)·e^{-(k+K_u)·t}}[/math]

Y comparándola con el enunciado:

[math] \begin{equation*} T(t)=B_2-B_1·F_1(t)+C·e^{-K_1·t} \begin{cases} K_1=k+K_u \\ B_2=\frac{K_u·T_D+k·M_0+H_0}{K_1} \\ B_1=\frac{k·B}{K_1} \\ F_1 (t)=\frac{cos⁡(ϖ·t)+\frac{ϖ}{K_1}·sen(ϖ·t)}{1+(\frac{ϖ}{K_1})^2} \\ C=T_0-B_2+F_1 (0) \end{cases} \end{equation*} [/math]

Nos queda de la misma forma.

4.2 Demostración[math]: B_2≅T_{med}≅T_D[/math] cuando desaparece el término exponencial

Se nos pide comprobar que [math]B_2[/math] es aproximadamente la temperatura media diaria y está muy próxima a [math]T_D[/math] si a los 30 minutos el término exponencial desaparece.

[math]T_{med}=\frac{1}{nº horas}·\sum_{i=1}^nT=\frac{1}{24}·\int_0^{24}T(t)·dt=\frac{1}{24}·\int_0^{24}[B_2-B_1·F_1(t)+C·e^{-k·t}]·dt=\frac{1}{24}·\int_0^{24}[B_2-B_1·F_1(t)]·dt[/math]

Para que sea más fácil de calcular y seguir, calcularemos la integral separando cada una de los términos:

[math]\frac{1}{24}·\int_0^{24}[B_2-B_1·F_1(t)]·dt=\frac{1}{24}·(\int_0^{24}B_2·dt+\int_0^{24}-B_1·F_1(t)·dt)[/math]

• Término [math]B_2[/math]:

[math]\int_0^{24}B_2·dt=B_0\int_0^{24}dt=B_2[t]_0^{24}=24·B_0=24·[\frac{K_u·T_D+k·M_0+H_0}{K_1}][/math]

[math]\boxed{\int_0^{24}B_0·dt=24·[\frac{K_u·T_D+k·M_0+H_0}{K_1}]}[/math]

• Función [math]F(t)[/math]:

[math]\int_0^{24}-B·F(t)·dt=-\frac{k·B}{K_1·(1+(\frac{ω}{K_1})^2)}\int_0^{24}cos(ω·t)+\frac{ω}{K_1}sen(ω·t)·dt[/math]

[math]=\biggl\{k=-\frac{k·B}{K_1·(1+(\frac{ω}{K_1})^2}\biggr\}[/math]

[math]=K\biggl[\int_0^{24}cos(ω·t)+\int_0^{24}\frac{ω}{K_1}sen(ω·t)·dt\biggr][/math]

[math]=K\biggl[\frac{1}{ϖ}\int_0^{24}ω·cos(ω·t)-\frac{1}{K_1}\int_0^{24}-ω·sen(ω·t)·dt\biggr][/math]

[math]=K·(\frac{1}{ϖ} [sen(ϖ·t)]_0^{24}+\frac{1}{K_1} [cos⁡(ϖ·t) ]_{0^24})=\{ϖ=\frac{π}{12}\}=[/math]

[math]=K·(\frac{1}{ϖ} (sen(2·π)-sen(0))+\frac{1}{K_1} (cos(2·π)-cos(0)))[/math]

[math]=K·(\frac{1}{ϖ} (0-0)+\frac{1}{K_1} (1-1))=0[/math]

[math]\boxed{\int_0^{24}-B_1·F_1(t)·dt=0}[/math]

Vamos a comparar las constantes:

• [math]k=[/math]Transmisividad de los muros
• [math]K_u=[/math]Transmisividad de los radiadores

La transmisividad de los muros es mucho menor que los de los radiadores, ya que se quiere conseguir que la temperatura exterior no varíe la interior y a su vez se pretende que los radiadores puedan dejar fluir todo el calor que sea posible. Resumiendo:

[math]K_u≫k⟹K_1=k+K_u⟹K_1≫k⟹K_1≅K_u[/math]

[math]T_{med}=\frac{K_u·T_D+k·M_0+H_0}{k+K_u}=\frac{K_u·T_D}{k+K_u}+{k·M_0}{k+K_u}+\frac{H_0}{k+K_u}≅T_D+0+0[/math]

[math]T_{med}≅T_D[/math]

4.3 Representación por método numérico e interpretación

Vamos a resolver numéricamente el probelma de Cauchy definido anteriormente, para ello tomaremos una temperatura interior inicial igual a 24°C y supondremos que la temperatura exterior varía con el tiempo adoptando la forma de una función senoidal [math]M(t)=4-5·cos(\frac{π}{12}·t).[/math] Además, tendremos un calor inicial producido por las personas y otros elementos del edificio de 3°C. La representaremos mediante los métodos de Euler y de Runge-Kutta de orden 4 con longitudes de paso h=0.1, 0.01 y 0.001 para un lapso de tiempo de 24h.

Variación de la temperatura con el tiempo para distintos pasos usando calefacción y aire acondicionado

% VARIACIÓN DE TEMPERATURA EN EL EDIFICIO CON CALEFACCIÓN O AIRE
% ACONDICIONADO
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc; clf; clear all;
% DATOS DE PARTIDA:
t0=0; tN=24; % Intervalo de tiempo en horas.
T0=input('Introduzca el valor de la temperatura interna inicial: ');       
h=input('Introduzca vector de longitudes de paso: ');
for i=1:length(h)
    h2=h(i);
% DISCRETIZACIÓN DE LA VARIABLE INDEPENDIENTE "t":
t=t0:h2:tN;           % Vector de tiempos.
f=inline('1/3*(4-5*cos((pi/12)*t)-T)+3+(7/8)*(22-T)','t','T');
figure(i)
T=zeros(1,length(t)); % Matriz "vacía" de la solución numérica.
T(1)=T0;
% FUNCIÓN TEMPERATURA EXTERIOR
M=4-5*cos((pi/12)*t);
% CÁLCULO NUMÉRICO Y GRÁFICAS:
%-------------------------- Método de Euler -------------------------------
TE=T;
for j=1:length(t)-1;
    TE(j+1)=TE(j)+h2*f(t(j),TE(j));
end
subplot(2,1,1) % Primera gráfica: Método de Euler.
hold on
plot(t,M,'b')  % Gráfica de temperaturas exteriores.
plot(t,TE,'r') % Gráfica de temperaturas interiores.
k=h(i);
legend('Exterior temperature','Indoor temperature','Location','best')
title(['Euler Method with  ' num2str(k) ''])
xlabel('Time (hours)'); ylabel('Temperature (ºC)')
hold off
%------------------------- Método de Runge-Kutta --------------------------
TRK=T;
for j=1:length(t)-1
  K1=f(t(j),TRK(j));
  K2=f(t(j)+h2/2,TRK(j)+(h2/2)*K1);
  K3=f(t(j)+h2/2,TRK(j)+(h2/2)+K2);
  K4=f(t(j)+h2,TRK(j)+h2*K3);
  TRK(j+1)=TRK(j)+(h2/6)*(K1+2*K2+2*K3+K4);
end
subplot(2,1,2)  % Segunda gráfica: Método de Runge-kutta de orden 4
hold on
plot(t,M,'b')   % Gráfica de temperaturas exteriores.
plot(t,TRK,'r') % Gráfica de temperaturas interiores.
legend('Exterior temperature','Indoor temperature','Location','best')
title(['Runge Kutta Method with ' num2str(k) ''])
xlabel('Time (hours)'); ylabel('Temperature (ºC)')
hold off
end

5 Planteamiento del edificio como dos compartimentos

El apartado D nos dice que el edificio ahora tiene dos zonas diferenciadas como se ve en la figura.

Modelo de edificio con dos zonas diferenciadas.

Cada zona tiene unas características:

• [math]T_a(t)=[/math]Temperatura en el interior de la zona A
• [math]H_a(t)=[/math]Calor producido por las personas y máquinas de la zona A
• [math]U_a(t)=[/math]Calentamiento producido por la calefacción de la zona A
• [math]k_e1=[/math]Constante de transmisividad de la zona A al exterior
• [math]T_b(t)=[/math]Temperatura en el interior de la zona B
• [math]H_b(t)=[/math]Calor producido por las personas y máquinas de la zona B
• [math]U_b(t)=[/math]Calentamiento producido por la calefacción de la zona B
• [math]k_e1=[/math]Constante de transmisividad de la zona B al exterior
• [math]k_i=[/math]constante de transmisividad entre las zonas A y B

5.1 Determinación del problema de Cauchy

La zona A sóla tendría una ecuación diferencial del tipo:

\begin{equation*} \begin{cases} T_a'(t)=k_{e1}(M(t)-T_a(t))+H_a(t)+U_a(t)\\ T_b(0)=T_(b,0) \end{cases} \end{equation*}

La zona B sóla tendría una ecuación diferencial del tipo:

\begin{equation*} \begin{cases} T_b'(t)=k_{e2}(M(t)-T_b (t))+H_b(t)+U_b(t)\\ T_b (0)=T_(b,0) \end{cases} \end{equation*}

Ahora las zonas no están separadas sino que interactuan entre sí, de la misma manera que una zona interactua con el exterior, es decir, que la zona A tendrá ahora un término nuevo como [math]k_i·(T_b(t)-T_a(t))[/math]. De la misma manera, la zona B tendra otro nuevo término [math]k_i·(T_a (t)-T_b(t))[/math].

Queda un sistema de ecuaciones:

\begin{equation*} \begin{cases} T_a'(t)=k_{e1}(M(t)-T_a(t))+H_a(t)+U_a(t)+k_i·(T_b(t)-T_a(t))\\ T_b'(t)=k_{e2}(M(t)-T_b(t))+H_b(t)+U_b(t)+k_i·(T_a(t)-T_b(t))\\ T_a(0)=T_(a,0); T_b(0)=T_(b,0) \end{cases} \end{equation*}

Donde:

[math]M(t)=2-7·cos⁡(\frac{π·t}{12}); y k_i=1/2[/math]

Zona A:

[math]H_a(t)=10\frac{t}{1+t}[/math]

[math]U_a(t)=\frac{5}{3}(22-T_a(t))[/math]

[math]k_{e1}=\frac{1}{4}[/math]

Zona B:

[math]H_b(t)=4·cos⁡(\frac{π·t}{6}[/math]

[math]U_b(t)=\frac{13}{7}(23-T_b(t)[/math]

[math]k_{e2}=\frac{1}{5}[/math]

Sustituimos y queda un sistema:

\begin{equation*} \begin{cases} T_a'(t)=-(\frac{5}{3}+k_{e1}+k_i)·T_a(t)+k_i·T_b(t)+k_{e1}(2-7·cos⁡(\frac{π·t}{12})+\frac{110}{3}+\frac{10·t}{1+t}) \\ T_b'(t)=-(\frac{13}{7}+k_{e2}+k_i)·T_b(t)+k_i·T_a+k_{e2}·(2-7·cos⁡(\frac{π·t}{12})+\frac{299}{7}+4·cos⁡(\frac{π·t}{6})) \\ T_a(0)=20 \\ T_b(0)=18 \end{cases} \end{equation*}

5.2 Representación por método numérico e interpretación

Para nuestro programa de Matlab, el sistema queda de la forma:

[math]T_a'(t)=-\frac{29}{12}·T_a(t)+\frac{1}{2}·T_b(t)+\frac{1}{4}(2-7·cos⁡(\frac{π·t}{12})+\frac{110}{3}+\frac{10·t}{1+t}[/math]

[math]T_b'(t)=-\frac{179}{70}·T_b(t)+\frac{1}{2}·T_a+1/5·(2-7·cos⁡(\frac{π·t}{12})+\frac{299}{7}+4·cos⁡(\frac{π·t}{6}))[/math]

Bucle de Euler:

La función que debemos meter en el programa es el vector:

[math]f= \begin{pmatrix} -\frac{9}{12}·T_a(t)+1/2·T_b(t)+\frac{1}{4}(2-7·cos⁡(\frac{π·t}{12})+\frac{110}{3}+\frac{10·t}{1+t})\\ -\frac{179}{70}·T_b(t)+\frac{1}{2}·T_a+\frac{1}{5}·(2-7·cos⁡(\frac{π·t}{12})+\frac{299}{7}+4·cos⁡(\frac{π·t}{6})\\ \end{pmatrix} [/math]

Siendo el bucle:

[math]T(:,k+1)=T(:,k)+h·f(t(k),T(:,k))[/math]

Bucle de Euler implícito:

La función que debemos meter en el programa es el vector:

[math]f= \begin{pmatrix} -\frac{9}{12}·T_a(t)+1/2·T_b(t)+\frac{1}{4}(2-7·cos⁡(\frac{π·t}{12})+\frac{110}{3}+\frac{10·t}{1+t})\\ -\frac{179}{70}·T_b(t)+\frac{1}{2}·T_a+\frac{1}{5}·(2-7·cos⁡(\frac{π·t}{12})+\frac{299}{7}+4·cos⁡(\frac{π·t}{6})\\ \end{pmatrix} [/math]

Pasamos a matricial:

[math]T(t)=A·T+B= \begin{pmatrix} -(\frac{5}{3}+k_{e1}+k_i) & k_i\\ k_i & -(\frac{13}{7}+k_{e2}+k_i)\\ \end{pmatrix} · \begin{pmatrix} T_A\\ T_B\\ \end{pmatrix}+ \begin{pmatrix} k_{e1}·(2-7·cos⁡(\frac{π·t}{12})+\frac{110}{3}+\frac{10·t}{1+t})\\ k_{e2}·(2-7·cos⁡(\frac{π·t}{12})+\frac{299}{7}+4·cos⁡(\frac{π·t}{6})\\ \end{pmatrix} [/math]

Calculamos el bucle:

[math]T_{n+1}=T_n+h·f(t_{n+1},T_{n+1})⟷T_{n+1}=T_n+h·[A·T_{n+1}+B_{n+1}][/math]

[math]→T_{n+1}-h·A·T_{n+1}=T_n+h·B_{n+1}[/math]

[math]→T_{n+1}·(I-h·A)=T_n+h·B_{n+1}[/math]

[math]→T_{n+1}=(T_n+h·B_{n+1})/((I-h·A))[/math]

Quedando:

[math]T(:,k+1)=(I-h·A)\\(T(:,k)+h·B(:,k+1))[/math]

Bucle de Runge-Kutta:

La función que debemos meter en el programa es el vector:

[math]f= \begin{pmatrix} -\frac{9}{12}·T_a(t)+1/2·T_b(t)+\frac{1}{4}(2-7·cos⁡(\frac{π·t}{12})+\frac{110}{3}+\frac{10·t}{1+t})\\ -\frac{179}{70}·T_b(t)+\frac{1}{2}·T_a+\frac{1}{5}·(2-7·cos⁡(\frac{π·t}{12})+\frac{299}{7}+4·cos⁡(\frac{π·t}{6})\\ \end{pmatrix} [/math]

Siendo el bucle:

[math]K1=f(t(k),T(:,k));[/math]

[math]K2=f(t(k)+\frac{1}{2}*h,T(:,k)+1/2*K1*h);[/math]

[math]K3=f(t(k)+{1}{2}*h,T(:,k)+1/2*K2*h);[/math]

[math]K4=f(t(k)+h,T(:,k)+K3*h);[/math]

[math]T(:,k+1)= T(:,k)+\frac{h}{6}*(K1+2*K2+2*K3+K4);[/math]

% EDIFICIO CON DOS ZONAS DIFERENCIADAS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc; clf; clear all;
% DATOS DE PARTIDA:
ne=2;        % Número de ecuaciones
t0=0; tN=24; % Intervalo de tiempo en horas.
T0=[20,18];  % Vector temp. iniciales de las zonas A y B respectivamente.
ke1=1/4;     % Transmisividad térmica externa en la zona A.
ke2=1/5;     % Transmisividad térmica externa en la zona B.
ki=1/2;      % Transmisividad térmica interna entre zonas.
    
% PLANTEAMIENTO DE LA FUNCIÓN DERIVADA:
fun=input('Introduzca el vector función y´=f(t,T)= ','s');
% DISCRETIZACIÓN DE LA VARIABLE INDEPENDIENTE "t":
h=input('Introduzca el tamaño de paso h: ');
for i=1:length(h)
    h2=h(i);
 N=round((tN-t0)/h2);   % Número de subintervalos.
 t=t0:h2:tN;            % Vector de tiempo.
 f=inline(fun,'t','T');
 T=zeros(ne,length(t)); % Matriz "vacía" de la solución numérica.
 T(:,1)=T0';
 figure(i)
%--------------------------------------------------------------------------
% FUNCIÓN TEMPERATURA EXTERIOR:
  M=2-7*cos(pi*t/12);
% CÁLCULO NUMÉRICO Y GRÁFICAS:
%--------------------------- Método de Euler ------------------------------
for k=1:N
    T(:,k+1)=T(:,k)+h2*f(t(k),T(:,k));
end
subplot(2,2,1)     % Primera gráfica: Método de Euler.
hold on
plot(t,M,'b')      % Gráfica de temperaturas exteriores.
plot(t,T(1,:),'r') % Gráfica de temperaturas en la zona A.
plot(t,T(2,:),'g') % Gráfica de temperaturas en la zona B.
j=h(i);
legend('Exterior temperature','Temperature A','Temperature B','Location','best')
title(['Zonal temperature (Euler Method) with ' num2str(j) ''])
xlabel('Time (hours)'); ylabel('Temperature (ºC)')
hold off
%---------------------- Método de Euler implícito -------------------------
I=eye(ne);
for k=1:N
    A=[-(5/3+ke1+ki),ki;ki,-(13/7+ke2+ki)];
    B=[ke1*(2-7*cos(pi*t/12))+110/3+10*t./(1+t);
       ke2*(2-7*cos(pi*t/12))+299/7+4*cos(pi*t/6)];
    T(:,k+1)=(I-h2*A)\(T(:,k)+h2*B(:,k+1));
end
subplot(2,2,2)     % Segunda gráfica: Método de Euler Implícito.
hold on
plot(t,M,'b')      % Gráfica de temperaturas exteriores.
plot(t,T(1,:),'r') % Gráfica de temperaturas en la zona A.
plot(t,T(2,:),'g') % Gráfica de temperaturas en la zona B.
j=h(i);
legend('Exterior temperature','Temperature A','Temperature B','Location','best')
title(['Zonal temperature (Backward Euler Method) with ' num2str(j) ''])
xlabel('Time (hours)'); ylabel('Temperature (ºC)')
hold off
%----------------------- Método de Runge-Kutta ----------------------------
for k=1:N
    K1=f(t(k),T(:,k));
    K2=f(t(k)+1/2*h2,T(:,k)+1/2*K1*h2);
    K3=f(t(k)+1/2*h2,T(:,k)+1/2*K2*h2);
    K4=f(t(k)+h2,T(:,k)+K3*h2);
    T(:,k+1)=T(:,k)+h2/6*(K1+2*K2+2*K3+K4);
end
subplot(2,2,3)     % Tercera Gráfica: Método de Runge-Kutta de orden 4.
hold on
plot(t,M,'b')      % Gráfica de temperaturas exteriores.
plot(t,T(1,:),'r') % Gráfica de temperaturas en la zona A.
plot(t,T(2,:),'g') % Gráfica de temperaturas en la zona B.
j=h(i);
legend('Exterior temperature','Temperature A','Temperature B','Location','best')
title(['Zonal temperature (Runge-Kutta Method) with ' num2str(j) ''])
xlabel('Time (hours)'); ylabel('Temperature (ºC)')
hold off
end
%--------------------------------------------------------------------------

Visualización de la temperatura en ambas zonas con distintos métodos y paso 0.1

Visualización de la temperatura en ambas zonas con distintos métodos y paso 0.01

Visualización de la temperatura en ambas zonas con distintos métodos y paso 0.001