Computacional program for the 3D simulation
of the indoor temperatura distribution 

José A. Dopazo, Nastia Almao y José Rincón 

Laboratorio de Simulación Computacional, Departamento de Energía,
Escuela de Ingeniería Mecánica, Facultad de Ingeniería. Universidad del Zulia.
Maracaibo, Venezuela. E-mail: jdopazo@luz.ve 

Abstract 

This work presents the development of a computational program called EVITA 3D written in Fortran 77, for the simulation of the indoor temperature distribution. The program solves the unsteady form of continuity, momentum and energy equations using a Piso-like scheme on a non-staggered grid and with a bounded high order treatment for the convection terms. The 3-D physical domain is divided into several zones to include the different materials of the building (e.g. walls, floor, ceiling). The energy equation is solved using transient boundary conditions such as outdoor temperature, horizontal and vertical surfaces global solar radiation, long wave radiative exchange between external surfaces and sky, floor temperature and mean wind velocity. The movement of air inside the building is simulated using the Boussinesq approach. EVITA 3D results were compared with experimental data, showing good agreement. 

Key words: 

Energy simulation, Indoor temperature distribution, CFD. 

Algoritmo computarizado para la simulación tridimensional de la distribución de temperaturas en el interior de una vivienda 

Resumen 

Este trabajo presenta el desarrollo del código EVITA 3D, elaborado con el propósito de simular la distribución de temperaturas en el interior de una vivienda. Para su desarrollo, se utilizó el método de los volúmenes finitos en mallas no desplazadas y tratamiento acotado de alto orden de los términos convectivos. Se adoptó un esquema similar al algoritmo PISO para la solución de las ecuaciones acopladas de Continuidad, Cantidad de Movimiento y Energía, bajo condiciones transitorias. El dominio de cálculo se modela en un ambiente cartesiano tridimensional divido en zonas, sobre las que se construye la malla de manera que se pueda considerar cualquier material que forme parte de la edificación a evaluar. Se establecen condiciones de contorno transitorias para la solución de la ecuación de Energía tales como: la temperatura ambiente, radiación solar global sobre superficies exteriores, intercambio radiativo de onda larga entre las superficies externas y el cielo, temperatura del suelo y la velocidad promedio del viento. El movimiento del aire en el interior de la edificación se modela utilizando la aproximación de Boussinesq. Para evaluar el desempeño del modelo desarrollado fueron realizadas diferentes simulaciones de casos con información conocida experimentalmente, obteniéndose resultados satisfactorios, lo que permite considerar el código una herramienta útil. 

Palabras clave: 

Simulación computacional, distribución de temperaturas, comportamiento térmico. 

Recibido el 05 de Marzo de 2003

En forma revisada el 22 de Marzo de 2004 

Introducción 

El tema de optimización energética es, hoy por hoy, materia primordial en cualquier comunidad a nivel mundial. Uno de los focos de consumo identificado lo constituye el consumo de energía dedicada a lograr el acondicionamiento térmico de viviendas. Se hace necesaria la elaboración de herramientas que permitan realizar la evaluación del comportamiento térmico de edificaciones con el propósito de determinar la mejor manera de climatizarlas adecuadamente al menor costo. En esta dirección han sido desarrollados proyectos de investigación que incluyen, tanto experimentación, como elaboración de códigos computacionales que han permitido comprobar que es posible la reducción de la temperatura interior con la utilización de diferentes sistemas pasivos de enfriamiento. En este trabajo se presenta el desarrollo del código computacional EVITA 3D, el cual es una extensión del código EVITA [1]. Con el código desarrollado se puede simular la distribución de temperaturas en el interior de una vivienda y constituye un modelo tridimensional basado en el método de los volúmenes finitos, con mallas no desplazadas y tratamiento acotado de alto orden para el manejo de los términos convectivos [2], en donde el algoritmo utilizado para la solución en condiciones transitorias de las ecuaciones de Continuidad, Cantidad de Movimiento y Energía en forma acoplada, es similar al algoritmo PISO [3]. Para la validación del código desarrollado se hicieron simulaciones de diferentes casos comparando los resultados con datos experimentales registrados en módulos construidos a escala natural para la evaluación y caracterización de sistemas pasivos de enfriamiento [4]. Del análisis de los resultados se puede afirmar que el código desarrollado puede utilizarse como una herramienta útil y de alta confiabilidad para la obtención de la distribución de temperaturas en el interior de edificaciones construidas con diferentes materiales y sometidas a condiciones climáticas no permanentes. Para propósitos de investigación en relación a diferentes propuestas para “mejorar” las condiciones climáticas en el interior de edificaciones que podrían ser aplicables en Venezuela, resulta de mayor conveniencia tener un código computacional propio, lo cual permite adaptarlo a diferentes situaciones propuestas, en contraposición a códigos comerciales de usos restringidos. 

