Por otra parte, para determinar la absorción fotónica dentro de la placa plana, se debe resolver la ecuación de transferencia radiante (ETR) en el sistema:

En donde I representa la intensidad de radiación fotónica en el intervalo de longitudes de onda de interés, Δλ; s es la dirección y Ω el ángulo a través de los cuales viaja el haz de fotones. Por su parte, el primer sumando del lado derecho de la ecuación representa la cantidad de energía que se absorbe en el volumen de control, en el intervalo de longitudes de onda de interés, siendo к el coeficiente de absorción; el segundo sumando indica la cantidad de energía que se disipa por el fenómeno de dispersión hacia fuera (o scattering) del volumen de control, siendo а_дл el coeficiente de dispersión en el intervalo de longitudes de onda de interés; y el tercer sumando indica la cantidad de energía que gana el sistema por la dispersión hacia dentro del sistema, siendo p(Ω’→Ω) una función de fase que indica la probabilidad de incidencia de los fotones sobre el elemento de dirección ds.

La ETR se puede resolver mediante el bien conocido método de los seis flujos [24-26], el cual tiene los siguientes supuestos:

El fluido no absorbe radiación de la longitud de onda de interés.

Figura 3

Estos supuestos se han validado previamente de manera satisfactoria en el modelado global de reactores fotocatalíticos de diferentes configuraciones [7, 27] y, teniendo en cuenta el planteamiento original del SFM, son especialmente válidos para reactores de placa plana como el que se estudia en el presente trabajo [24].

Las ecuaciones del SFM son:

En donde LVRPA es la velocidad local de absorción de fotones; I₀ (en unidades de W/m²) representa la intensidad de radiación que incide sobre la película en la placa; λ_ωсоrr representa un tipo de longitud de extinción característica; x es la coordenada de dirección de los fotones entrantes; y ω_corr simboliza el albedo de dispersión corregido, el cual se define como:

Donde

Donde ω es el albedo de dispersión y se define por medio de la Ec. (22); mientras que p_f, p_b y p_s representan las probabilidades de dispersión de los fotones hacia adelante, atrás y hacia los lados en el sistema, cuyos valores, respectivamente, se tomaron como 0,11; 0,710 y 0,045, considerando una función de fase del tipo esfera grande con reflexión difusa.

En donde σ y к son los coeficientes másicos de dispersión y absorción promediados en la zona del espectro electromagnético en la que catalizador es activo (385-295 nm para el TiO₂ Aeroxide-P25), los cuales tienen los siguientes valores: σ = 1295,75 m²/kg, к = 174,75 m2/kg [27]. γ es un parámetro adimensional, definido por la Ec. (23), el cual resulta de la agrupación de términos como el albedo de dispersión y el camino óptico aparente, en la solución de la ecuación diferencial resultante del SFM:

Los otros parámetros del SFM se estiman como se muestra a continuación:

En donde C_cat es la concentración del catalizador; т y т_app representan el espesor óptico y el espesor óptico aparente; mientras que δ_SFM es el espesor de la película líquida.

La LVRPA debe aparecer explícitamente en la ecuación cinética fotocatalítica, ya que la generación de los radicales oxidantes responsables de la degradación de los contaminantes, depende de la absorción fotónica por parte del catalizador.

En el modelo se consideraron dos sistemas: la tubería y la placa plana, ambos se ilustran en las figuras 4 y 5. La tubería, de 1 pulgada nominal de diámetro (0,028 m), incluía una sección de entrada de 7 cm de largo, un codo de 90°, una sección de tubería de 14 cm de largo ubicada después del codo, y una tubería perforada de 50 cm, la cual tenía diez agujeros de 1/4" de diámetro (0,00635 m). Sólo se modeló medio sistema porque el plano xy es un plano de simetría (ver figura 4).

Figura 4

Figura 5

Cabe mencionar que las perturbaciones que se generan por la presencia de un codo y una tee (T) en la sección precedente al distribuidor perforado pueden afectar notablemente el campo de velocidades del fluido dentro del distribuidor. De esta manera, para modelar acertadamente la distribución del campo de velocidades en la entrada del mismo, fue necesario realizar un análisis en tres dimensiones en donde se consideraron estos accesorios. Además, debido a que al emplear un modelo de dos dimensiones se considera que la tubería es infinita en la dirección que va fuera del plano (fuera del papel), se tendrían problemas de continuidad en los agujeros del distribuidor, teniendo en cuenta que su área transversal es diferente a la de la sección de entrada del flujo.

