Ernesto Javier Muñoz Suárez

Escuela de Procesos y Energía, Facultad de Minas, Universidad Nacional de Colombia.
ejmunoz0@unal.edu.co

Fecha Recepción: 19 de junio de 2012
Fecha Aceptación: 6 de diciembre de 2012

En este trabajo se plantea un modelo Euleriano - Lagrangiano mejorado con el método de los elementos discretos y se simula un campo de flujo formado por una fase líquida y una sólida dispersa cuyas densidades son ligeramente diferentes, el objeto de la simulación es evaluar el comportamiento de la fase dispersa formada por 7050 partículas distribuidas según la distribución normal en 10 tamaños diferentes. De la simulación fue posible encontrar la expansión del lecho, como se distribuye la fase dispersa y una evaluación de los fenómenos que producen la tendencia a la segregación y/o al entremezclado producto de la respuesta diferente que hacen las partículas de distinto tamaño a la fuerza de arrastre consecuencia del flujo ascendente.

Palabras clave: Euleriano-Lagrangiano, DEM, flujo multifase, distribución de partículas, segregación, entremezclado.

This paper presents a simulation of flow field formed by a liquid phase and a solid dispersed which densities are slightly different, the simulation was done using a Eulerian - Lagrangian model improved with the method of discrete elements, the simulation was performed to evaluate the dispersed phase particles behavior formed by 7050 particles distributed according to the normal distribution in 10 different sizes. From the simulation was possible to find the bed expansion, the distribution of particles from dispersed phase into control volume, and an assessment of the phenomena that produces trend toward segregation and/or intermixing of particles with different size like consequence drag force differences by result of upward flow.

Keywords: Eulerian-Lagrangian, DEM, multiphase flow, particles distribution, segregation, intermixing.

Los flujos multifásicos están presentes en muchas de las aplicaciones de ingeniería, lo cual, constituye un campo de estudio bastante amplio para los investigadores. Sin embargo, el análisis teórico de estos flujos suponen grandes desafíos, y en respuesta a ellos se han desarrollado varias técnicas que vienen desde la Dinámica de Fluidos Computacionales (CFD) tales como, el modelo Euleriano-Euleriano, el modelo Euleriano- Lagrangiano y el modelo VOF (volume of fluid) [1-4]. Con el eventual crecimiento en los avances de las computadoras que permiten mayor velocidad de procesamiento, ha cobrado importancia desarrollar técnicas que permitan el estudio en el nivel discreto y de este modo modelar los fenómenos que ocurren a dicha escala, de allí que el Método de Elementos Discretos (DEM) [5] cobre importancia. El DEM es un esquema de modelado computacional que tiene en cuenta los desplazamientos y rotaciones de los elementos separadamente, reconociendo nuevos contactos y calculando su progreso en sistemas de muchos cuerpos (manybodies) [6]. El siguiente paso lógico en el camino de modelar los flujos multifásicos nace al acoplar las técnicas de CFD con la técnica del DEM desarrollando un esquema Euleriano-Lagrangiano que permite modelar la interacción fluido-partículas y partículapartícula o interacción en las cuatro vías [7]. El acople entre estas dos técnicas han creado una poderosa herramienta que ha permitido abordar muchos de los problemas que involucran campos de flujos con fases dispersa granulares o de burbujas [8-11].
En la operación de los reactores donde coexisten una fase continua y otra dispersa, es decir un campo de flujo multifásico, es importante tener la habilidad de predecir la expansión del lecho y su tendencia a la segregación o al entremezclado, dada su importancia en el desempeño del intercambio de masa entre las fases y la reacción que suceda entre las mismas [12], además otro comportamiento que se caracteriza en estos flujos son los aglomerados y clúster de partículas que son formados por las fuerzas de cohesión y por los efectos del campo de flujo, respectivamente [7]. La hidrodinámica del flujo multifásico ha sido ampliamente abordada, sin embargo, bajo las condiciones de baja velocidad donde las partículas tienen una densidad ligeramente mayor que la fase continúa ha sido estudiada poco, puesto que este es de especial interés en el estudio de reactores anaeróbicos de flujo ascendente de biopelícula [13,14]. En los reactores anaeróbicos de biopelículas los microorganismos se aglomeran en colonias que tienen una frontera definida [15], lo que configura una fase dispersa dentro del campo de flujo con dicha condición. Basado en lo anterior, en este trabajo se busca modelar un campo de flujo como se presenta en los reactores anaeróbicos usando el esquema Euleriano-Lagrangiano potenciado con el DEM que permita modelar la interacción en las cuatro vías (interacción fluido-partículas e interacción partícula-partícula), en un campo de flujo ascendente en régimen laminar con una fase dispersa que tiene una distribución normal en cuanto al tamaño de partícula. Con el objeto de predecir la expansión del lecho, la distribución de las partículas al interior del reactor, las tendencias a la segregación o al entremezclado por causa de las fuerzas de arrastre.