Detalles Numéricos
del Código EVITA 3D 

Las ecuaciones gobernantes de los fenómenos de transferencia de calor y flujo de fluidos, son las ecuaciones de conservación de la masa, cantidad de movimiento y energía, aplicadas a un fluido viscoso e incompresible, con propiedades físicas constantes, en proceso no isotérmico, y pueden ser expresadas por la ecuación diferencial general [5] 

(1) 

la cual representa cualquier principio de conservación en un volumen de control diferencial, en donde el primer término representa la variación de la variable de estudio (f) con respecto al tiempo (t), el segundo representa la variación del flujo convectivo, el tercero corresponde a la variación del flujo difusivo y el último sería un término llamado fuente (S) en el que se incluye cualquier factor que no se pueda incluir en los términos anteriores (r: es la densidad; G, representa el coeficiente difusivo; u, la velocidad) El método de solución se basa en la aproximación de los volúmenes finitos [5] en mallas no desplazadas, y las ecuaciones se reducen a un conjunto de ecuaciones algebraicas discretizadas, cuya forma general es 

                                                                          (2) 

donde “a” son los coeficientes y f la variable a evaluar en cada posición nodal de la malla. Los subíndices “_E,W,N,S,T,B ”, indican cada uno de los nodos vecinos posibles en un sistema cartesiano tridimensional. Los términos convectivos son discretizados utilizando un esquema de alto orden, desarrollado con asistencia del diagrama de Variables Normalizadas. Los términos difusivos son aproximados usando diferencias centradas. Para la solución del conjunto de ecuaciones algebraicas que resulta de la discretización de la ecuación diferencial general, se utiliza un método iterativo que incluye: el TDMA (Algoritmo de la matriz tridiagonal), combinado con un barrido línea por línea, un esquema de corrección en bloque [6], y un procedimiento completamente implícito de marcha en el tiempo. En la solución en forma acoplada de las ecuaciones de cantidad de movimiento y conservación de la masa, se implementó un algoritmo adicional de solución similar al algoritmo PISO [3]. 

Condiciones de contorno 

En la evaluación del comportamiento térmico de edificaciones y sistemas pasivos de enfriamiento, se tienen como condición de contorno: Flujos de calor en todos los bordes, excepto en el piso, donde se especifica el valor de la temperatura del suelo [1]. En paredes y techo se considera el flujo de calor debido a la radiación solar, flujo de calor por convección debido al viento, el intercambio por radiación de onda larga entre las superficies y el medio ambiente, y para el sistema pasivo considerado, se tiene adicionalmente sobre el techo, transferencia de calor por evaporación o condensación de agua e intercambio por radiación entre dos placas paralelas. Se considera un modelo simple de turbulencia para el flujo de aire interior y no se toma en cuenta el intercambio de calor por radiación entre las superficies interiores de paredes, piso y techo. Como solo se está considerando la carga térmica a través del piso y superficies exteriores, y el movimiento del aire interior por diferencias de densidades, las condiciones de contorno involucran: 

• Funciones dependientes del tiempo, de la temperatura ambiente e irradiancia global sobre superficies horizontales y verticales orientadas, cuyo modelado se realiza dependiendo de que existan o no datos promedio horarios medidos. 

En el caso de que existan datos promedios horarios, el programa interpola entre los valores dados para cada hora; de lo contrario, estas funciones se modelan de acuerdo a lo propuesto por Almao [7]. 

• Valores promedio de la velocidad del viento y valores promedio horario de la humedad relativa, correspondientes al mes o día durante el cual se realiza la simulación. 