Por otra parte, como el espesor de película para la operación de este tipo de reactores generalmente está entre 1 y 1,5 mm [28], se decidió hacer la simulación de la placa en dos dimensiones, considerando que, debido al flujo turbulento, las capas límites de velocidad y concentración que se encuentran de manera ortogonal a la placa tendrán un perfil más plano; mientras que a lo ancho de la placa es más difícil lograr un mezclado completo. Se tuvo en cuenta, además, que la velocidad y la fracción volumétrica del catalizador no son uniformes en los agujeros del distribuidor que alimentan a la placa, como se puede constatar con los resultados obtenidos (ver, por ejemplo, tabla 1).

Tabla 1

Orificio	Velocidad (m/s)	Fracción volumétrica x 104
1	0,455	2,46
2	0,462	2,72
3	0,476	2,82
4	0,486	2,88
5	0,488	2,90
6	0,494	2,93
7	0,502	2,96
8	0,505	2,97
9	0,511	3,01
10	0,514	3,14

Vale la pena mencionar que con la simulación 2D no es explícita la condición de no deslizamiento en la superficie de la placa ni los esfuerzos cortantes que pueda generar el viento encima de la placa. Sin embargo, como se mencionó anteriormente, el régimen de flujo turbulento hace que el perfil de velocidades en el eje y (aquel asociado al espesor de la película líquida) no tenga gradientes tan pronunciados como en un perfil laminar, lo cual está acorde con lo reportado por otros investigadores [11, 29, 30]. De esta manera, y teniendo en cuenta que existe una mayor fuente de perturbación a lo ancho de la placa debido a la presencia de agujeros en el distribuidor que alimenta el fluido a la placa, con su respectiva variación de flujo; la física a lo ancho de la placa se considera más relevante.

Ahora bien, adicional al hecho de que con la simulación 2D se ahorran muchos esfuerzos computacionales, otros investigadores que han trabajado en películas descendentes han promediado la velocidad en la dirección y (aquella ortogonal a la superficie de la placa) y analizan el comportamiento de esta velocidad promediada a lo ancho de la placa, debido a que existen mayores gradientes en tal dirección [30]. Esto ocurre a pesar de que el mecanismo de alimento del fluido a la placa es más homogéneo en comparación con el nuestro, el cual tiene una distribución a partir de tubería perforada. Por ello, en este trabajo el análisis de la velocidad es especialmente importante a lo ancho de la placa. Finalmente, en cuanto a los esfuerzos cortantes del viento, vale la pena mencionar que se asume una atmósfera estancada para la realización de las simulaciones, condición que ha sido aceptable en el modelado de reactores foto-catalíticos [11].

Por lo tanto, la geometría para la placa plana consistió en una zona rectangular de 50 cm de ancho y 1,25 m de largo. En la parte superior cuenta con diez segmentos de líneas por medio de las cuales se asumió que ingresaba el fluido. En este caso, también se modeló medio sistema porque en x = 0 se tenía una línea de simetría (ver figura 5).

Zona de entrada:

En donde u, к y s simbolizan el vector velocidad, la energía cinética turbulenta y la velocidad de disipación de la energía cinética turbulenta, como se mencionó anteriormente; U₀ es la velocidad promedio en la entrada de la tubería o en cada una de las entradas de la placa; n es el vector normal unitario; I_T y L_T representan la intensidad de la turbulencia y la longitud de escala de la turbulencia, las cuales se definen más adelante; φ_d es la fracción volumétrica de la fase dispersa y φ_d0 es la fracción volumétrica en la entrada.

Para el cumplimiento del objeto del presente trabajo, primero fue necesario estimar la distribución de la velocidad y de la concentración de la fase dispersa (el catalizador) en los orificios del distribuidor, la cual serviría como condición de frontera de entrada para la película descendente en la placa plana. Estas variables corresponden respectivamente a los valores de u y φ_d, mostrados anteriormente en el sistema de ecuaciones del modelo de mezclas (Ec. (1)-(16)).