Se construirá un modelo usando la aproximación Euleriana-Lagrangiana en una geometría de dos dimensiones. Para la solución se optará por el método CFD para el caso de la fase Euleriana y con el DEM para la fase Lagrangiana, donde el agua se constituirá en la fase Euleriana comportándose como un fluido newtoniano en un régimen de flujo laminar a temperatura constante. Se construirá una fase Lagrangiana con 7050 partículas de un material de densidad ligeramente mayor a la del agua (1070kg/m³), las partículas que formarán dicha fase tendrán una distribución de tamaño normal discreta formada por 10 clases de partículas (Tabla 1), que estarán presentes dentro del reactor y no habrá intercambio con los alrededores, además no se permitirá el intercambio de masa entre partículas y el agua.

El acoplamiento entre la fase continua y la fase dispersa será en cuatro vías para el caso del intercambio de momentum, es decir, se consideran las fuerzas de arrastre entre las dos fases, las fuerzas de elevación que actúan sobre las partículas causadas por los gradientes de presión en el agua y las interacciones que ocurren entre las partículas con su pares.

Modelo fase Euleriana
La fase continua Euleriana estará descrita por el balance de masa que puede ser expresado según la Ecuación 1 sí se consideran las restricciones expuestas anteriormente.

Donde α_L es la fracción de agua residual y ρ_L la densidad de agua residual, es el vector velocidad. El balance de momentum sobre el dominio de la fase continúa se expresa en la Ecuación 2 (ecuación de Navier-Stokes para flujos newtonianos incomprensibles).

Donde P es la presión y g es la gravedad respectivamente, μ es la viscosidad, V_Co es el volumen de control infinitesimal de la celda k que contiene las partículas, es el término de acoplamiento que describe la interacción debida a las fuerzas de arrastre y de elevación que ocasionan el intercambio de momentum entre las dos fases.

Modelo fase Lagrangiana
La fase dispersa Lagrangiana será descrita por la ecuación de movimiento de Newton, la ecuación que tiene en cuenta la rotación de las partículas y la ecuación de posición para cada una de las partículas individualmente, siendo éstas las Ecuaciones 3, 4 y 5 que forman el método de elementos discretos.

Donde m_n y son la masa y la velocidad translacional de la partícula, es el término de acoplamiento que describe la interacción debida a las fuerzas de arrastre y de elevación que ocasionan el intercambio de momentum entre las partículas y el agua, es la aceleración de la gravedad la cual configura la fuerza de cuerpo que actúa sobre las partículas. Finalmente es el término asociado a las interacciones que se dan entre dos o más partículas, este último término es lo que le da importancia y la robustez al DEM, puesto que permite evaluar la interacción individual de cada partícula con el resto existente en el dominio, y de este modo conformar el comportamiento de conjunto que se observa macroscópicamente [6].

La posición en el espacio de las partículas es dada por la Ecuación 4 y la Ecuación 5 describe su rotación, esto se debe a que las partículas tienen extensión espacial, y por ende dicha rotación influye en las trayectorias de las partículas [16].

Donde I_n y son el momento de inercia y la velocidad angular de las partículas, y es el torque generado por fuerzas externas no centradas que someten a las partículas y están asociados a las fuerzas tangenciales producto de las colisiones. El término , el cual considera la interacción fluido-partícula se define en la Ecuación 6.