• Propiedades ópticas y térmicas constantes de los materiales de construcción. 

• Temperatura del subsuelo. 

Flujos de calor en la superficie
del techo y paredes 

En ausencia de sistemas pasivos de enfriamiento, en la superficie del techo el flujo de calor viene dado por 

Q_TECHO = Q_CONV + Q_SOL + Q_RSC     (3) 

donde Q_CONV es la transferencia de calor por convección debido al viento, Q_SOL es el flujo de calor radiante que absorbe la superficie expuesta y Q_RSC es el intercambio radiativo superficie-cielo. Cada uno de estos flujos está determinado como sigue: 

Q_CONV = h_w (T_amb(t) – T_TECHO(t))     (4) 

Q_SOL = a I_TECHO         (5) 

Q_RSC = e_s s ((T_C⁴(t) – T_TECHO⁴ (t))     (6) 

donde h_w es el coeficiente convectivo debido al viento determinado según Clark [8], T_amb es la temperatura ambiente horaria, T_TECHO es la temperatura del borde del dominio en la superficie techo, I_TECHO es intensidad de radiación solar instantánea que incide sobre la superficie, a es la absortividad de la superficie, e_s es la emisividad de la superficie, s es la constante de Stefan-Boltzman y T_C es la temperatura promedio horaria del cielo que es función de la emisividad del cielo, que a su vez, depende de la temperatura de rocío del aire atmosférico [1]. 

Combinando las ecuaciones (4), (5) y (6) con (3) se tiene: 

Q_TECHO = h_w(T_amb(t) – T_TECHO(t)) + a I_TECHO +
         e_s s ((T_C⁴(t) – T_TECHO⁴(t))     (7) 

Esta ecuación se linealiza utilizando el método de la tangente a la curva de transferencia de calor vs. temperatura, recomendado por Patankar [6] adoptando la misma estructura del termino fuente: 

S = S_P f_P + S_C         (8) 

y en analogía con la ecuación (8) anterior, 

S_c = a I_TECHO + HSC         (9) 

S_p = – HTC         (10) 

donde 

HSC = e_s s(3T_TECHO⁴ – T_C⁴) + h_wT_amb     (11) 

HTC = h_w + 4e_s s T_TECHO^{* 3}     (12) 

representando (*) el valor de la variable calculado en la iteración anterior. 

De forma análoga, para las paredes los coeficientes respectivos se determinan mediante: 

S_c = a I_PARED + h_wT_amb + R_C     (13) 

S_p = – h_w + R_P         (14) 

donde 

R_C = e_PARED s (3T_PARED⁴ + T_C⁴ )/(1+e_PARED )    (15) 

R_P = – 4 e_PARED s T_PARED³/(1+e_PARED )    (16) 

Cuando se esté evaluando un sistema pasivo de enfriamiento, donde exista, además de enfriamiento por radiación nocturna, enfriamiento por evaporación de agua hacia la atmósfera, debe considerarse una intensidad de flujo de calor adicional Q_evap que viene dado por [1]: 

Q_evap = h_evap (P_s@tambHR-P_s@TB)   (17) 

donde 

h_evap=0.622 h_W h_fg@Tagua/(P_atm(C_P+ w C_Ps ))    (18) 

siendo h_fg@Tagua la entalpía de vaporización a la temperatura del agua en la superficie; C_P el calor específico del aire; ù la humedad específica y C_Ps el calor específico del vapor de agua. 

El valor de h_fg@Tagua se determinó mediante [1]: 

h_fg@Tagua= 3148.445 – 2.3696 T    (T en K)     (19) 

Este flujo de calor no es una función lineal de la temperatura en la superficie del agua, por lo que se hace necesario el mismo tratamiento dado al flujo de calor por radiación. La temperatura en la superficie del agua T_B se determina en función de la temperatura en el nodo vecino T_P al igualar el calor neto transferido a través de la superficie de agua con aquel transferido a través de los volúmenes de control adyacentes. Como el sistema pasivo de enfriamiento utilizado incluye control solar, la configuración del modelo difiere entre el período de asoleamiento y el nocturno: durante el día se tiene el panel aislante de poliestireno, el cual es removido durante la noche para permitir el intercambio radiativo nocturno y el enfriamiento evaporativo del agua del pozo. Con el fin de utilizar la misma malla en ambos períodos, durante el período nocturno son bloqueados los nodos correspondientes al aire sobre el pozo de agua y al panel aislante de poliestireno que tapa el contenedor de agua, esto se hace asignándole el valor de cero al coeficiente de difusión respectivo, lo que se traduce en un flujo de calor nulo. La intensidad de flujo de calor total sobre el techo Q_TECHO se suministra como un término fuente extra para los volúmenes de control adyacentes (agua) al contorno bloqueado mediante: 

Q_TECHO = Q_CONV + Q_RSC + Q_evap     (20) 

y sustituyendo cada uno de los flujos de calor por sus expresiones correspondientes: 

Q_TECHO = h_w(T_amb(t) – T_B(t)^*)+ e_s s((T_C⁴ (t) –
         T_B^{4 *} (t))  h_evap(P_S@TambHR – P_S@TB)
            (21) 

La ecuación resultante está en función de los valores de temperatura ambiente, temperatura de la superficie del agua y de la temperatura del cielo, los cuales son actualizados en cada paso de tiempo por ser función de éste. Para el período soleado, la condición de contorno en el techo es exactamente la misma que se utiliza en ausencia de sistemas pasivos de enfriamiento, por lo que se desactiva el bloqueo utilizado durante el período nocturno. 

Resultados y Análisis 

Para evaluar el desempeño del programa EVITA 3D, se desarrolló la simulación de la distribución de temperaturas en los módulos prototipos construidos en la Facultad de Arquitectura de L.U.Z. por González [4] de los cuales se tienen datos experimentales obtenidos por González [4], y simulados con otros programas obtenidos por Almao [1] y por Hernández [9]. 

La Figura 1 muestra detalles de los módulos prototipos (módulo de referencia y módulo experimental) construidos para la evaluación de sistemas pasivos de enfriamiento y de componentes constructivos [4]. 

Se puede apreciar que las edificaciones son idénticas a excepción de la construcción del techo en donde se ha establecido un sistema pasivo de enfriamiento en el módulo experimental. 

Para la simulación de cada módulo, se utiliza la misma malla con 31 nodos en el eje X, 31 nodos en el eje Y, y 27 nodos en el eje Z. Para la selección de esta malla, se tomaron como referencia las experiencias previas de Alamao [1]. La Figura 2 muestra la malla utilizada. 

La Figura 3 presenta las temperaturas del aire interior medidas y simuladas a 1.5 m de altura del piso y a 1.5 m de separación de cada pared del módulo de referencia, para los meses de Enero, Febrero, Marzo y Agosto. Se incluye el valor de la temperatura ambiente como referencia. 

Se puede apreciar que para las cuatro simulaciones realizadas, el comportamiento cualitativo es muy similar en comparación con las mediciones experimentales. El error relativo máximo fue de 1.75% correspondiente al mes de Enero (hora 2). Para ese mismo mes se obtuvo los máximos valores en el error porcentual absoluto, 1.3, y la desviación estándar de 1.4 (todos calculados con temperaturas en grados Celsius). En todos los días simulados las temperaturas máximas y mínimas simuladas se corresponden en horario a las medidas a excepción del mes de Marzo, en donde la temperatura mínima tiene un desfase de una hora (en adelanto), en relación al la medida. 

En las Figuras 4 y 5 se muestran los resultados de las temperaturas simuladas en las superficies exteriores e interiores del Módulo de Referencia para el mes de Febrero 1997. En la Figura 4 se puede apreciar que la mayor temperatura se alcanza en el techo cerca de la hora 14. Así mismo, se observa que la pared Este es la primera que registra incrementos de temperatura, la pared Oeste es la última que se caliente mientras que las paredes norte y sur registran un comportamiento similar de incremento de temperatura, siendo en todo momento, la temperatura de la pared sur superior a la de la pared norte; en ambos casos, la máxima temperatura se alcanza alrededor de la hora 14. Para la pared este, la máxima temperatura se alcanza cerca de la hora 12, mientras que para la pared oeste, la temperatura máxima se registra alrededor de la hora 16.