En el reactor de placa plana de la Universidad del Valle, que se tomó como reactor real para análisis, se suele emplear un flujo de 20 L/ min para garantizar un régimen de flujo turbulento. Por lo tanto, fue necesario emplear un módulo de COMSOL apropiado para este régimen. La tubería por la cual se conduce el fluido a la placa plana es de 1 pulgada nominal (0,028 m de diámetro real); por lo tanto, el área de la sección transversal es de 6,15 x 10^-4 m² y la velocidad promedio en la entrada se fijó en 0,54 m/s. Se consideró que el perfil de la velocidad en la entrada del sistema estaba completamente desarrollado debido a que la sección de la entrada se ubicó aproximadamente a 50 cm del último accesorio de la tubería que genera perturbaciones de flujo importantes (codo). Cabe recordar que para tuberías circulares por las que circula un flujo turbulento, la longitud de entrada es aproximadamente diez veces el diámetro de la misma [31]; es decir, en nuestro caso, es de aproximadamente 28 cm. Por lo tanto, se puede esperar que en la zona de entrada el perfil de velocidad esté completamente desarrollado.

Las propiedades del fluido se tomaron como las del agua a la temperatura ambiente, debido a que la suspensión estaba diluida. De esta manera, la densidad se fijó en 1000 kg/m³ y la viscosidad dinámica en 1 x 10^-3 kg/(m x s). Para tuberías, el número de Reynolds (Re_T) se calcula como sigue:

En donde V_prom es la velocidad promedio; p es la densidad; ц es la viscosidad y D es el diámetro interno de la tubería. Con la Ec. (40) se obtuvo que el número de Re para la tubería es de 15.157, el cual está dentro del rango de flujo turbulento.

Por otra parte, las propiedades de la fase dispersa (TiO₂ P-25) necesarias para ejecutar la simulación son: densidad: 3800 kg/m³ y diámetro de partícula: 1,1 x 10^-6 m a pH neutro [32]. Cabe anotar que este diámetro corresponde a las partículas aglomeradas. Se empleó este valor para considerar los efectos de atracción por fuerzas electrostáticas y de Van Der Waals que ocurren entre las partículas de TiO2 [33]. Se consideró una carga inicial de TiO2 de 1 g/L, aplicados a 20 L de la suspensión, teniendo en cuenta las condiciones de carga del catalizador en la sección previa a la tubería, de manera que la fracción volumétrica de la fase dispersa se tomó como 2,6316 x 10^-4.

Para lograr una convergencia más rápida, primero se resolvió el problema de la tubería perforada, considerando el flujo como homogéneo y laminar. Posteriormente, los resultados obtenidos con estas condiciones fueron usadas como valores de arranque (o iniciales) para nuestro caso de interés (flujo turbulento bifásico).

Por otra parte, el valor de la escala de longitud de la turbulencia se tomó como 1,96 x 10^-3, teniendo en cuenta que equivale al 7% del diámetro en una tubería circular [34, 35]. Este parámetro es una cantidad física que describe el tamaño de los remolinos que se generan en el flujo. De igual manera, la intensidad de la turbulencia, definida como la relación entre el promedio de las componentes fluctuantes de la velocidad y la velocidad media, se fijó en 5%, como se recomienda para flujos turbulentos [34]. Tanto la escala de la longitud de la turbulencia y la intensidad de la turbulencia son útiles para estimar las propiedades de la turbulencia en las fronteras de entrada, tales como la energía cinética turbulenta к y su disipación ε.

Además, se hicieron varias simulaciones para determinar la independencia de los resultados respecto al tipo de malla que se empleó. Para ello, se analizaron tanto la velocidad de la mezcla como la fracción volumétrica del catalizador en la salida de cada uno de los orificios del distribuidor. Se seleccionaron estas dos variables para realizar este análisis de sensibilidad, porque el posterior estudio de la placa dependía directamente del valor que alcanzaban en los orificios. Vale la pena mencionar que el enmallado se hace más riguroso en zonas en donde la perturbación del fluido era importante, tales como los orificios en el distribuidor y la sección de ingreso del fluido en la placa.

En la figura 6 se muestra el detalle de la malla empleada para modelar el distribuidor perforado. El análisis de sensibilidad mostró que con una malla de 331.613 elementos la solución era un 2,33% diferente de la obtenida con una de 190.567 elementos de la misma clase. De esta manera, se tomaron los resultados de la malla de 331.613 elementos (conformada por 238.571 elementos tetraédricos, 89.545 en forma de prisma y 3497 piramidales), para establecer las condiciones de entrada de la placa plana.

Figura 6

De esta manera, los perfiles de velocidad y distribución del catalizador en la tubería que se obtuvieron con el módulo de mezcla de COMSOL fueron los descritos en las figuras 7 y 8, respectivamente.

Figura 7

Figura 8