Donde es el gradiente de presión de la partícula n-esima ubicado en la posición y da cuenta de la fuerza de elevación, v_pn es el volumen de las partículas, y es el coeficiente local de transferencia de momentum entre la partícula n-esima y el agua que tiene en cuenta el aporte de la fuerza de arrastre, y se modela usando la correlación experimental propuesta por Syamlal y O'brien 1987, expresada por la Ecuación 7.

Donde V_ter y d_n son la velocidad terminal y el diámetro de las partículas, C_Ds es el coeficiente de arrastre asociado al tipo de geometría de las partículas y R_es es el número de Reynolds de partícula. La velocidad terminal es calculada empleando la ecuación de Garside y Al Dibouni 1977 (Ecuación 8). Adicionalmente, el coeficiente de arrastre (Para un numero de Reynolds característico del campo de flujo) y el número de Reynolds de partícula se evalúan utilizando las Ecuaciones 9 y 10, respectivamente.

Con

Modelo interacción partícula - partícula
La interacción partícula-partícula presente en el campo de flujo se modelará usando el modelo de esfera blanda [5], cuyo objetivo es determinar la fuerza neta que actúa sobre la partícula como consecuencia de contacto con otra partícula aproximando dicho contacto a un sistema de resorte y amortiguador, esta aproximación permite que sea posible evaluar el contacto producido entre más de dos partículas simultáneamente.
Este modelo requiere de propiedades tales como el coeficiente de rigidez asociados al resorte, coeficiente de amortiguamiento asociado al amortiguador y el coeficiente de fricción de Coulomb, que tienen que estar relacionados con las propiedades de las partículas.
El modelo de esfera blanda tiene una fuerte dependencia de la dinámica de la colisión dada su aproximación del contacto a un sistema de resorteamortiguador [17] (Figura 1).

En la Figura 1 se observa lo que sucede durante una colisión entre las partículas n y m, de donde resultan cantidades de interés tales como el traslape en la dirección normal δ_n y su velocidad relativa V_nm.

La Ecuación 11 describe el traslape δn entre partículas donde D_n y D_m son los diámetros de las partículas, r_n y r_m son las posiciones de las partículas para las partículas n y m, respectivamente.

La Ecuación 12 describe la velocidad relativa en las partículas n y m, donde y son sus velocidades angulares, L_n y L_m están descritas en las Ecuaciones 12a y 12b, respectivamente.
El análisis del contacto entre partículas dada su geometría esférica se hace más simple si se realiza en las direcciones normal y tangencial, por tanto, es necesario encontrar dichas componentes de la velocidad relativa y , siendo esta determinas por las Ecuaciones 13 y 14, respectivamente.

Donde es el vector unitario en la dirección normal.
En la Figura 2 se observa una representación esquemática del contacto entre partículas en la dirección normal y tangencial, respectivamente.

Las Ecuaciones 15 y 16 determinan las componentes normal y tangencial de la fuerza que se genera en el contacto entre la partícula n-esima y m-esima respectivamente, donde representa la fuerza conservativa asociada al resorte y representa la fuerza disipativa asociada al amortiguador en la dirección normal al contacto. Y representa la fuerza conservativa asociada al resorte y representa la fuerza disipativa asociada al amortiguador en la dirección tangencial al contacto.
Donde las componentes conservativas se expresan como la ley de Hooke para resortes, en ambas direcciones, siendo esto expresado en las Ecuaciones 17a y 17b.

Donde k_n y k_t es el coeficiente de rigidez (o de resorte) en la dirección normal y tangencial, respectivamente, con siendo el traslape en dirección tangencial que se produce al momento de iniciar el contacto entre partículas n-esima y m-esima, expresado en la Ecuación 18.

Una vez iniciado el contacto es necesario tener en cuenta su historia, puesto que el traslape tangencial se acumula a medida que la partícula intenta introducirse en la otra, por lo tanto para un tiempo t+Δt se define en la Ecuación 18a.

Debido al desplazamiento de las partículas los planos tangenciales donde se encuentran y , no necesariamente van a coincidir por lo que se hace necesaria una corrección del vector , con el fin de que corresponda al plano tangencial de la nueva posición sustrayéndole la componente que se extiende hacia la dirección normal según la Ecuación 18b.

