Fís. Sergio Ramos Herrera
e-mail:
Septiembre de 2022

Villahermosa, Tabasco, México

 

La ecuación de balance de masa

Para una compartimento simple, dónde al aire interior está bien mezclado, la ecuación de balance de masa para un contaminante interior establece lo siguiente:

\[Acumulación\ de\ masa = flujo\ de\ masa\ que\ entra - flujo\ de\ masa\ que\ sale + fuente-sumidero\]  

Cada término de la ecuación de balance tiene una forma matemática:

  1. La acumulación de masa es la derivada de la masa del contaminante en el interior con respecto al tiempo, \(dm/dt\). Si \(V\) representa el volumen del compartimento y \(C\) la concentración en un momento dado, entonces la masa \(m\) del contaminante en el aire interior se obtiene como \(C\cdot V\) y \(dm/dt=V\cdot dC/dt\) si el volumen \(V\) es constante.

  2. El flujo de masa que entra es igual al producto del flujo de aire que entra al compartimento \(Q\) por la concentración entrante \(Ce\), esto es, \(Q\cdot Ce\).

  3. El flujo de masa que sale es el producto de \(Q\) y la concentración \(C\) en el aire interior, esto es \(Q\cdot C\).

  4. El flujo volumétrico de aire que entra y sale del compartimento se puede estimar como \(Q=ACH\cdot V\), dónde ACH se conoce como la Tasa de Intercambio de Aire y tiene dimensiones de \(t^{-1}\).

  5. La fuente es cualquier proceso que emita al aire interior el contaminante de interés a una tasa o flujo de masa \(E\), con dimensiones de \(mt^{-1}\).

  6. Si el contaminante es conservativo y no existen filtros que eliminen el contaminante de interés, el sumidero es igual a cero.

 

Con estas definiciones, la ecuación queda como:

\[\begin{equation} \tag{ecuacion 1} V\frac{dC}{dt}=QCe +E-QC \end{equation}\]

 

La ecuación 1, es una ecuación diferencial lineal de primer orden. Su solución se puede consultar en cualquier libro de ecuaciones diferenciales ordinarias. Si se supone que la concentración inicial del contaminante en el interior del compartimento es Co, las soluciones son:

\[\begin{equation} \tag{ecuacion 2} C(t)=Co + \left(\frac{E}{V}\right)*t \end{equation}\]

\[\begin{equation} \tag{ecuacion 3} C(t)= \left(Ce + \frac{E}{ACH\cdot V}\right)\cdot \left[1-e^{-ACH\cdot t}\right] + Co\cdot e^{-ACH\cdot t} \end{equation}\]

 

La ecuación 2 es la solución cuando el compartimento esta cerrado, es decir, Q=0. La ecuación 3 es la solución cuando el compartimento esta abierto al exterior (puertas y/o ventanas abiertas) o esta cerrado pero tiene una tasa de infiltración de aire (por las rendijas en puertas y ventanas) significativa. En ambas ecuaciones, Co y Ce pueden estar en \(ppm\) o \(\mu g/m^3\) , V en \(m^3\), E en \(L/h\) o \(mg/h\), ACH en \(1/h\) y \(t\) en horas.

 

Cuando el “contaminante” es el dióxido de carbono, \(CO_2\), y el espacio está ocupado por N personas, la tasa de emisión se puede estimar como \(N\cdot E\), dónde E es ahora la tasa de generación de \(CO_2\) por persona, ya que el \(CO_2\) es un producto de la exhalación en el proceso de la respiración humana. En términos generales, se estima que un adulto realizando actividades de oficina genera 20 litros de \(CO_2\) por hora.

 

Figura 1. Esquema de entradas y salidas de un contaminante de interiores en un compartimento simple

 

Simulación con R