El promedio de la velocidad y la fracción volumétrica del catalizador que se alcanzó en la salida de cada agujero del distribuidor se muestra en la siguiente tabla, en la cual el orificio uno es el que está en el borde derecho de la figura 7 (a la entrada del sistema), y el orificio diez el que está en el borde izquierdo (extremo final):

Como se nota en la figura 7, existe un gradiente pronunciado de la velocidad en la entrada del distribuidor. Esto se debe a la presencia del accesorio de bifurcación ubicado justo antes del mismo. Vale la pena anotar que con la geometría en 3D se obtuvieron valores razonables de la velocidad de flujo a la salida de los orificios de los distribuidores (ver tabla 1), de manera que se satisfizo la ecuación de continuidad. Esto se debe al principio de conservación de la materia, dado que el flujo volumétrico total saliente de los orificios (y el flujo másico, tomando la densidad constante) corresponde con razonable aproximación con el flujo de entrada al distribuidor. En la aproximación 2D que se realizó (resultados no presentados aquí), los valores de la magnitud de la velocidad eran de varios órdenes de magnitud menor. Sin embargo, como se nota en las figuras 7 y 8, la orientación de los accesorios y los orificios de la tubería no permiten realizar una aproximación completamente acertada a un modelo 2D.

Además, de acuerdo con la tabla 1 y la figura 7, el porcentaje de mala distribución definido como la variación porcentual del flujo entre el primer y último orificio del distribuidor es de aproximadamente 13%. Para el diseño de este tipo de accesorios usualmente se emplea un valor nominal de 5% para este parámetro, con el fin de que la distribución sea lo suficientemente uniforme en el sitio de interés [36]. De esta manera, se evidencia que debería mejorarse el diseño del mismo para aumentar la homogeneidad en el reactor de placa plana.

Después de haber resuelto el campo de velocidades y concentración del catalizador en la tubería, se tomó el promedio de la velocidad de la mezcla bifásica y la fracción volumétrica del catalizador, en cada uno de los orificios del distribuidor, como condiciones de frontera de entrada para el estudio en la placa plana.

Para películas descendentes que fluyen a través de una placa, el número de Reynolds (Rep) se calcula como sigue [36, 37]:

En donde Q es el flujo volumétrico y a es el ancho de la placa (aproximadamente 1 m). Por su parte, ρ y μ simbolizan respectivamente la densidad y la viscosidad, como se mencionó anteriormente. Aplicando la Ec. (41) se obtuvo que el número de Reynolds para el caudal de trabajo es aproximadamente 1370. En la literatura se muestra que con este valor el flujo se encuentra en régimen turbulento [37]. Por lo tanto, para el estudio de la placa también se empleó el módulo de mezcla para flujo turbulento de COMSOL.

En este caso, para las entradas del sistema, la intensidad de la turbulencia tuvo un valor de 5% y la escala de la longitud de la turbulencia se fijó en 0,09, tal como se recomienda para este tipo de configuraciones [34].

De igual forma, se realizó un examen de sensibilidad respecto al tipo de malla que se empleó para resolver el sistema. Se analizó el comportamiento del promedio de la concentración del catalizador, debido a que en el caso 2D esta era la variable de mayor interés. El promedio de la fracción volumétrica se realizó sobre 11 líneas de corte que se definieron en el modelo, las cuales atravesaban la totalidad del ancho de la placa y estaban ubicadas a ciertas alturas de la misma, como se indica en la figura 9.

Figura 9

Los resultados obtenidos con los tres tipos de mallas son similares en cuanto a la forma de la curva; sin embargo, los valores de la fracción volumétrica que arrojan las mallas de 1596 y 1074 elementos triangulares son muy similares tanto en forma como en magnitud. De hecho, la diferencia promedio entre los mismos es de 1,67%. De esta manera, se empleó la malla de 1596 elementos triangulares para realizar los cálculos subsiguientes. El detalle de la placa con el enmallado se muestra en la figura 10. Por otra parte, en la línea de corte ubicada en y = 1,25 se nota que la fracción volumétrica es mucho más baja. Esto se debe a que en esta zona se encuentran varias secciones que no son cubiertas por el catalizador, lo cual produce un descenso en la concentración promedio del mismo.

Figura 10

Los perfiles de velocidad y concentración del catalizador en toda la placa se muestran a continuación, en los cuales se ha proyectado la solución hacia el extremo izquierdo de la placa, teniendo en cuenta que en la línea x = 0 estaba la condición de frontera de simetría.