Las componentes disipativas de la fuerza en la dirección normal y tangencial deben cumplir la ley de Stokes, dado que estas están asociadas a un amortiguador, por tanto las Ecuaciones 19 y 20 describen las componentes de la fuerza en la dirección normal y tangencial para el tiempo t.

Donde η_n y η_t son los coeficientes de amortiguamiento en dirección normal y tangencial, respectivamente, y es la velocidad de deslizamiento (slip) entre las partículas n-esima y m-esima que se determina mediante la Ecuación 21.

Ya determinadas las componentes conservativas y disipativas de la fuerza en la dirección normal y tangencial, la fuerza tangencial debe cumplir con la restricción que impone la ley de fricción de Coulomb (Ecuación 22).

Donde μ es el coeficiente de fricción entre las partículas. Si se alcanza la condición expresada en la Ecuación 22 se establece el deslizamiento y es necesario usar las siguientes condiciones.

Una vez definida la fuerza que se genera en el contacto entre la partícula n-esima y las partículas en su vecindad, es necesario sumar todas estas contribuciones a la fuerza y de este modo determinar la fuerza neta que actúa sobre dicha partícula debida a las interacciones partículapartícula, siendo esto descrito por la Ecuación 23.

Donde N es el número máximo de partículas con la cual entra en contacto la partícula n-esima en un tiempo t del mismo modo se expresa el momento que hace que la partícula n-esima gire sobre su propio eje la cual es determinada por la Ecuación 24.

Donde L_n es la distancia definida por la Ecuación 12a.

Solución, condiciones de frontera e iníciales
La solución de las ecuaciones asociadas a la fase Euleriana (agua residual), se resuelven en estado transiente con un algoritmo predictorcorrector cuyo paso es Δt_Euleriano, teniendo como fundamento el algoritmo Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) propuesto por Patankar y Spalding 1972 para resolver las ecuaciones, incluyendo los términos transientes. Este algoritmo predictor-corrector se desarrolla sobre una malla escalonada donde se almacena los datos de velocidad en las caras de los volúmenes de control y el valor de las propiedades en los centros de los volúmenes de control. La discretización de las ecuaciones eulerianas fue de primer orden, estimando los valores de las propiedades en las caras con el método Up-Wind. El método de integración del modelo Lagrangiano, se basa en el algoritmo de Euler que consisten en truncar la serie de Taylor de la ecuación de movimiento en el término Δt², siendo este un algoritmo de primer orden. Es deseable que el algoritmo de integración sea estable, eficiente y preserve la energía [17], para ello dicho algoritmo requiere de un paso de tiempo pequeño, siendo este la 50ª del tiempo de colisión entre las partículas n-esima y m-esima, Δt_Lagrangiano=min(t_nm^col/50), siendo el tiempo de colisión dado por la Ecuación 25 [18].

Las escalas de tiempo distintas entre los fenómenos que suceden en la fase Euleriana y la Lagrangiana imponen que el algoritmo acoplado sea segregado, es decir, se avanza la solución Lagrangiana hasta que el tiempo acumulado por sus avances igualen el paso de tiempo de la solución Euleriana (∑Δt_Lagrangiano = Δt_Euleriano), alcanzado este punto se pone en marcha la solución Euleriana y se continua nuevamente con las solución Lagrangiana [17].
En la Tabla 2 se relacionan las condiciones de frontera que se aplican para la fase Euleriana.

Las fronteras para el modelo Lagrangiano se toman como pared todas las fronteras de la geometría, para no permitir la salida o entrada de partículas al interior del mismo.

En modelo matemático planteado permitió simular el flujo, y por ende evaluar los parámetros tales como altura del lecho y el proceso de entremezclado y segregación de las partículas para la distribución de tamaño de partículas.