Es común que no se lleva un registro de cuantas personas hay en una hora determinada del día en un espacio interior. Si la condición térmica no es adecuada, la tendencia es abrir puertas y ventanas durante cierto tiempo o de manera permanente. También es común que no se conozcan el valor exacto de ciertos parámetros o que incluso estos varíen durante el periodo de estudio, como por ejemplo, el valor exacto de la tasa de generación de \(CO_2\) de cada ocupante. Esto es lo que llamaremos uso del espacio en condiciones reales. Para ganar idea del comportamiento del sistema en condiciones de uso real, se puede suponer que algunos parámetros varían de forma aletaoria siguiendo ciertro tipo de distribución, como la distribución normal.

 

El siguiente código emplea las ecuaciones 2 y 3 para calcular la concentración de \(CO_2\) en una oficina que puede estar ocupado como máximo por tres personas (\(N_\max =3\)) durante una jornada laboral de ocho horas, de 7 am a 3 pm. Debido a sus actividades, los ocupantes salen y entran de la oficina, de manera alternada o al mismo tiempo. El volumen de la oficina es de 173 \(m^3\) y a veces se mantiene ventilada naturalmente abriendo las puertas y ventanas y otras veces los ocupantes optan por mantenerlas cerradas.

 

#es una grafica que considera todos los casos posibles
con=matrix(0,nrow=1000, ncol=9);C=0
for (k in 1:9){
  for (i in 1:1000){
      cond=sample(c(0,1),1)#0=cerrado; 1=abierto
      N=sample(c(1,2,3),1)#no. personas dentro
      Ce=runif(1, min=380, max=400)#Ce varia entre 380 y 400 ppm
      if(k==1){Co=runif(1, min=364, max=400)}#varia entre 380 y 400 ppm
      ACH=runif(1, min=0.2, max=1)#ACH varia entre 0.2 y 1.0
      Q=ACH*173 #flujo en m3/s
      Resp=runif(1,min=12, max=16) #l/h
      t=seq(0, 1, length.out = 13)
      if (cond==0){
          C=Co+((N*Resp*1e6/173000)*t)}else{
             Cest=Ce+(Resp*1e06*N/(ACH*173000))
             C=Cest+((Co-Cest)*exp(-ACH*t))
      } 
      con[i,k]=median(C)
  }#fin del for i
  Co=median(con[,k])
}#fin del for k

 

Para simular el comportamiento del \(CO_2\) en la oficina en condiciones de uso real se programaron las ecuaciones uno y dos tomando en cuenta lo siguiente:

  1. Se supuso que los parámetros Ce, Co, ACH y E tienen una distribución uniforme (runif en R) dentro de cierto rango. Por ejemplo, Ce y Co pueden tomar cualquier valor distribuido uniformemente entre 380 ppm y 400 ppm.

  2. La función sample asigna al azar un valor a N, la cantidad de personas que puede haber en una hora determinada del día en la oficina. A veces N puede tomar el valor de uno, otras veces dos y otras veces tres .

  3. Las puertas y ventanas pueden permanecer cerradas o abiertas permanentemente o durante un periodo determinado al azar. Para simular este efecto se generó la variable cond que toma un valor al azar entre 0 y 1 que corresponden a cerrado y abierto respectivamente.

  4. El if(cond==0)… calcula la concentración con la ecuación 1 o la ecuación 2, según el valor de la variable cond.

  5. El for k … permite calcular las posibles concentraciones en los intervalos de tiempo de duración de una hora: 7am-8am, 8am-9am,…, 12pm:1pm, 1pm:2pm,2pm:3pm.

  6. Con el for i… se ejecutan 1000 cálculos de concentración bajo diferentes condiciones para cada intervalo de una hora: 7am-8am, 8am-9am,…, 12pm-1pm, 1pm-2pm,2pm-3pm.

  7. Finalmente, el programa, guarda en la matriz con los resultados de las 1000 ejecuciones para cada hora del día desde las 7 am a las 3 pm.

 

