Objeto de estudio

Se usó un tubo de PVC de 3.0 m de longitud y 15.5 cm de diámetro (D _L ). En la entrada se conectó un sistema de alimentación de aire consistente de un ventilador centrífugo y una válvula de mariposa, donde se marcaron cinco posiciones para el control del flujo. En la salida se colocó una boquilla reductora de diámetro de 7 cm (D _e ), donde se evaluó la velocidad lineal del aire. A 2.0 m de distancia del punto de entrada se colocaron dos rejillas de alambre separadas entre sí por 30.0 cm. El espacio entre las mallas se llenó con esferas de unicel y representó un medio macroporoso semejante a una unidad de empaque de frutos colocada en un proceso de enfriamiento, con ingreso lateral de aire. Se propuso que la caída de presión ((P, Pa) puede representarse con el modelo de Ergun (Ecuación 1):

donde ( es densidad (kg·m^-3) del aire a la presión promedio entre la entrada y la salida; G ₀ es flux másico (kg·s^-1·m^-2) del fluido, expresado como el producto [(·v ₀ ]; v ₀ es la velocidad superficial, que representa la velocidad lineal promedio (m·s^-1) del fluido que habría en la caja si no hubiera producto; D _p es el diámetro (m) de la unidad de relleno (una esfera en este caso), expresado por el cociente (6·a _v ^-1 ); a _v es la superficie específica del relleno, dado como el cociente de la superficie total del relleno (sumando todas las unidades) entre el volumen del relleno; L es la longitud (m) del lecho poroso en la dirección del flujo; Ɛ es la porosidad o fracción hueca (adimensional) y Re es el número de Reynolds (adimensional). Adicionalmente, k ₁ y k ₂ son constantes características de la ecuación de Ergun, con valores de 150 y 1.75, respectivamente.

Estructura de investigación

Se desarrollaron dos etapas de investigación. En la primera, se evaluó la eficacia de la ecuación de Ergun para estimar la caída de presión ((P _est ) a través de la unidad de empaque. En la segunda, se calibraron las constantes k ₁ y k ₂ de la Ecuación (1) para optimizar la predicción (P _est .

Evaluación de la ecuación de Ergun

Se realizó un experimento a dos condiciones térmicas, una en el rango de 10 a 12 °C y la otra de 24 a 26 °C. Se probaron cinco tamaños de objetos empacados consistentes de esferas de unicel, las cuales se denominaron genéricamente como partículas u objetos de relleno. En cada caso se colocó en el tubo un número suficiente de partículas para conseguir el llenado equivalente al que se tendría normalmente con material biológico. Así, los diámetros (D _p ) usados y el número de objetos alojados (N) fueron: 2, 3, 4, 5 y 7 cm, con 772, 215, 86, 42 y 14 objetos, respectivamente. Para cada tamaño de partícula se evaluó la Ɛ, mediante la Ecuación 2, donde V _L es volumen del lecho poroso (0.00566 m³) y V _p el volumen de un objeto de relleno (Ecuación 3).

Mediante la Ecuación 4 y con base en la Figura 1, se determinó el diámetro hidráulico (D _h ) del área de flujo de aire en el lecho poroso para cada tamaño de partícula.

Figura 1 Identificación del área de flujo en el lecho poroso de partículas. 

Adicionalmente, para cada D _p se probaron cinco condiciones de flujo de aire en el tubo (F ₁ , F ₂ , F ₃ , F ₄ y F ₅ ), controladas mediante las posiciones de válvula de mariposa. Esta estructura dio como resultado 25 condiciones de experimentación, que representaron tratamientos (D _p -F _j ). Cada tratamiento se sometió, por triplicado, a medición experimental de caída de presión (∆P _exp ) a través de un sistema de adquisición de datos (LabQuest, Vernier Software & Technology, LLC, USA) equipado con un sensor de medición de presión manométrica. Adicionalmente, se midió la velocidad del aire (v _e , m·s^-1) en la salida de la unidad de empaque con un anemómetro (Kestrel 4000 Pocket Weather Meter, USA). Se calculó el caudal (V˙, m³·s^-1) con la Ecuación 5, la velocidad superficial (v ₀ , m·s^-1) con la Ecuación 6, el flux másico (G ₀ , kg·s^-1·m^-2) con la Ecuación 7 y el número de Reynolds (Re) con la Ecuación 8; donde D _L es el diámetro del lecho (0.155 m), ( densidad del aire (kg·m^-3) y μ viscosidad (Pa·s) del mismo fluido, estas últimas evaluadas en función de la temperatura (^{Kundu, Cohen, & Dowling, 2012}).

Finalmente, se usó el modelo de Ergun (Ecuación 1) para obtener un valor estimado de la caída de presión ((P _est ), considerando los valores originales de k ₁ = 150 y k ₂ = 1.75.

Calibración de la ecuación de Ergun

Con base en la información experimental, se realizó la calibración de los parámetros k ₁ y k ₂ de la Ecuación 1, en forma numérica, empleando un enfoque de mínimos cuadrados no lineales ordinarios; para ello se usó la subrutina lsqnonlin.m, disponible en el paquete de herramientas de optimización de MatLab^® (^{The Mathworks Inc., 1990-2006}), con el algoritmo de optimización de Levenberg-Marquardt. En este proceso, se tomaron como referencia los valores nominales k ₁ = 150 y k ₂ = 1.75. Para ambos parámetros se definieron las restricciones de dominio dadas por las Ecuaciones 9 y 10, respectivamente. Asimismo, la función f(x) a minimizar fue la mostrada por la Ecuación 11, donde (P _exp corresponde a los datos experimentales y (P _est a los datos calculados con la Ecuación 12:

Análisis de datos

La evaluación experimental de la caída de presión se hizo en forma congruente con un diseño factorial 2 × 5 × 5, donde los factores de variación fueron la temperatura (rangos de 10 a 12 y 24 a 26 °C), el diámetro de elementos de empaque (2, 3, 4, 5 y 7 cm) y el control del flujo, dado por cinco posiciones de apertura de válvula en la alimentación. Los datos se sometieron a análisis de varianza y pruebas de comparación de medias mediante el estadístico de Tukey (α = 0.05). Los datos experimentales de la caída de presión se compararon con las estimaciones del modelo de Ergun a través de pruebas de comparación de medias hechas con el estadístico de Tukey (α = 0.05). Todo el análisis se realizó con el paquete estadístico Statistical Analysis System, ver. 8.0 (^{SAS, 1999}).

Fracción hueca (Ɛ) y diámetro hidráulico (D _h )

La Ɛ varió entre 0.43 y 0.56, y se incrementó de forma lineal conforme aumentó el diámetro de partícula (R² = 0.9873, Figura 2). En el caso del D _h , la variación se ubicó en el rango de 0.55 a 1.91 cm y también se tuvo un aumento lineal con D _p , aunque en este caso los valores se basan en el arreglo teórico de las partículas (Figura 1, Ecuación 4); por lo cual el coeficiente de determinación (R²) fue igual a la unidad.

Figura 2 Relación de la fracción hueca (Ɛ) y el diámetro hidráulico (D _h ) con el diámetro de partícula de relleno (D _p ) en el lecho empacado. 

Flujo volumétrico

El flujo volumétrico varió entre 8.8 y 35.0 L·s^-1 y las posiciones de la válvula de control permitieron flujos que fueron significativamente distintos entre ellos (P ≤ 0.05, Cuadro 1). A temperatura mayor se tuvo 8.5 % más flujo. En el aire, a medida que aumenta la temperatura la densidad se reduce, pero la viscosidad aumenta (^{Kundu et al., 2012}). El cociente de la viscosidad entre la densidad representa una difusividad de momento (δ, m²·s^-1) y, con base en la Ley de Newton de la Viscosidad (^{Bird, Stewart, & Lightfoot, 2006}), el incremento en la temperatura debió corresponder a mayor transporte de momento entre el fluido y las partículas empacadas, dificultando el flujo. Sin embargo, la evidencia empírica mostró que en realidad hubo mayor velocidad, lo que sugirió el reacomodo de las partículas al pasar el sistema de experimentación de la condición de baja a la de alta temperatura.

En forma similar, el diámetro de las partículas de relleno afectó el flujo, y a medida que el tamaño de las mismas aumentó el caudal también fue mayor (P ≤ 0.05). Esta situación era esperada, ya que conforme el diámetro de partícula fue mayor, incrementaron la porosidad y el diámetro hidráulico (Figura 2). La tortuosidad es una medida de la falta de trayectorias lineales y falta de paralelismo de corrientes en el flujo de un fluido en un lecho empacado (^{Matyka, Khalili, & Koza, 2008}). Así, a medida que se incrementaron Ɛ y D _h , el fluido en movimiento encontró áreas mayores de flujo, lo que significó menor tortuosidad, favoreciendo con ello el incremento del caudal. Por otro lado, hubo interacción significativa (P ≤ 0.05) entre los niveles de apertura de válvula y el diámetro de partícula (Cuadro 1); pues con la apertura mayor, el diámetro más grande de las esferas de relleno permitió mayor flujo, pero a menor apertura, el diámetro de los objetos de relleno no causó variación del caudal (Figura 3).

Figura 3 Efecto de interacción entre la apertura de válvula y el diámetro de partícula en la afectación del flujo volumétrico. Los símbolos *, ** y ns indican modificación significativa (P ≤ 0.05), altamente significativa (P ≤ 0.01) y no significativa (P > 0.05), respectivamente. 

Número de Reynolds

El número de Reynolds varió entre 2,054.5 y 8,615.7; sin embargo, en columnas de relleno, la transición de flujo laminar a flujo turbulento ocurrió con valores de entre 10 y 1,000 (^{Bird et al., 2006}), lo que indica que hubo flujo turbulento en todas las condiciones probadas. Al igual a como ocurrió con el flujo volumétrico, el número de Reynolds se afectó por la temperatura y fue mayor en la condición de 24 a 26 °C, aunque el contraste fue de únicamente 7.5 % en relación con el valor observado en la condición de 10 a 12 °C, y puede no tener importancia práctica. Asimismo, el número de Reynolds se afectó por el tamaño de las partículas empacadas, y se tuvieron valores significativamente mayores con las partículas grandes que con las más pequeñas (P ≤ 0.05, Cuadro 1); lo que se explica por los diámetros hidráulicos más grandes y la mayor porosidad. De igual manera, en forma congruente con el mayor flujo obtenido mediante una mayor apertura de válvula, el número de Reynolds incrementó a medida que el caudal fue más grande.

Caída de presión experimental

La caída de presión varió entre 50 y 555 Pa. Esta pérdida fue mayor de manera proporcional con el caudal, condición en que hubo mayor turbulencia; la cual se estableció por números de Reynolds más grandes. Sin embargo, la caída de presión se redujo a medida que se incrementó el tamaño de las partículas del lecho empacado, y esto confirmó que hubo menor tortuosidad en la medida en que aumentó el diámetro hidráulico. La caída de presión también se afectó por la temperatura. Los mayores valores ocurrieron con la menor condición térmica, y se cree que esto ocurrió por un reacomodo de las partículas que permitió mayor flujo y menor caída de presión en la temperatura más alta.

Predicción de la caída de presión

La caída de presión es uno de los elementos más importantes en el diseño de un proceso de enfriamiento, pues el dimensionamiento de ventiladores en sistemas con aire forzado se basa en dicho parámetro (^{Amanlou & Zomorodian, 2011}). Aunque eventualmente esta variable puede determinarse experimentalmente, resulta de mayor utilidad disponer de relaciones matemáticas que permitan una predicción adecuada. En el presente trabajo se evaluó el desempeño del modelo de Ergun (Ecuación 1) con las constantes originales k ₁ = 150 y k ₂ = 1.75, para realizar la predicción de la caída de presión en el sistema empacado utilizado en la fase experimental. Los datos de predicción se compararon con la información experimental, y se encontró que en las dos condiciones térmicas hubo discrepancia significativa (P ≤ 0.05) entre ambos conjuntos de datos y, de acuerdo con la Figura 4, el modelo de Ergun en su formato original claramente subestimó las mediciones, con errores promedio de predicción del orden de 72.2 %.

Figura 4 Representación 1:1 entre los datos de caída de presión calculados con la ecuación original de Ergun y los datos experimentales. Los recuadros A y B corresponden a las condiciones térmicas de 10 a 12 y 24 a 26 °C, respectivamente. Los símbolos ●, ○ y □ indican distintas repeticiones. La diagonal de referencia muestra la condición de similitud entre datos experimentales y calculados. 

En adición, la aplicación de rutinas de comparación de medias mediante el estadístico de Tukey (α = 0.05) mostró que, en tanto que a nivel experimental se registró una caída de presión promedio de 266.0 Pa, el modelo original de Ergun predijo un valor significativamente distinto de 72.0 Pa. Al someter los datos a rutinas de comparación de medias se obtuvo una diferencia mínima significativa honesta (DMSH) igual a 49.5 Pa (Cuadro 2), que resultó menor al error de predicción; lo que confirmó un desempeño pobre del modelo original para reproducir la información experimental. Esta situación fue congruente con el hecho de que la ecuación de Ergun, en su formato original, describe la región de flujo laminar y transicional en lechos de partículas empacadas (^{Bird et al., 2006}). No obstante, en el presente trabajo se usaron condiciones de flujo en régimen turbulento (ver subsección de número de Reynolds), lo que explica la discrepancia entre los valores estimados y calculados.

En tal sentido, resultó necesario desarrollar rutinas de calibración de las constantes k ₁ y k ₂ de la Ecuación 1. La aplicación del enfoque de mínimos cuadrados no lineales ordinarios a través de la subrutina lsqnonlin.m del paquete de herramientas de optimización de MatLab^® (^{The Mathworks Inc., 1990-2006}) produjo los valores de k ₁ = 10,628.67 (± 4,636.97) y k ₂ = 3.44 (± 0.96). Cuando estos parámetros calibrados se usaron en la Ecuación 1, se obtuvo una estimación mejor de la caída de presión en el lecho empacado (Figura 5), con valores promedio de 250.9 Pa; los cuales, aunque difirieron aún en un 29.8 % respecto de los datos experimentales, en la aplicación de rutinas de comparación de medias se observó que con la ecuación de Ergun calibrada se logró la predicción de la caída de presión que fue estadísticamente igual al valor promedio de 266.0 Pa obtenido en forma experimental (Cuadro 2). Lo que se confirmó con un valor de DMSH (49.5 Pa) mayor al error de predicción. Sin embargo, la discrepancia que aún se observa sugiere que puede haber elementos del sistema estudiado que no se han tomado en cuenta, lo cual es una situación común; ya que por lo general cualquier modelo matemático puede representar algunas características del fenómeno pero no todas (^{Vargas-Castañeda, Ibáñez-Castillo, & Arteaga-Ramírez, 2015}).

Figura 5 Representación 1:1 entre los datos de caída de presión calculados con la ecuación de Ergun calibrada y los datos experimentales. Los recuadros A y B corresponden a las condiciones térmicas de 10 a 12 y 24 a 26 °C, respectivamente. Los símbolos ●, ○ y □ indican distintas repeticiones. La diagonal de referencia muestra la condición de similitud entre datos experimentales y calculados. 

Adicionalmente, para verificar los resultados anteriores se usó el algoritmo de evolución diferencial (^{Price, Storn, & Lampinen, 2005}), el cual es considerado como un método de optimización global conocido por su eficacia y eficiencia para resolver problemas multimodales. Los valores obtenidos fueron nuevamente k ₁ = 10,628.67 y k ₂ = 3.44, que son los mismos que los obtenidos mediante mínimos cuadrados no lineales. Esto permite deducir que la calidad de la estimación fue correcta. Por tanto, la dinámica de flujo en el sistema tubular usado como dispositivo de enfriamiento se comporta de acuerdo con la Ecuación 13, que constituye una ecuación de Ergun calibrada:

No obstante, con base en la Figura 5, a medida que la caída de presión fue mayor se observó un incremento en la discrepancia entre datos experimentales y calculados; lo cual fue más evidente a partir de 400 Pa, por lo que este valor se consideró el límite de utilidad del modelo propuesto.