La altura del lecho fue estimada tomando el promedio de la posición vertical de las celdas con una fracción volumétrica de líquido de 0,999 para tres posiciones transversales a 2,5, 6,5 y 12,5cm respectivamente, el resultado de esto se visualiza en la Figura 3 indicando que el lecho del reactor se expande una vez arranca el flujo, empezando a mostrar una altura estable pasados los 240s, luego de trascurrir 80s el lecho alcanza un máximo de 23cm de altura a los 320s. Este comportamiento sugiere que se ha alcanzado un estado estacionario en este parámetro, esto es de gran importancia, debido que la altura del lecho es un parámetro importante para calcular el tiempo en que las dos fases permanecen en contacto dentro del reactor (tiempo de residencia). Es importante observar que este resultado es consecuente con la premisa de que las partículas se mantienen siempre dentro del reactor y no son expulsadas por la salida.

En la Figura 4 se observa la fracción volumétrica de partículas en la dirección axial promediando la extensión transversal del reactor, y como un promedio temporal desde los 240 y 440s rango en el cual la hidrodinámica alcanza un estado estacionario. Dicha figura exhibe una tendencia inicial que se extiende hasta la altura aproximada de 10cm con un valor de la fracción de volumen de partículas que oscila alrededor de 0,21, en segunda instancia se encuentra una tendencia decreciente que va desde los 10cm hasta los 13cm de altura aproximadamente hasta que la fracción de partículas alcanza el valor de 0,16, lo anterior implica que la mayor presencia de partículas se encuentra hasta los 13cm de altura formando una mezcla densa (no diluida); por último, a partir de los 15 y hasta los 23cm aproximadamente hay una disminución progresiva de la fracción volumétrica de partículas hasta llegar a cero, en esta franja se constituye una zona donde el flujo multifase esta diluido, puede observarse en la Figura 5k que a 15cm de altura las partículas están más separadas y finalmente a partir de los 23cm de altura hacia adelante se encuentra la zona libre de partículas.

Lo anteriormente expuesto evidencia que se forman tres zonas, una densa, una diluida y por último una zona libre de partículas, tal como se puede encontrar en las descripciones de los flujos que se presenta en reactores anaeróbicos de flujo ascendente encontrados en la literatura [13].

La formación de las tres zonas siendo la más densa en cuanto a partículas la que se encuentra hacia la base del dominio de simulación, indica que la mayor cantidad de sólidos se encuentran en los primeros 13cm de altura, dando cuenta que la gravedad tiene un efecto muy importante sobre las partículas, como consecuencia de la baja velocidad del flujo ascendente y por lo tanto no genera el arrastre suficiente para vencer dicha fuerza de cuerpo.

La presencia de una fase dispersa con una distribución de tamaños de partículas como la mostrada en la Tabla 1, hacen que las características de expansión del lecho sean más complejas, debido a sus diferencias en las velocidades de asentamiento, esto genera un comportamiento que lleva a las partículas pertenecientes a las diferentes clases a entremezclarse o segregarse [12].

En la Figura 5 se observa que la presencia de las 10 clases de partículas exhiben un comportamiento que va desde el entremezclado hasta la segregación mostrado en el instante de 440s. Esto se debe al comportamiento diferente de las fuerzas de arrastre en cada una de las partículas pertenecientes a cada una de las clases, El comportamiento hacia el entremezclado y a la segregación de una clase respecto de la otra dependerán del tamaño [12].