Visualización de resultados con ggplot2

Para graficar los resultados de la variable con, es necesario cargar a la sesión de trabajo de R las librerías ggplot2 y dplyr, darle formato de data frame a la variable con que es una matriz con nueve columnas y mil filas. La función stack() redefine la data frame para que tenga solo dos columnas: valores de concentración y hora del día. Esta estructura permite obtener un gráfico de cajas como el que se muestra en la figura 2.

 

library(ggplot2)
library(dplyr)
datosmod <- stack(as.data.frame(con))
levels(datosmod$ind)=c("7", "8", "9","10", "11", "12","13", "14", "15")

 

De acuerdo a la figura 2, la simulación muestra que la concentración de \(CO_2\) puede variar entre 400 y 900 ppm y que no ze alcanza la condición de estado estacionario durante el periodo que dura la simulación.

 

datosmod%>% 
  ggplot(mapping = aes(x =ind, y = values)) + 
  geom_boxplot(fill="gray",alpha=0.2)+
  geom_jitter(mapping = aes(x =ind, y =values),col="orange",alpha = .13)+
  stat_boxplot(geom = "errorbar",width = 0.2)+
  scale_y_continuous(breaks = seq(0, 1000, 200))+ #clave para modificar rotulos
  ylab("C(ppm)")+xlab("Hora")+
  labs(caption="Figura 2. Comportamiento temporal del CO2 empleando una distribución uniforme")+
  theme(plot.caption = element_text(hjust = 0.5,size=11))

 

El siguiente código es parecido al anterior solo que se supone una distribución normal para los parámetros antes mencionados, en lugar de una distribución uniforme.El resultado se puede observar en la figura 3.

 

con=matrix(0,nrow=1000, ncol=9);C=0
for (k in 1:9){
  for (i in 1:1000){
      cond=sample(c(0,1),1)#0=cerrado; 1=abierto
      N=sample(c(1,2,3),1)#no. personas dentro
      Ce=rnorm(1, mean=390, sd=5)#Ce varia entre 380 y 400 ppm
      if(k==1){Co=rnorm(1, mean=390, sd=5)}#varia entre 380 y 400 ppm
      ACH=rnorm(1, mean=0.6, sd=0.2)#ACH varia entre 0.2 y 1.0
      Q=ACH*173 #flujo en m3/s
      Resp=rnorm(1,mean=14, sd=1) #l/h
      t=seq(0, 1, length.out = 13)
      if (cond==0){
          C=Co+((N*Resp*1e6/173000)*t)}else{
             Cest=Ce+(Resp*1e06*N/(ACH*173000))
             C=Cest+((Co-Cest)*exp(-ACH*t))
      } 
      con[i,k]=median(C)
  }#fin del for i
  Co=median(con[,k])
}#fin del for k

 

datosmod <- stack(as.data.frame(con))
levels(datosmod$ind)=c("7", "8", "9","10", "11", "12","13", "14", "15")

datosmod%>% 
  ggplot(mapping = aes(x =ind, y = values)) + 
  geom_boxplot(fill="gray",alpha=0.2)+
  geom_jitter(mapping = aes(x =ind, y =values),col="brown",alpha = .13)+
  stat_boxplot(geom = "errorbar",width = 0.2)+
  scale_y_continuous(breaks = seq(0, 1000, 200))+ #clave para modifucar rotulos
  ylab("C(ppm)")+xlab("Hora")+
  labs(caption="Figura 3. Comportamiento temporal del CO2 empleando una distribución normal")+
  theme(plot.caption = element_text(hjust = 0.5,size=11))

 

Referencia:

  1. Ramos, H.S.;Magaña, V.E. & Carrera, V.J.M. (2015). Introducción a la Modelación del Agua, del Aire y del Transporte de Contaminantes en el Suelo. Universidad Juárez Autónoma de Tabasco. Villahermosa, Tabasco, México. 128 p.