Desarrollo de un modelo espacio-temporal de la temperatura del aire en el interior de un edificio de múltiples zonas

Equation Chapter 1 Section 1 Trabajo de Fin de Máster Máster en Ingeniería Industrial Desarrollo de un modelo espacio-temporal de la temperatura del aire en el interior de un edificio de múltiples zonas Dpto. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Autor: Alfonso Domínguez Galindo Tutora: Amparo Núñez Reyes Sevilla, 2021 iii Proyecto Fin de Carrera Ingeniería Industrial Desarrollo de un modelo espacio-temporal de la temperatura del aire en el interior de un edificio de múltiples zonas Autor: Alfonso Domínguez Galindo Tutora: Amparo Núñez Reyes Profesora titular Dpto. de Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 v Proyecto Fin de Carrera: Desarrollo de un modelo espacio-temporal de la temperatura del aire en el interior de un edificio de múltiples zonas Autor: Alfonso Domínguez Galindo Tutora: Amparo Núñez Reyes El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2021 El Secretario del Tribunal vii A mi familia ix Resumen La eficiencia energética en edificios ha cobrado significado conforme se han desarrollado herramientas computacionales capaces de desarrollar la matemática que rige las ecuaciones que gobiernan las transferencias de calor. De la mano de los controladores predictivos basados en modelo, se plantea en este trabajo la concepción de un simulador, por el cual se pueda obtener una predicción espacial de la temperatura interior en un volumen de aire. Ciertas hipótesis y simplificaciones serán tenidas en cuenta a la hora de abordar el comportamiento físico y termodinámico del conjunto considerado. Mediante el software de cálculo matemático Matlab, se desarrollarán dos modelos para el mismo modelo físico: uno simple, donde una única temperatura define el volumen interior a considerar, y uno tridimensional, donde se ha aplicado al volumen en cuestión una malla tridimensional, con el fin de poder caracterizar el aire espacialmente y valorar resultados entre modelos. Finalmente, tras los resultados se demuestra la conveniencia de la caracterización espacial del elemento a considerar sujeto a un análisis térmico. 1 INTRODUCCIÓN 1.1. Motivación y objetivos La temperatura del aire interior de los edificios residenciales, oficinas o cualquier otro tipo de cubiertas donde la presencia humana pueda ser considerada, comprende un tema ampliamente estudiado por la comunidad científica con anterioridad. Variables como la propia temperatura, ya sea media o local, el confort térmico o la humedad son factores que están estrechamente relacionados con el consumo energético del ser humano desde que se dispone de sistemas de climatización en interiores. Es por ello que se persigue la obtención de simuladores y controladores robustos que puedan garantizar una correcta optimización en tiempos de uso y potencia en dichos sistemas de aclimatamiento de interiores, a fin de reducir la potencia consumida para usos comunes y aumentar tiempos de residencia de la energía en edificios. En definitiva, se busca aumentar en la medida de lo posible la eficiencia energética en el conjunto sistema, que comprende la estructura en sí con su inercia térmica junto con los sistemas de calefacción y/o climatización. El comportamiento térmico de edificios y análisis térmico de materiales aislantes en las envolventes se vuelve crucial en términos energéticos, debido a que permite acumular la energía almacenada en interiores en climas fríos, o sesga la radiación solar que incide sobre envolventes transparentes como, por ejemplo, son las ventanas de doble acristalamiento bajo emisivo para otros climas más cálidos, según el I.D.A.E. (2007). Por otra parte, una elección adecuada de los sistemas de climatización de interiores comprende otro factor importante en el tema que se aborda, puesto que con configuraciones distintas de espacio y consumo esperado pueden ser convenientes sistemas de calefacción portátiles, suelo radiante u otros sistemas que ofrecen gradientes de temperatura y flujos netos de calor distintos. Los datos expuestos a continuación dan un atisbo de la importancia de la gestión eficiente de la energía en los hogares, la cual está comprendiendo un aumento paulatino debido al calentamiento global y a los fenómenos meteorológicos cada vez más frecuentes. Acorde a Lin y Lin (2017), la calefacción de espacios residenciales ocupa la mayor parte de la energía consumida por los edificios en China. En el norte de China, el consumo de energía en la calefacción de espacios representa alrededor del 40% del consumo energético total en edificios (Tsinghua University, 2011). El suministro de calefacción es una de las necesidades más básicas de los habitantes del norte de China durante el invierno, y la demanda energética en este ámbito se ha visto incrementada en la última década si comparamos esta subida con la que se ha producido en la industria de este país. En el año 2011 se apuntaba el consumo de energía para calefacción de hogares en un 30% sobre el total de energía consumida en China, notándose una previsión alcista en años posteriores. 2 Figura 1-1. Suministro de calefacción central en China durante 1985-2014. (Lin y Lin, 2017) El autor relaciona esta subida de consumo energético en calefacción en hogares con factores como el crecimiento económico del país (PIB nacional), aumento del poder adquisitivo del individuo, mejora de la calidad de vida y a la densidad de población, entre otros. Por otra parte, gracias a los datos ofrecidos por la Oficina de Estadística de la Unión Europea (2020), sabemos que el consumo energético en España durante el año 2018 en hogares, comercios y otros servicios afronta el 31,93% de la energía total consumida en el país durante este período. En este sentido, se están llevando a cabo en la actualidad programas de fomento de mejoras de eficiencia energética en viviendas unifamiliares recogidas en el Real Decreto 106/2018, de 9 de marzo, por el que se regula el Plan Estatal de Vivienda 2018-2021, contenidas en el Plan Estatal de Vivienda 2018-2021, con el fin de reducir la factura energética de las familias, y las emisiones de gases de efecto invernadero. En efecto, el término de eficiencia energética es la que ha impulsado principalmente los estudios al respecto, enfocando el trabajo en la obtención de modelos predictivos donde se puedan implementar técnicas de control sobre los sistemas de calefacción en edificios residenciales principalmente. Estos se derivan de la necesidad de reducir el consumo de energía, integrar suministros de energía limpia y reducir el impacto medioambiental (Clarke y Hensen, 2015). Se pueden diferenciar dos disciplinas a priori en lo que respecta al objetivo final del control de estos sistemas: el modelo para simulación y los métodos de control aplicados al sistema. 3 1.1.1 Building Performance Simulation En lo referente a la simulación de variables de comfort en espacios residenciales, existe el término BPS, acrónimo de Building Performance Simulation, el cual hace referencia a la utilización de sistemas computacionales para simular el comportamiento de edificios en múltiples escenarios y su control. Este campo lleva desarrollándose desde los 60, siendo el principal enfoque en la última década el cálculo de cargas y el análisis térmico (Spitler, 2006). Con el paso del tiempo y la enorme evolución de las tecnologías computacionales se han podido alcanzar límites más ambiciosos en el campo de la simulación, en términos de número de variables simultáneas y su integración en modelos, como pueden ser aproximaciones más exigentes de transferencia de calor y masa, mejores aproximaciones en la definición de corrientes de aire interiores o condiciones climatológicas externas (Clarke y Hensen, 2015). Como apuntan Clarke y Hensen (2015), el objetivo fundamental del BPS es reflejar de la forma más integral y aproximada el comportamiento dinámico de estos sistemas, que por lo general son procesos con variables fuertemente dependientes entre sí y sujetos a no linealidades. Los autores remarcan los tres aspectos más importantes de la disciplina: • Integridad física de las variables involucradas. En la mayoría de ocasiones se criban algunas variables a fin de no mermar tiempos computacionales o la efectividad de los modelos en ciertos campos. La complejidad geométrica de elementos presentes, que dificulta su modelado espacial, o la aproximación estocástica de la ocupación humana son otros ejemplos. Estos errores (en su sentido matemático) deben ser entendidos y analizados en la obtención de los resultados de las simulaciones. • Simultaneidad de diferentes dominios. Se refiere al avance en el acoplamiento de dominios como la temperatura y la iluminación, o el calor y el flujo de aire. Estos avances son lentos y aún no se ha desarrollado un modelado integral a grandes rasgos. Se tiende a simplificar en este tipo de modelos en términos de transferencia de calor, siendo recurrentes los análisis térmicos en régimen permanente, o la transferencia de calor y masa en una sola dimensión al hablar de paredes. • Integración del proceso de diseño. Se busca de manera continuada la posibilidad de importación y exportación de datos con otros softwares especializados de CAD, estimación de costes o análisis estructural, lo que hace ganar valor al simulador de cara al cliente/usuario. Figura 1-2. Simulación térmica de edificio en software especializado. (PRWeb, 2021) 4 1.1.2 Control Predictivo Basado en Modelo La metodología del Modelo de Control Predictivo es un término extendido y amplio que de por sí no establece una estrategia de control definida, si no una gran variedad de posibilidades en cuanto a métodos de control se refiere, en base a un modelo predefinido (Camacho y Bordons, 2000). En el tema de estudio, esta metodología juega un papel importante, ya que conociendo el comportamiento del sistema mediante un modelo (campo del BPS), se pueden obtener señales de control minimizando una función objetivo. La referencia a este método de control se debe a que el modelo objeto de este trabajo está enfocado a la implementación en este tipo de metodología de control (véase objetivos). Camacho y Bordons (2000) ofrecen definiciones y descripciones básicas al principio de su libro, en lo referente a la estructura de esta metodología y la estrategia de control seguida en el MPC. 1.1.2.1 Estructura del MPC Esta metodología se compone de los siguientes bloques: • Un modelo matemático para predecir el comportamiento o las variables salida del proceso en instantes de tiempo futuros, previa definición de un horizonte de control. • Una función objetivo a minimizar. Rossiter (2003) declara al respecto: “si la función de costes es correcta, la estabilidad y el ajuste se cuidan solos, ya que, por definición, se está optimizando un coste que sólo puede ser pequeño para un buen rendimiento” (p. 25). • Una estrategia control con horizonte deslizante, por la cual el horizonte de control se va desplazando al futuro en cada instante de tiempo. Aquí, las señales de control calculadas en el primer paso serán aplicadas para calcular la sucesiva predicción, y las demás serán deshechadas. 1.1.2.2 Estrategia del MPC La metodología de los controladores de esta familia comprende los siguientes pasos: 1. En cada instante de tiempo t se predicen las salidas del proceso y(t + k | t) sobre un horizonte de tiempo N donde k=1…N. Para ello se hace uso de un modelo del proceso, donde se tienen en cuenta las entradas y salidas anteriores al instante t, y a las señales de control futuras u(t + k | t) donde k=0… N-1. Estas señales de control futuras son enviadas al sistema. 2. Las señales de control futuras u(t + k | t) son calculadas bajo un criterio de optimización en aras de mantener las salidas del proceso lo más cerca posible de una trayectoria de referencia w(t + k), la cual puede ser el propio setpoint o una aproximación al mismo, apuntan los autores. Estos criterios de optimización suelen tomar la forma de funciones cuadráticas de los errores entre la trayectoria de referencia y la salida calculada. 3. La señal de control u(t | t) se envía al proceso, mientras que las señales calculadas en instantes de tiempo posteriores son rechazadas, puesto que el proceso en este momento se repetirá desde el primer paso, y otro conjunto nuevo de señales de control serán calculadas yo la salida del proceso y(t + 1) es ya sabida. Debido a esta estrategia de retroalimentación, las señales de control para un mismo instante futuro t serán distintas en cada iteración, puesto que se irán recalculando. 5 Figura 1-3. Estructura básica del MPC (Camacho y Bordons, 2000) En la Figura 1-3 se ve claramente la metodología en cuestión. El modelo recibe en primer lugar las entradas y salidas del paso anterior (o set-point a su defecto), y obtiene una predicción futura en base al cálculo matemático. Esta predicción se coteja con una trayectoria de referencia y se obtiene por ende cierto error. Tras esto, se aplica una estrategia de control definida en base a este error, una función objetivo y restricciones a considerar, para mandar de nuevo al modelo una nueva señal de control y volver a repetir el proceso. 1.1.2.3 Principales ventajas y desventajas de la metodología Camacho y Bordons (2000) señalan de forma clara los principales puntos a favor del control basado en modelos frente a otros controladores: • Los conceptos implementados son intuitivos y resulta asequible para personas con nociones básicas de control. • Puede ser usado en una gran variedad de procesos de diferentes complejidades. Algunas aplicaciones se encuentran en el campo de la automoción para el control en trenes de potencia (Raković y Levine, 2018) o en el sector petroquímico como pueden ser las columnas de destilación de crudo (Buck et al., 2011). • Permite manejar sistemas tanto de una sola variable como multivariables, principal ventaja respecto a los controladores PID, cuya implementación en sistemas MIMO (‘multi-input-multi-output’) es complicada. • La introducción del control anticipativo compensa perturbaciones medibles. • Trabaja muy cercanamente a los límites reales (restricciones) impuestos al proceso, permitiendo ello maximizar rendimientos. Entre las principales desventajas frente a controladores más simples, encontramos que la derivación de la ley de control es más compleja que en los clásicos controladores PID. El tiempo computacional será mayor, más si tenemos en cuenta que el control adaptativo necesita un esfuerzo computacional en cada instante de tiempo, puesto que la dinámica del sistema es cambiante; y, por otra parte, las restricciones del proceso son implementadas. Así mismo, se hace evidente que una buena obtención matemática del proceso a modelar requiere conocimiento previo sobre el mismo, lo cual cuantificará de forma directa y significativa las discrepancias entre el proceso real y el modelo. 6 2 ESTADO DEL ARTE Y OBJETIVOS Se describen a continuación algunos artículos científicos en los que se basa el presente trabajo, a fin de conocer qué hay hecho hasta el momento y la dirección que toma el mismo. En el artículo de Mendes et al. (2001) se realiza un modelo matemático aplicado al análisis térmico de edificios y al diseño de sistemas de control del mismo. Se implementa en el software Matlab/Simulink de forma práctica un modelo dinámico y multinodal el transitorio de la temperatura media del aire interior de un edificio en un día frío en el centro de Brasil. Las temperaturas de trabajo que caracterizan volúmenes de masa completos son temperaturas medias, y la implementación de las paredes es capacitiva. La temperatura exterior del aire se modeló de forma sinusoidal, incluyendo un parámetro de temperatura equivalente aire-sol, donde el efecto de la radiación solar en la envolvente exterior está intrínseco. Puesto que el artículo se centra principalmente en la obtención del modelo, no se indaga profundamente en metodologías de control, y sólo hay implementado un actuador on-off sobre la fuente de calor del espacio, un radiador convectivo de aceite de 5kW. Li y Zang (2020) en su artículo realizan un modelo matemático para caracterizar la temperatura interior de un edificio multizona de forma similar al anterior. En éste, la radiación no está contemplada, y la fuente de calor reside en la instalación de radiadores de agua caliente, con la caracterización fluidomecánica que conlleva el sistema. El artículo se centra en el desarrollo de la metodología MPC para el control de la temperatura interior, mediante algoritmos de optimización híbridos para edificios de una zona, y MPC distribuido posteriormente para un edificio multizona, declarando diferencias sustanciales frente al MPC descentralizado. Los autores Moroşan et al. (2010) presentan en su artículo una estructura de control predictivo para la regulación de la temperatura en edificios. De nuevo se vuelve a enfocar la diferencia entre los modelos de zona única versus multizona, y las diferencias entre la estrategia del MPC distribuido y el centralizado, siendo este último el más conveniente cuando existen paredes intermedias donde se produce transferencia de calor, los cuales son nodos de acoplamiento. El punto más importante desde el punto de vista de este trabajo es la implementación de un horario de ocupación humana en las zonas a controlar. Los autores utilizan estos perfiles de ocupación preestablecidos para reducir el consumo de energía sin afectar al comfort térmico en estas franjas horarias. 2.1. Objetivos En los modelos consultados, se simplifica el cálculo del transitorio del aire interior de los edificios en términos de una sola variable, es decir, una sola temperatura; la temperatura media la cual caracteriza el volumen completo de aire que comprende el espacio interior. A priori, el emplazamiento en el espacio de las fuentes de calor, el de las ventanas y otros elementos constructivos; elementos disipadores de calor o las personas, hacen pensar que la influencia sobre el gradiente térmico en los espacios será importante. Esta caracterización espacio-temporal de los modelos será la principal vía de desarrollo para el presente trabajo. El objetivo principal es desarrollar un modelo matemático que pueda implementarse en futuros controladores, donde esta propiedad del espacio esté contemplada. De esta forma, se podrán obtener los datos de temperatura en ciertos puntos concretos del espacio donde estén colocados sensores o similares. Además, se busca también una correcta caracterización del gradiente térmico en el espacio para construir de forma aproximada set-points con datos distribuidos espacialmente. 7 3 MODELOS: CONSIDERACIONES GENERALES En este capítulo se explicarán las consideraciones comunes a los modelos creados, referentes a diversas cuestiones diferentes entre sí. El hecho de caracterizar el comportamiento de un fluido compresible como el aire supone recurrir constantemente a suposiciones e hipótesis en lo que respecta, en este caso, a su carácter térmico y fluidomecánico. Es por ello que la coherencia en la declaración de premisas a la hora de preceder a un modelo de estas características es importante. Por otro lado, cabe mencionar que cualquier simplificación que se considere se verá traducido, en mayor o menor medida, en una reducción de la carga computacional del programa, o, dicho con otras palabras, en una reducción del tiempo de simulación empleado. Previamente a listar estas consideraciones conviene declarar la idea de los modelos que se generarán a fin de entender mejor el por qué de estas hipótesis. Se desarrollarán principalmente dos modelos sobre el mismo edificio: • Un modelo simple, donde las temperaturas a calcular e involucradas en el cálculo serán únicas para el volumen que representan. Una sola temperatura (media) gobierna el estado de un elemento a considerar. • Un modelo 3-D, donde se aplica un mallado al volumen de aire, y su correspondiente extrapolación a las envolventes que lo contienen. El primer modelo se desarrolla a modo de ejemplificar el comportamiento del sistema y poder tener una referencia, ya que no se dispone de una trayectoria de referencia o la posibilidad de hacer mediciones sobre un modelo físico. De este modo, se utilizará el modelo simple principalmente para validar las ecuaciones de transferencia de calor que gobierna el modelo, pudiendo hacer comprobaciones matemáticas de forma rápida en comparación al otro modelo. 3.1. Caracterización del aire interior En primera instancia, se han considerado todos los volúmenes de material existentes en el modelo como medios isotrópicos y de propiedades físicas constantes. Términos como el poder calorífico Cp, la densidad ρ, la conductividad térmica k o el coeficiente de película convectiva h no se verán afectados por la temperatura que se esté evaluando en el instante de la simulación. En primera instancia, serán tomados, por ejemplo en el aire, las propiedades físicas a temperatura ambiente para cualquier estado. Es una consideración que es entendible en términos del simulador que se plantea. Las temperaturas exterior e interior que se esperan manejar en estos modelos no excederán en gran medida tanto por exceso como por defecto estos valores, lo que no influirá significativamente en el resultado previsiblemente. Asimismo, una consideración importante en este tema ha sido no establecer flujos macroscópicos del aire interior (o corrientes). Es una decisión que será justificada conforme se exponga el modelo 3-D, y que será discutida y corregida con posterioridad. 8 3.2. Mecanismos de transferencia de calor Por teoría de transferencia de calor conocemos los tres principales mecanismos de transferencia que existen en la materia: radiación, convección y conduccón (o transferencia de masa). Las siguientes definiciones y expresiones matemáticas están representadas en régimen permanente a fin de facilitar la comprensión de las mismas. No obstante, en lo que prosigue a este apartado, se utilizarán términos diferenciales a fin de modelar el transitorio. • Conducción. Es el principal mecanismo en el gradiente existente en las paredes, y, por otro lado, en el gradiente de temperatura del volumen de aire interior en el modelo 3-D que se plantea. Consiste en la transferencia de calor por diferencia energética en partículas adyacentes en un mismo medio. 𝑄󰇗𝑐𝑜𝑛𝑑𝑢𝑐𝑐𝑖ó𝑛=𝑘𝐴(𝑇2−𝑇1) 𝐿 Donde k es el coeficiente de transferencia de calor en W/(mK), A es el área transferencia de calor, en m2. Tn es temperatura del nodo n en K ó ºC (nótese que, al ser una diferencia, es posible utilizar cualquiera de las dos unidades para este cálculo). L la distancia medible que atraviesa el flujo de calor entre los dos puntos evaluados, y finalmente 𝑄󰇗𝑐𝑜𝑛𝑑𝑢𝑐𝑐𝑖ó𝑛 es la cantidad de calor en W transmitida. • Convección. Es el mecanismo de transferencia de calor presente entre los medios sólidos de las envolventes (paredes, techo y suelo) y el aire interior y exterior. Este mecanismo implica efectos combinados entre conducción y cinemática del fluido en cuestión. De forma general podemos expresarlo como: 𝑄󰇗𝑐𝑜𝑛𝑣𝑒𝑐𝑐𝑖ó𝑛=ℎ𝐴(𝑇2−𝑇1) Donde h es el coeficiente de película convectiva en W/(m2K). A es el área compartida entre el sólido y el fluido, o el área de transferencia de calor, en m2. Tn es temperatura del cuerpo n en K ó ºC, y finalmente 𝑄󰇗𝑐𝑜𝑛𝑣𝑒𝑐𝑐𝑖ó𝑛es la cantidad de calor en W transmitida. • Radiación. La radiación es la energía emitida por la materia en forma de ondas electromagnéticas (o fotones) como resultado de los cambios en las configuraciones electrónicas de los átomos o moléculas. El presente trabajo desarrolla modelos puramente conductivo-convectivos, siendo la radiación omitida en la obtención del transitorio de la temperatura. En primera instancia, no implementar la radiación supone ahorrar tiempo computacional considerablemente, debido a las no linearidades y exponentes de cuarto orden que suponen las ecuaciones que rigen este mecanismo. Asimismo, la obtención de los factores de forma en espacios residenciales con diversos elementos (muebles, calefactores, tabiques, etc.) no es factible a la hora de hacer una aplicación de estas características. El autor Y. Cengel (2002) cita al respecto: “La radiación suele ser importante en relación con la conducción o la convección natural pero insignificante en relación con la convección forzada. Por lo tanto, la radiación en aplicaciones de convección forzada no es tenida en cuenta, especialmente cuando las superficies implicadas tienen bajas emisividades y temperaturas bajas o moderadas” (p.29). Los materiales típicos de las paredes de los edificios, como pueden ser el yeso, el mortero o el ladrillo, suelen tener emisividades en torno al 0.8, lo que se considera relevante. Además, la convección planteada en los modelos es natural, y los sistemas de calefacción implementados toman temperaturas altas. Todo ello hace pensar que estamos cometiendo un error considerable al no tener en cuenta esta radiación. Por otra parte, las temperaturas moderadas que encontramos en el aire interior y en las envolventes, sumado a sus lentas inercias térmicas, contrarrestan estos puntos. 9 Mendes et al. (2001), por su parte, analizaron la contribución de la radiación implementada en su modelo en términos de consumo energético e inercia térmica del sistema, concluyendo que los efectos de la radiación suponían pérdidas de energía pequeñas, representadas en retrasos en la evolución térmica del aire interior, y en un aumento de energía consumida en torno al 2% con respecto al radiador puramente convectivo. Todo ello hace pensar en lo despreciable de la radiación en el sistema que se pretende analizar, y el error cometido teniendo en cuenta esta hipótesis. Como forma alternativa e implementable en revisiones futuras, se propone implementar la transferencia de calor por radiación térmica modificando el coeficiente de película convectiva de los materiales, lo cual es un enfoque recurrente y validado muchos autores como el propio Cengel. La radiación de una superficie envuelta por un gas como en este caso, ocurre de forma paralela a la convección, por lo que es posible linealizar la expresión de la radiación desarrollando un coeficiente convectivo-radiante que cuantifique ambas transferencias de calor (Cengel, 2002). 16 12(𝜌1𝑒1𝐶1)𝜕𝑇1 𝜕𝑡=ℎ𝑎𝑖𝑟𝑒(𝑇𝐵𝑖−𝑇1𝑖)+𝐾1(𝑇2𝑖−𝑇1𝑖) 𝑒1 12(𝜌1𝑒1𝐶1)𝜕𝑇2 𝜕𝑡=ℎ𝑎𝑖𝑟𝑒(𝑇𝐶𝑖−𝑇2𝑖)+𝐾1(𝑇1𝑖−𝑇2𝑖) 𝑒1 5.1.3 Radiadores La variación de la temperatura en el radiador se ha obtenido de manera simplificada, tomando como volumen de control la masa de aceite que interactúa en el dispositivo. De esta manera, se descartan los transitorios en el armazón metálico y suponemos un aprovechamiento completo de la potencia neta del radiador, quedando de la forma: 𝜌𝑎𝑐𝑒𝑖𝑡𝑒𝑉𝑎𝑐𝑒𝑖𝑡𝑒𝐶𝑎𝑐𝑒𝑖𝑡𝑒𝜕𝑇𝑅𝑎𝑑𝐵 𝜕𝑡 =ℎ𝑎𝑖𝑟𝑒𝐴𝑅𝑎𝑑𝐵(𝑇𝐵𝑖−𝑇𝑅𝑎𝑑𝐵 𝑖)+𝑄󰇗𝑅𝑎𝑑𝐵 Donde 𝐴𝑅𝑎𝑑𝐵 es el área convectiva de transferencia de calor del radiador con el aire que lo envuelve, 𝑇𝑅𝑎𝑑𝐵 𝑖es la temperatura del volumen de aceite del radiador B en el instante i, y 𝑄󰇗𝑅𝑎𝑑𝐵 es la potencia consumida por el radiador, en este caso desde la red eléctrica, en W. 5.3. Método de derivación Para este modelo simple, donde el número de ecuaciones no supone un problema y se presenta como un comprobante de la funcionalidad del modelo matemático, se ha optado por implementar la función ode15s propia del software Matlab, que garantiza resultados más fiables en la derivación numérica, pudiendo configurar tiempos de paso más estrechos con tolerancias más ajustadas. 5.4. Implementación en Matlab En lo referente a este aspecto, la secuencia de procesos del simulador generado consiste en los siguientes pasos ordenadamente: 1. Definición de variables por el usuario. La primera parte del código consiste en la definición del tiempo de simulación total u horizonte de predicción, las variables de entrada o set-point, y la dimensionalización de las variables involucradas en el código, que comprende la reserva de memoria y declaración de las matrices de temperatura para su futura utilización. 2. Bucle de cálculo matemático. Se calcularán aquí los incrementos de temperatura que se irán sumando a cada intervalo de tiempo predefinido, hasta llegar al horizonte de control. Para ello, el programa sobreescribe los valores de las temperaturas para hacer estos cálculos. En cada iteración, se declaran variables auxiliares que guardan estos valores para tener un registro de la evolución de las variables. 3. Guardado automático de matrices de temperatura al instante final y del historial de temperaturas medias a lo largo de la simulación. 17 6 MODELO 3-D En este capítulo se aplica una discretización en el espacio para el cálculo del transitorio de las temperaturas en el edificio del modelo. Los fundamentos matemáticos y las ecuaciones que gobiernan el comportamiento de este modelo son similares al modelo univariable, salvo que, para este caso, se implementa la cinética de transferencia de masa en el aire. Por otra parte, se han acomodado los términos matemáticos implícitos en cada iteración del bucle para reducir tiempos computacionales, y se han desarrollado algoritmos adicionales para conectar matrices de temperatura y realizar bucles de posicionamiento. 6.1. Discretización del espacio 6.1.1 Aire interior En este modelo se utilizará como unidad elemental un “cubo” de aire. Las dimensiones de este cubo están declaradas previa ejecución de la simulación por el usuario por las variables dx, dy y dz, que definen las dimensiones (en m) del cubo infinitesimal en sus tres ejes cartesianos. De esta forma, si la distancia de un habitáculo en el eje X midiera 3 metros y fijáramos un cubo de dimensiones 0.05 m de lado, el número de nodos que cabrían en esta dimensión quedaría definido por el número 𝑁𝑜𝑑𝑜𝑠 𝑋=𝐷𝑖𝑠𝑡𝑎𝑛𝑐𝑖𝑎/𝑑𝑥−1=3/0.05−1=59 Y de igual forma para las restantes dos dimensiones, quedando el volumen de aire definido por 𝑁𝑜𝑑𝑜𝑠 𝑋∗𝑁𝑜𝑑𝑜𝑠 𝑌∗ 𝑁𝑜𝑑𝑜𝑠 𝑍=𝑛º 𝑑𝑒 𝑛𝑜𝑑𝑜𝑠 𝑑𝑒𝑙 𝑣𝑜𝑙𝑢𝑚𝑒𝑛 De esta forma queda discretizado el volumen completo de aire interior de cada compartimento del edificio según unas dimensiones infinitesimales dadas, dando como resultado un conjunto de nodos o volúmenes diferenciales que componen en su conjunto el mismo volumen contemplado en el modelo simple. Por consiguiente, a cada extremo de la dimensión discretizada quedaría sin contemplar una distancia igual a dx/2, las cuales colindan con las paredes del edificio. En estos puntos, por otra parte, aplicaremos las pertinentes condiciones de contorno aire-pared, donde dicha diferencia en un mallado que pretende ser lo más infinitesimal posible no conllevaría problemas de continuidad, siendo este espacio representado a una temperatura igual a la de la última capa de la pared colindante. 6.1.2 Envolventes Para las paredes del edificio, se han discretizado las dimensiones del plano de la pared de igual forma que el volumen de aire, quedando iguales áreas de sección transversales que tendrían los cubos de aire colindantes. Esto implica una correcta transmisión y correspondencia entre variables termodinámicas, quedando una distribución de temperaturas a lo largo de cada pared. Estando una pared contenida en el eje x del espacio cartesiano, quedará subdividida en su superficie en una matriz de nodos de dimensiones de superficie dx·dz. Asimismo, este conjunto de nodos se repetirá de igual forma por cada capa de distinto material que componga el elemento pared. En este caso, todas las paredes exteriores están compuestas por 3 capas, cuyo espesor de capa varía en cada una de ellas. De esta forma, las paredes quedarían compuestas por n nodos de superficie en cada capa con espesor (o longitud transversal) ex. 18 6.1.3 Radiadores Para la caracterización de los radiadores en el modelo en 3D, se ha procedido a modelar el aporte de la transferencia de calor del radiador, como un conjunto de planchas verticales sin entidad física. Pese a la complejidad de generar sólidos en un espacio mallado y robusto, esta aproximación genera el resultado deseado en términos energéticos, si bien no hemos contemplado el efecto de la radiación en otros posibles elementos sólidos que pudiesen existir en el interior del edificio. 6.1.4 Ventanas De forma similar a las paredes, las ventanas se han discretizado siguiendo la sección de la malla del aire, para facilitar las ecuaciones de transmisión térmica en el simulador. Si bien las ventanas no se han modelado como un volumen con inercia térmica, si no como una resistencia térmica al paso de energía calorífica entre las temperaturas interior y exterior, el efecto del mallado en este modelo no supone algún cambio con respecto al modelo simple. 6.2. Código: definición de variables La principal diferencia con respecto al modelo simple, es la necesidad de un trabajo previo adicional en la dimensionalización de las matrices que van actuar en el programa. Para ello, es necesario dimensionalizar una serie de variables previas que comprenden. • Dimensiones físicas del edificio y distancias infinitesimales dx, dy, dz. • Matrices de temperaturas en las habitaciones A, B y C. Cada elemento i, j, k de la matriz representará un cubo diferencial de posicionamiento único en el espacio. • Matrices de temperaturas en las paredes, cuyos dos primeros índices determinarán la posición del elemento diferencial en el plano de la pared, y el tercer índice determinará el límite de capa que se esté considerando. • Matrices de posicionamiento de radiadores. Se trata de una matriz tridimensional donde se le dará, de forma manual por el usuario, un valor 1 a modo de bandera que identifique que en la posición i, j, k existe un radiador emplazado físicamente. Las dimensiones serán las mismas que las matrices A, B y C, aumentadas cada subíncide en 2. Por ejemplo, si la matriz A quedara definida inicialmente como una matriz A=zeros(20,20,20), la matriz de posicionamiento de radiadores para esta habitación quedaría de la forma RadA=zeros(20+2, 20+2, 20+2). Esto se debe a una decisión en la ordenación y tratamiento de datos que se explicará más adelante, y que simplifica la estructura del código. • Matrices de posicionamiento de ventanas. Matrices de iguales características que las matrices de posicionamiento de radiadores. • Matrices auxiliares de procesamiento. Estas matrices servirán para ordenar las submatrices de temperaturas en la ejecución matemática, a fin de simplificar el proceso y reducir tiempos computacionales. Sirven para poder comparar en una misma matriz todo el conjunto de elementos que intervienen en el análisis térmico, donde están incluidas las anteriores matrices de posicionamiento y de temperaturas paredes y aire. 19 6.3. Código: definición y funciones En este apartado se pretende declarar, de manera intuitiva, el recorrido de la información a través del código, definiendo las funciones principales de cada función que se ha creado para el procesado de los datos en el mismo. El modelo 3-D se compone de las siguientes etapas principales: 1. Definición de variables y condiciones iniciales. 2. Bucle: obtención de matrices de posicionamiento, derivación numérica y obtención de Temperaturas 𝑇𝑡+𝑑𝑡 3. Generación de banco de resultados. 6.2.1 Definición de variables Las variables a definir serán principalmente las nombradas en el apartado 6.2, más una serie de variables auxiliares que serán necesarias declarar para ejecutar el código, como el horizonte de predicción deseado N y el intervalo de tiempo diferencial dt, datos que marcarán el número de iteraciones en el bucle. 6.2.2 Bucle En el instante t=0, el bucle recibe el conjunto de matrices de temperaturas en las condiciones iniciales que ha marcado el usuario. Este paquete comprende las matrices de temperaturas de aire interior, envolventes, las matrices de posicionamiento de ventanas y radiadores, y, por último, la variable de temperatura exterior. Tomando la matriz de temperaturas de la habitación A como ejemplo, el siguiente paso será generar la matriz de posicionamiento de temperaturas para este volumen de aire encerrado en el cuarto. Para ello, se ha diseñado la función “GeneraMatriz”: MatrizA=GeneraMatriz(ANx,ANy,Nz,A,Toldpared1A,Toldpared4A,ToldparedAB,Toldpar edAC,ToldtechoA,Tsuelo); Siendo A la matriz de temperaturas en el instante t, ANx, ANy, y Nz las dimensiones de la matriz A, y el resto de variables las matrices de temperaturas en el instante t de las envolventes que envuelven el espacio A. La función devuelve por salida una matriz “MatrizA” de dimensiones (ANx+2, ANy+2, Nz+2). Nótese que estos dos nuevos elementos en cada dimensión corresponden a las temperaturas interiores de las envolventes que rodean la matriz A por cada lado en cada eje cartesiano. Una vez generada la matriz de posicionamiento, se introducirá en la función “GradienteMod”: A=GradienteMod(ANx,ANy,Nz,dx,dy,dz,Aold,MatrizA,VentanaA,SensorVentanaA,Kwin, ewin,hi,hext,rho,Cp,dt,Kaire,Text,RadA,TeqA,TmedA,1); Siendo VentanaA y RadA las matrices de posicionamiento de las ventanas y radiadores de la habitación A, SensorVentanaA una variable booleana que identifique si la ventana está abierta o cerrada, Kwin la resistencia térmica equivalente de la ventana (en W/mK), ewin el espesor total de la ventana (en m), hi y hext los coeficientes de película interno y externo respectivamente, Cp la capacidad calorífica del airea presión constante (en J/KgK), Kaire la conductividad térmica (W/m2K), Text la temperatura exterior, TeqA la temperatura del radiador A, TmedA la temperatura media del volumen de aire A, y el último número corresponde a una variable auxiliar que define el tipo de malla a utilizar entre dos modelos prestablecidos. En esta función se obtendrá mediante el método de derivación de las diferencias finitas las temperaturas en el instante t+dt. Para ello, en cada elemento nodal, el programa calculará para cada lado del cubo la transferencia 20 de calor asociada al elemento colindante en cada dirección de los ejes cartesianos. Se ha dispuesto el cubo como una variable struct, donde cada lado viene identificado con una variable única predefinida de la a a la f, tal y como se muestra en la figura. Figura 6-1. Esquema de evaluación posicional para cubo diferencial de aire. Suponiendo que en el elemento 𝑀𝑎𝑡𝑟𝑖𝑧𝐴𝑖,𝑗,𝑘 tengamos: • Lado a: elemento colindante ventana. 𝑎=(𝑇𝑒𝑥𝑡−𝑀𝑎𝑡𝑟𝑖𝑧𝐴𝑖,𝑗,𝑘 2 ℎ𝑎𝑖𝑟𝑒+𝑒𝑤𝑖𝑛 𝑘𝑤𝑖𝑛 )/𝑑𝑦 • Lado b: elemento colindante radiador. 𝑏=ℎ𝑎𝑖𝑟𝑒(𝑇𝑅𝑎𝑑𝐴−𝑀𝑎𝑡𝑟𝑖𝑧𝐴𝑖,𝑗,𝑘)/𝑑𝑥 • Lado c: elemento colindante pared. 𝑐=ℎ𝑎𝑖𝑟𝑒(𝑀𝑎𝑡𝑟𝑖𝑧𝐴𝑖−1,𝑗,𝑘−𝑀𝑎𝑡𝑟𝑖𝑧𝐴𝑖,𝑗,𝑘)/𝑑𝑦 • Lados d, e, f: elemento colindante aire 𝑑=𝐾𝑎𝑖𝑟𝑒(𝑀𝑎𝑡𝑟𝑖𝑧𝐴𝑖,𝑗−1,𝑘−𝑀𝑎𝑡𝑟𝑖𝑧𝐴𝑖,𝑗,𝑘)/𝑑𝑥2 𝑒=𝐾𝑎𝑖𝑟𝑒(𝑀𝑎𝑡𝑟𝑖𝑧𝐴𝑖,𝑗,𝑘−1−𝑀𝑎𝑡𝑟𝑖𝑧𝐴𝑖,𝑗,𝑘)/𝑑𝑧2 𝑓=𝐾𝑎𝑖𝑟𝑒(𝑀𝑎𝑡𝑟𝑖𝑧𝐴𝑖,𝑗,𝑘+1−𝑀𝑎𝑡𝑟𝑖𝑧𝐴𝑖,𝑗,𝑘)/𝑑𝑧2 21 Obtendremos finalmente la predicción de la temperatura 𝑀𝑎𝑡𝑟𝑖𝑧𝐴𝑖,𝑗,𝑘 𝑡+𝑑𝑡 aplicando: 𝐴𝑖−1,𝑗−1,𝑘−1 𝑡+𝑑𝑡 =𝑑𝑡 𝜌𝐶𝑝∗(𝑎+𝑏+𝑐+𝑑+𝑒+𝑓)−𝐴𝑖−1,𝑗−1,𝑘−1 𝑡 Nótese que se han acomodado los términos de áreas y volúmenes para poder ejecutar el conjunto de funciones. Asimismo, cabe mencionar que los elementos i, j, k de MatrizA correspondería al elemento i-1, j-1, k-1 de la matriz de temperaturas A, deshaciendo de esta forma el proceso de posicionamiento anteriormente citado. De esta forma, se obtiene la predicción de la temperatura en t+dt, por el método de las diferencias finitas, dependiendo del instante predecesor t, y de las aportaciones caloríficas representadas anteriormente. Tras la obtención de las temperaturas en el aire interior, el siguiente paso en el programa reside en la obtención de las temperaturas correspondientes a las distintas capas de las envolventes. El método de derivación seguirá siendo el de las diferencias finitas debido a la cuantía de operaciones en cada iteración del programa. Para tal fin se ha diseado la función EnvolventeNewton, que implementa las ecuaciones matemáticas que se estipulaban en el 5.1.2., siendo necesarios algunos datos adicionales como la orientación de la pared a considerar con respecto al volumen de aire que contiene (identificador 1, 2, 3 ó 4). [Tpared1A]=EnvolventeNewton(ANx,Nz,dt,e,d,c,K,hi,hext,Aold,Toldpared1A,Text,1 ); [Tpared4A]=EnvolventeNewton(ANy,Nz,dt,e,d,c,K,hi,hext,Aold,Toldpared4A,Text,4 ); [TparedAB]=Envolvente2Newton(ANy,Nz,dt,e,d,c,K,hi,Aold,Bold,ToldparedAB,3); [TparedAC]=Envolvente2Newton(ANx,Nz,dt,e,d,c,K,hi,Aold,Cold,ToldparedAC,2); Con estas 4 llamadas a la función se completaría el cálculo de las temperaturas en las paredes. Estas 4 sentencias corresponden a las envolventes situadas en cada orientación correspondiente según los identificadores N, S, E, y O que hemos prestablecido con números del 1 al 4. Nótese que la función Envolvente2Newton es la correspondiente a las paredes interiores, compuestas de dos elementos nodales por cada posición que recorre la matriz del plano, según las ecuaciones estipuladas en la sección 5.1.2. Tomando como ejemplo la pared en la orientación Norte (identificador 1) que colinda con el aire interior del cuarto A, la obtención matemática de las envolventes quedaría de la forma: 𝑇𝑗,𝑘,1 𝑡+𝑑𝑡=𝜕𝑡 0.5𝜌1𝑒1𝐶1(ℎ𝑎𝑖𝑟𝑒(𝐴1,𝑗,𝑘 𝑡−𝑇𝑗,𝑘,1 𝑡)+𝐾1(𝑇𝑗,𝑘,2 𝑡−𝑇𝑗,𝑘,1 𝑡) 𝑒1)+𝑇𝑗,𝑘,1 𝑡 𝑇𝑗,𝑘,2 𝑡+𝑑𝑡=𝜕𝑡 0.5(𝜌1𝑒1𝐶1+𝜌2𝑒2𝐶2)(𝐾1(𝑇𝑗,𝑘,1 𝑡−𝑇𝑗,𝑘,2 𝑡) 𝑒1+𝐾2(𝑇𝑗,𝑘,3 𝑡−𝑇𝑗,𝑘,2 𝑡) 𝑒2)+𝑇𝑗,𝑘,2 𝑡 𝑇𝑗,𝑘,3 𝑡+𝑑𝑡=𝜕𝑡 0.5(𝜌2𝑒2𝐶2+𝜌3𝑒3𝐶3)(𝐾2(𝑇𝑗,𝑘,2 𝑡−𝑇𝑗,𝑘,3 𝑡) 𝑒2+𝐾3(𝑇𝑗,𝑘,4 𝑡−𝑇𝑗,𝑘,3 𝑡) 𝑒3)+𝑇𝑗,𝑘,3 𝑡 𝑇𝑗,𝑘,4 𝑡+𝑑𝑡=𝜕𝑡 0.5𝜌3𝑒3𝐶3(ℎ𝑎𝑖𝑟𝑒(𝑇𝑒𝑥𝑡 𝑡−𝑇𝑗,𝑘,4 𝑡)+𝐾3(𝑇𝑗,𝑘,3 𝑡−𝑇𝑗,𝑘,4 𝑡) 𝑒3)+𝑇𝑗,𝑘,4 𝑡 Nótese que, en la pared norte, la variable que recorre el plano de la pared será el subíndice j de la matriz de temperaturas A, que junto con la altura expresada por el subíndice k solucionaría la correspondencia entre ambas matrices. De igual forma, para la pared interior que separa las habitaciones A y B, con orientación Este, tendríamos las ecuaciones: 22 𝑇𝑖,𝑘,1 𝑡+𝑑𝑡=𝜕𝑡 0.5𝜌𝑒𝐶(ℎ𝑎𝑖𝑟𝑒(𝐴𝑖,𝑒𝑛𝑑,𝑘 𝑡−𝑇𝑖,𝑘,1 𝑡)+𝐾(𝑇𝑖,𝑘,2 𝑡−𝑇𝑖,𝑘,1 𝑡) 𝑒)+𝑇𝑖,𝑘,1 𝑡 𝑇1,𝑘,2 𝑡+𝑑𝑡=𝜕𝑡 0.5𝜌𝑒𝐶(ℎ𝑎𝑖𝑟𝑒(𝐵𝑖,1,𝑘 𝑡−𝑇𝑗,𝑘,1 𝑡)+𝐾1(𝑇𝑖,𝑘,1 𝑡−𝑇𝑗,𝑘,2 𝑡) 𝑒)+𝑇𝑖,𝑘,2 𝑡 En este caso, el subíndice que recorre la pared sería el i para las matrices de temperaturas de los aires interiores A y B. Una vez calculado el paquete de temperaturas predictivas en el instante t+dt, el programa actualiza asimismo la temperatura exterior según la función diseñada en el apartado 4.4. Una vez realizados los nuevos cálculos, el bucle se posicionará en el siguiente instante de tiempo, reiniciíandose el conjunto de operaciones con la siguiente condición inicial: 𝑇𝑡=𝑇𝑡+𝑑𝑡 23 7 SIMULACIÓN Y CONCLUSIONES En el presente capítulo se analizarán los modelos ejecutados bajo un mismo conjunto de condiciones iniciales y horizontes de predicción, a fin de conocer y analizar el comportamiento del programa bajo un conjunto de premisas. 7.1. Simulaciones El estudio comprenderá tres simulaciones distintas: • Modelo simple (monovariable), mediante el solucionador de ecuaciones diferenciales ordinarias ode15s. • Modelo mallado, con dimensión infinitesimal igual a 0.1m. • Modelo mallado, con dimensión infinitesimal igual a 0.05m. 7.2. Condiciones iniciales El conjunto de condiciones iniciales comunes se resume en: • Temperatura inicial del aire interior en las tres habitaciones: 13ºC. • Temperatura inicial del aceite de los tres radiadores: 13ºC. • Temperatura inicial de todas las capas en las envolventes: 11ºC. • Temperatura del suelo (constante, superficie adiabática): 13ºC. • Hora de incialicación del programa a las 00:00 horas, lo que corresponde a una temperatura exterior de 8ºC. • Horizonte de predicción: 6 horas. Temperatura exterior en el instante final: 5.59ºC. • Implementación de un controlador sencillo on-off, cuya variable de entrada es la temperatura media del aire interior de cada habitación, y la salida será una variable booleana que apague cada radiador, una vez superados los 25ºC de media, y lo encienda una vez descendida la temperatura por debajo de los 18ºC. 24 7.3. Resultados En las siguientes tres figuras se muestran las evoluciones de las temperaturas medias para las habitaciones A y C (A y B muestran comportamientos similares en esta simulación) según los tres modelos empleados. Figura 7-1. Evolución de las temperaturas medias. El resultado obtenido en el modelo simple difiere sustancialmente de los modelos tridimensionales. La dinámica de este modelo produce unas evoluciones muy exageradas en lo que a temperatura se refiere, alcanzándose, desde la temperatura inicial de 13ºC, los 25ºC en apenas 1,5 horas. Se trata de un sistema donde se evalúan constantemente en la ecuación de transmisión de calor la temperatura de la fuente de calor, en torno a 150 ºC, y la temperatura de la masa de aire monovariable, inicialmente en 13ºC. Esto produce un aporte calorífico prácticamente constante hasta que la masa total de aire alcance el nivel objetivo en el controlador. Por otra parte, en los modelos tridimensionales se obtienen unas dinámicas más suaves, donde la primera parte del transitorio, donde se inician los radidadores desde una temperatura igual a la del aire interior, tiene una 25 pendiente más acentuada. Esto se debe a que la diferencia de temperaturas entre la pared del radiador y los nodos colindantes es más dispar, produciéndose por ende una mayor transferencia de calor atendiendo a la ecuación tradicional de transferencia de calor por convección. Conforme los nodos de aire anexos al radiador toman temperatura, el factor que predomina en la dinámica del sistema es la capacidad de transferencia de calor del aire impuesta, habiendo modelado a tal efecto la masa de aire como un “sólido” compuesto por elementos nodales de igual volumen. De esta forma, la ecuación de difusión se aproxima a la transferencia de masa, quedando como variables predominantes en la velocidad del sistema las propiedades físicas del aire. En primera instancia también se puede observar la diferencia de comportamiento entre los distintos espacios A y C. Pese a tener ambos la misma superficie útil, la distribicuión del espacio no es similar. Atendiendo a las dimensiones en planta, la relación entre longitudes en el espacio A es de 4:3, mientras que la relación en el espacio C es de 2. Si espacialmente el radiador se encuentra en el centro de mayor dimensión, es evidente que la configuración alargada de la habitación C será menos eficiente a la hora de transmitir potencia calorífica por difusión a los extremos del espacio. En la siguiente imagen se muestra cómo se refleja el gradiente de temperaturas en el simulador dependiendo del espacio, sobre la habitación C en este caso, a modo de ejemplo. Figura 7-2. Gradiente de temperaturas habitación C.