En la Figura 5a y 5b se observa una tendencia fuerte al entremezclado entre la clase DP01 y la clase DP02, respecto a las clases DP03, DP04, DP05, DP06 y DP07 aún persiste la tendencia al entremezclado entre ellas, sin embargo, esta va disminuyendo a medida que se compara con las clases de mayor tamaño, como puede observarse cuando se compara el grupo de clases anterior con las clases DP08, DP09 y DP10, siendo clara esta tendencia para las clases DP09 y DP10, y se establece rotundamente entre las clases DP01 y DP10.
En la zona central de la distribución de tamaños se encuentran DP03, DP04, DP05, DP06 y DP07 se observa una tendencia clara al entremezclado (Ver Figuras 5c, 5d, 5e, 5f y 5g), pero al analizar de este grupo las clases DP03 y DP07 se observa una tendencia mixta o de transición entre el entremezclado y la segregación. Por último las clases DP08, DP09 y DP10 presentan una tendencia al entremezclado entres ellas (ver Figuras 5h, 5i y 5j), mientras que muestran una tendencia clara a la segregación respecto las clases DP01, DP02 y DP03, mostrando un comportamiento de transición con las clases de la zona central de la distribución. Este comportamiento permite concluir que el flujo establece una clasificación por tamaño de las diferentes clases involucradas.
El entremezclado y la segregación han sido observados experimentalmente encontrando que sucede segregación parcial cuando la relación entre los diámetros de partículas de las clases están entre 1,24-1,58 [19] siendo este el comportamiento observado por DP03 y DP07, DP05 y DP10. Pruden y Epstein en otro experimento encontraron que partículas sólidas con diámetros menores a 1,3 se entremezclan como ocurre para las clases DP04, DP05 y DP06, y para relaciones mayores las fases se segregan incluyendo en este grupo las clases DP03 y DP07. Un estado de completa segregación es posible cuando participan clases cuya relación de tamaños de partículas superen el valor de 2,10 [20] sucediendo esto en las clases extremas DP01, DP02, DP09 y DP10. Lo anteriormente expuesto muestra concordancia entre los resultados experimentales y los resultados calculados por el modelo presentado, debido a que la tendencia a entremezclarse o a segregarse se da en concordancia con la descripción y la relación de diámetros reportados por los autores mencionados (Tabla 3).

En la Figura 6 se muestran las velocidades promedio en dirección axial de las clases promediadas en el tiempo entre 240s y 440s, lo que explica desde el comportamiento de las clases de partículas el proceso de entremezclado y de segregación siendo éste consecuencia de la diferencia de sus velocidades promedios en la dirección axial. De la figura se establecen tres grupos de clases; el primero formado por DP01 y DP02 siendo estas clases las que más afectan la fuerza de arrastre lo que lleva a que su velocidad de ascenso sea mayor, por lo que tienden a ubicarse en la zona superior del lecho.
El segundo formado por DP03, DP04, DP05, DP06 y DP07 cuya variación en la velocidad promedio no es tan marcada respecto del grupo anterior coincidiendo en la ubicación central del lecho en una tendencia de entremezclado. Sin embargo, los dos extremos DP03 y DP07 tienen una diferencia marcada en su velocidad por lo que éstas dos clases se encuentran en la condición mixta (entremezclado - segregación).
Por último el grupo formado por las clases DP08, DP09 y DP10 que muestran una velocidad promedio negativa lo que implica que estas clases se desplacen hacia el fondo del lecho, explicando el porqué se establece la segregación entre las clases de los extremos DP01, DP02, DP09 y DP10. Los promedios de velocidad por clases y temporal obtenidos por el modelo Lagrangiano avalan teóricamente los resultado obtenidos experimentalmente por Pruden, Epstein, Al - Dibouni y Garside. En las Figuras 5a, 5b, 5g, 5h, 5i y 5j se observa una fuerte presencia de las clases DP01, DP02, DP07, DP08, DP09 y DP10 hacia el sector del dominio entre 7cm y 14cm en dirección X," es importante para que esta ultima parte recobre el sentido. siendo esto atribuible al promedio de la velocidad en la dirección transversal ver Figura 7.
En la Figura 7 se observa la velocidad promedio en dirección transversal de las clases promediadas en el tiempo entre 240 y 440s, siendo este resultado la base para explicar la mayor presencia de partículas de las clases DP01, DP02, DP07, DP08, DP09 y DP10 en el sector comprendido entre 7cm y 14cm dirección X, ya que su velocidad promedio es positiva, además siendo las clases DP09 y DP10 las que más rápido se mueven hacia esta dirección son las que mayor presencia muestran en las Figuras 5i y 5j.

El modelo permitió obtener una información más amplia y por ende darle una mejor explicación al fenómeno de entremezclado y segregación que ocurre en flujos dispersos con partículas con diferente diámetro, lo cual se debe a la respuesta diferente en relación al arrastre que tienen las partículas que se someten al mismo campo de flujo, siendo esto lo que consolida una velocidad promedio diferencial de cada una de las clases que participan en el flujo dándose así bien sea el entremezclado o la segregación.
El modelo desarrollado permitió determinar la formación de tres zonas dentro de la geometría, la zona inferior que se extiende hasta los 13cm se configura una zona densa, debido a su alto contenido de partículas. A partir de los 13cm hasta los 23cm, que es la altura que alcanza el lecho, se configura una zona diluida y por último una zona libre de partículas, la robustez del modelo Euleriano-Lagrangiano potenciado con el DEM permite reproducir bien dichas condiciones de flujo dado que se basa en el comportamiento discreto de las partículas.