Tanto en la figura 11 como en la figura 12se nota una zona en la que la magnitud de las variables es baja. Esta zona aparece porque justo en el plano central del distribuidor no hay agujeros; estos se encuentran a partir de, aproximadamente, 2,4 cm después de la bifurcación en la entrada del distribuidor.

Figura 11

Figura 12

En la figura 12 se evidencian varias zonas en las que la concentración del catalizador es diferente. En la figura 13 se muestran las siete divisiones que a juicio de los autores representan las zonas más importantes en donde la concentración del catalizador tiene un valor similar.

Figura 13

Con la opción de promedio sobre superficies de COMSOL, se calcularon los promedios de las fracciones volumétricas de las zonas que se muestran en la figura 13. Para ello, fue necesario definir un nuevo componente en el modelo de COMSOL y usar la opción "extrusión general", con el que se pudieron realizar los cálculos sobre la base de los resultados obtenidos en la figura 12. De esta manera, las fracciones volumétricas correspondientes a las zonas demarcadas de la figura 13son:

Tabla 2

Zona	Ancho (m)	Largo (m)	Fracción volumétrica de catalizador x 104
1	0,322	1,18	4,8843
2	0,101	1,18	4,612
3	0,077	0,678	4,5317
4	0,077	0,322	4,4287
5	0,077	0,25	3,963
6	0,322	0,07	4,624
7	0,101	0,07	4,3886

Para cada una de esas zonas se calculó el promedio de la absorción de fotones que se tendría para ciertas condiciones de intensidad de radiación.

Después de haber determinado los perfiles de concentración dentro de la placa plana, se procedió a realizar la estimación de la velocidad local de absorción de fotones (conocida como LVRPA por sus siglas en inglés) en cada una de las zonas delimitadas en la figura 13. Para ello, la ecuación de transferencia radiante se resolvió mediante el SFM en una rutina de simulación con Visual Basic.

Queda para una posterior investigación el acople simultáneo en COMSOL de la determinación de la distribución del catalizador en el reactor (resultados que se obtuvieron en el presente trabajo) con la determinación de la absorción y dispersión de fotones (scattering) dentro del sistema. Por lo pronto, para cumplir el objetivo del presente estudio fue satisfactorio emplear el SFM, con el cual se pueden obtener resultados comparables a los que arrojan otros modelos más sofisticados [24, 38]. En este caso, por efectos de estudio, se simuló solamente la componente de la radiación directa, con una intensidad de radiación de 30 W/m², la cual aparece como referencia en la mayoría de trabajos sobre fotocatálisis solar heterogénea [17, 27]. Por su parte, el valor de δ_SFM es aproximadamente un milímetro, conforme a lo reportado por Colina-Marquez et al. [28], dado que el flujo volumétrico considerado para el reactor de placa plana analizado fue de 20 L/min.

El catalizador seleccionado fue TiO₂ Degussa P-25 (80% anatasa, 20% rutilo), para el cual se considera una función de fase del tipo esfera grande con reflexión difusa, en la que las probabilidades de dispersión son: p_f = 0,11, pb = 0,710 y ps = 0,045 [24, 27]. Cabe anotar que en COMSOL se deben emplear polinomios de Legendre para definir la probabilidad de dispersión anisotrópica de las partículas del catalizador. No aparece la opción explícita de esfera grande con reflexión difusa.

De esta manera, se calculó tanto el promedio de la LVRPA como la absorción global de fotones (OVRPA) en cada zona demarcada de la figura 13. Para ello, se hizo uso del SFM y de las siguientes definiciones, en donde z es la longitud perpendicular a la placa plana del reactor; mientras que A_zona y L_zona, representan el ancho y el largo de las zonas delimitadas de la figura 13:

Los resultados se muestran en la tabla 3.

Tabla 3

Zona	LVRPAprom (Wm-3)	OVRPA(W)
1	3944,66	1,50
2	3742,97	0,45
3	3683,07	0,19
4	3605,95	0,09
5	3253,29	0,06
6	3751,91	0,08
7	3575,83	0,03

Como se puede observar en la tabla 3, existe una diferencia porcentual de más del 20% entre el mayor y el menor valor del promedio de la absorción fotónica en la placa. Debido a que el orden de la velocidad de reacción fotocatalítica respecto a la absorción fotónica suele tomar valores entre 0,5 y 1 [39], se puede esperar que la variación en la degradación de los contaminantes en estas zonas del reactor oscile entre el 10 y el 20%, lo cual implica que esta distribución afecta la eficiencia del reactor fotocatalítico.