Al grupo de desarrolladores de Multiphase Flow with Interphase Exchange MFIX del Departamento de Energía de EE.UU, Al grupo de investigación TAYEA y a la Universidad Nacional de Medellín.

[1] Pfleger D, Gomes S, Gilbert N, Wagner H-G. Hydrodynamic simulations of laboratory scale bubble columns fundamental studies of the Eulerian-Eulerian modelling approach. Chem. Eng. Sci. 1999;54(21):5091-9.

[2] Jaworski Z, Pianko-Oprych P. Two-phase laminar flow simulations in a kenics static mixer. Chem. Eng. Res. Des. 2002;80(8):910-6.

[3] Ataki A, Bart H-J. The use of the VOF-Model to study the wetting of solid surfaces. Chem. Eng. Technol. 2004;27(10):1109-14.

[4] Zhonghua W, Mujumdar AS. CFD modeling of the gas-particle flow behavior in spouted beds. Powder Technol. 2008;183(2):260-72.

[5] Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. Géotechnique. 1979;29(1):47-65.

[6] Stein E, Borst R de, Hughes TJR. Encyclopedia of computational mechanics. Chichester, West Sussex: John Wiley; 2004.

[7] Helland E, Occelli R, Tadrist L. Computational study of fluctuating motions and cluster structures in gas-particle flows. Int. J. Multiphase Flow. 2002;28(2):199-223.

[8] Zhang MH, Chu KW, Wei F, Yu AB. A CFD- DEM study of the cluster behavior in riser and downer reactors. Powder Technol. 2008;184(2):151-65.

[9] Grof Z, Cook J, Lawrence CJ, Štĕpánek F. The interaction between small clusters of cohesive particles and laminar flow: Coupled DEM/CFD approach. J. Pet. Sci. Eng. 2009;66(1-2):24-32.

[10] Pirker S, Kahrimanovic D, Kloss C, Popoff B, Braun M. Simulating coarse particle conveying by a set of Eulerian, Lagrangian and hybrid particle models. Powder Technol. 2010;204(2-3):203-13.

[11] Fries L, Antonyuk S, Heinrich S, Palzer S. DEM-CFD modeling of a fluidized bed spray granulator. Chem. Eng. Sci. 2011;66(11):2340-55.

[12] Reddy RK, Joshi JB. CFD modeling of solid-liquid fluidized beds of mono and binary particle mixtures. Chem. Eng. Sci. 2009;64(16):3641-58.

[13] Saravanan V, Sreekrishnan TR. Modelling anaerobic biofilm reactors-A review. J Environ Manage. 2006;81(1):1-18.

[14] Batstone D. Variation of bulk properties of anaerobic granules with wastewater type. Water Res. 2001;35(7):1723-9.

[15] Nicolella C, van Loosdrecht MCM, Heijnen SJ. Particle-based biofilm reactor technology. Trends Biotechnol. 2000;18(7):312-20.

[16] Crowe C. Multiphase flow handbook. Boca Raton FL: CRC;Taylor & Francis; 2006.

[17] Garg R, Galvin J, Li T, Pannala S. Opensource MFIX-DEM software for gas-solids flows: Part I-Verification studies. Powder Technol. 2012;220:122-37.

[18] Silbert L, Ertaş D, Grest G, Halsey T, Levine D, Plimpton S. Granular flow down an inclined plane: Bagnold scaling and rheology. Phys. Rev. E [Internet]. 2001;64(5).

[19] Pruden BB, Epstein N. Stratification by size in particulate fluidization and in hindered settling. Chem. Eng. Sci. 1964;19(9):696-700.

[20] Al-Dibouni M R, Garside J. Particle classification and intermixing in liquid fluidized bed Flow. T I Chem Eng-Lond. 1979;57:94-103.