 |
ISSN electrónico: 2145-9371 Fecha de recepción: 29 de febrero de 2008 |
Flujo de Poiseuille y la cavidad con pared móvil calculado usando el método de la ecuación de lattice Boltzmann
Poiseuille flow and the lid-driven cavity calculate using the Lattice Boltzmann equation method
Elkín G. Flórez S.*, Ildefonso Cuesta**, Clara Salueña**
* Magíster en Ingeniería Mecánica, candidato a Doctor de la Universidad Rovira i Virgili, Tarragona (España). Grupo de Investigación ECOMMFIT, profesor asistente, Programa de Ingeniería Mecánica, Universidad de Pamplona. eflorez@unipamplona.edu.co
**Doctor en Ingeniería Química, Grupo de Investigación ECOMMFIT, Universidad Rovira i Virgili. ildefonso.cuesta@urv.cat, clara.saluena@urv.catCorrespondencia: Avenida Pisos Catalans N° 24, Tarragona (España).
Resumen
El objetivo de este artículo es presentar los resultados de la aplicación del Método de Lattice Boltzmann (LBM) como una herramienta de solución en la dinámica computacional de fluidos. Después de una corta revisión de la teoría básica y utilizando el modelo bidimensional de 9 velocidades (D2Q9), el flujo de Poiseuille es simulado y los resultados son comparados con la solución analítica existente. También, es modelada la cavidad con pared móvil (Lid-driven) y los resultados obtenidos validados con datos existentes (Guía et al.). Las condiciones de frontera para pared estática y pared móvil son revisadas en el primer y segundo modelo, respectivamente. Los resultados indican la eficiencia del LBM para simular flujos de fluido incompresibles y laminares. También, que como efecto de incrementar el número de puntos en el lattice, mejora la convergencia computacional y reduce las oscilaciones espaciales de la solución cerca de puntos geométricamente singulares en el flujo.
Palabras claves: lattice Boltzman, flujo de Poiseuille, cavidad de tapa móvil.
Abstract
The aim of this article is to present the results of the lattice Boltzmann method (LBM) application as computational fluid dynamics solvers. After of short review of the basic theory and using the two-dimensional model with 9 velocities (D2Q9), the Poiseuille flow is modelled and validated the results with the analytical solutions. Also, the Lid-driven cavity is modelled and validated the results with existing data (Guía et al.). The boundary condition for static wall and moving wall are revised on the first and second model respectively. The results indicate the efficiency of LBM to simulate incompressible and laminar fluid flow. Also, that the effects of increment in the number of the lattice points, improve the computational convergence and reduce spatial oscillations of solution near geometrically singular points in the flow.
Key words: Lattice Boltzmann, Poiseuille flow, Lid Driven Cavity.
INTRODUCCIÓN
El movimiento de un fluido, visto como un continuo, está gobernado por las ecuaciones de continuidad
y las de Navier-Stokes (NE)
Donde u es la velocidad del fluido, p la densidad, P la presión, v la viscosidad cinemática y ζ viscosidad global (bulk viscosity). Las anteriores ecuaciones son de gran interés tanto por su complejidad matemática para los profesionales del ramo como por su amplio rango de aplicación para los ingenieros. La solución exacta o analítica de las mismas en muchos casos complejos no ha sido alcanzada. Por lo tanto, encontrar soluciones a problemas que involucran las ecuaciones 1 y 2 requiere de desarrollos experimentales o soluciones numéricas.
Desde que apareció la simulación numérica, esta ha ido reemplazando paulatinamente los métodos experimentales usados durante mucho tiempo en áreas bien probadas como en la ingeniería naval y aeroespacial. Ahora la dinámica de fluidos computacional (CFD, por sus siglas en inglés) se aplica a prácticamente cualquier área de la ingeniería como la mecánica, química de procesos, tecnología, etc.
La dinámica computacional de fluidos consiste en transformar el sistema de ecuaciones diferenciales parciales 1 y 2 en un sistema algebraico no lineal.
Así, con la transformación se hace discreta la variable tiempo por medio de métodos explícitos (fáciles de implementar, pero computacionalmente lentos) o implícitos (difíciles de implementar, rápidos y eficientes) tal y como se hacen discretas las variables espaciales por medio de cualquiera de los siguientes métodos: diferencias finitas, volúmenes finitos, elementos finitos, etc.
A diferencia de métodos utilizados en CFD tradicionales, el Método de Lattice Boltzman (LBM, por sus siglas en inglés) no resuelve las NE, sino que el modelo es construido sobre un espacio lattice (o reticulado) que contiene partículas de fluido. Cada una de estas partículas está dotada de un conjunto de velocidades discretas para moverse de un nodo a otro sobre la malla o dominio del problema. Las partículas son redistribuidas sobre cada nodo de acuerdo a un conjunto de normas que cubren el proceso dinámico (propagación y colisión) entre partículas [1].Los modelos lattice Boltzmann requieren de la total especificación de:
- (a) Un espacio lattice discreto donde reposan las partículas del fluido.
- (b) Un conjunto de velocidades discretas (que describen cómo pasan de un nodo a otro nodo vecino) para representar la advección de las partículas.
- (c) Un conjunto de reglas para la redistribución de las partículas que residen en un nodo para imitar el proceso de colisión de un fluido real.
1. MODELO BHATNAGAR-GROSS-KROOK EN EL MÉTODO LATTICCE BOLTZMANN
La ecuación de lattice Boltzmann está dada por [2-3-4]
Donde es la función de colisión, que se obtiene al establecer un promedio sobre el conjunto representativo de partículas y despreciar sus correlaciones.
El término de colisión debe satisfacer las leyes de conservación y ser compatible con la simetría del modelo. Este término puede simplificarse bastante utilizando la simple aproximación de Bhatnagar-Gross-Krook (BGK)[5], con un solo tiempo de relajación τ (SRT, por sus siglas en inglés), que es el método más sencillo utilizado para resolver la ecuación (3), y consiste en reemplazar el término de colisión por un término relajado, de la función f a la función de equilibrio de Maxwell-Boltzmann [2]
1.1. Lattice Boltzmann con el modelo de simple tiempo de relajación
El método LBM resuelve la ecuación cinética microscópica para la función de distribución de partículas f(x,υ,t), donde x y υ son la posición y velocidad de la partícula, respectivamente, en el espacio fase (x, υ) y tiempo t, donde las cantidades macroscópicas como velocidad y densidad se obtienen por medio de los momentos de f(x,υ,t). Si se utiliza el SRT con la aproximación BGK, la ecuación de lattice Boltzmann en el espacio de fase discreta, está dada por:
donde fi(x,t) y fi(eq)(x,t) son la función de distribución y la función de distribución de equilibrio de la i-ésima partícula, respectivamente, ei es el vector velocidad discreta y t es el tiempo de relajación el cual está relacionado con la viscosidad cinemática del fluido por medio de [2-3-4-6]
En el presente trabajo se utilizó el modelo D2Q9, mostrado en la figura 1, para simular el flujo de Poiseuille y la cavidad con pared móvil. Con este modelo y para flujos isotérmicos e incompresibles, la función de distribución de equilibrio viene dada por
donde wi es el factor de peso y u es la velocidad del fluido. Las velocidades discretas para este modelo son:
donde c=Δx/Δt=Δy/Δt y Δt, son la constante lattice y el tamaño del paso de tiempo, respectivamente. El valor de los factores de peso para cada velocidad discreta esta dado por
La densidad y la velocidad del flujo pueden ser obtenidas por medio de la integración de los momentos
La evolución del modelo lattice BGK consiste en dos pasos: propagación y colisión [7][8].
1.2. Implementación de condiciones de frontera
En los modelos LBM como en cualquier simulación de fluidos se requiere implementar las condiciones de frontera que impone el flujo [3][4]. En el presente trabajo se utilizaron tres tipos de condiciones: periódicas; para la entrada y salida de flujo en el problema de Poiseuille, condiciones de pared móvil en el problema de la cavidad con tapa-móvil y para ambos problemas se utilizo el bounce-back para las fronteras de pared fija.
1.2.1. Condiciones periódicas
Las condiciones periódicas son la forma más simple de aplicar condiciones de frontera. Consiste en mantener las condiciones (velocidad, presión, temperatura, etc.), impuestas a un extremo del dominio, iguales a las de su lado o frontera opuesta. Estas condiciones no tienen en cuenta cualquier perturbación que el flujo pueda presentar en sus fronteras; por lo tanto, son solo adecuadas para aplicar en fenómenos físicos donde los efectos de las superficies no jueguen un papel importante y donde no existan perturbaciones del flujo cerca de las fronteras.
Para el flujo de Poiseuille, se utilizó este tipo de condiciones donde se mantienen iguales la función de distribución, para los nodos entrada y nodos salida del flujo, al igual que para los nodos de la pared inferior y superior del canal (el perfil de velocidad y presión del flujo es simétrico en la dirección vertical). Lo anterior se escribe en forma explícita para un canal de longitud L y altura H
1.2.2. Condiciones de pared fija
Para tratar el paso de advección en presencia de fronteras de pared fija se utiliza la aproximación lineal propuesta por Bouzidi en [9], conocida como Bounce-back, y ubicando la pared en la mitad del camino entre el nodo sólido y el nodo fluido. Supongamos que rl es un nodo fluido tal que rl + ei es un nodo sólido. Llamando a ei la velocidad invertida de ei (ei = -ei), tenemos que
En el lado derecho de la ecuación (13), la f c es tomada después de la colisión y antes de la propagación La f ( ,t + 1) del lado izquierdo de la ecuación se usaría en valores después de la colisión y después de la propagación, la cual corresponde a un paso de tiempo LBM completo. Lo anterior indica que la dinámica de los nodos sólidos y los nodos fluido, vecinos de la pared, difieren únicamente en el paso de propagación [9]
1.2.3. Condiciones de pared móvil
Para fronteras con velocidad diferente a cero, se utiliza un argumento de primer orden que introduce el aporte en cantidad de movimiento de la pared móvil a las partículas fluido. Con referencia a la figura 2, se asume la velocidad macroscópica del fluido simulada bajo el esquema LBM como V, y que varía linealmente a lo largo de DA, y conocido el valor de la velocidad de la pared (V0) en B, con el punto C como origen para la coordenada y se obtiene
Para el primer orden en V los valores de equilibrio de las fi están dadas por
donde los coeficientes ai dependen del conjunto de velocidades discretas del modelo LBM utilizado. Para el caso D2Q9 los valores son:
Los valores de fi son el mismo dado a los factores de peso en la ecuación 7.
Para calcular los correspondientes valores de fi' en la ecuación 12 se debe agregar al lado derecho de la ecuación el término, debido al movimiento de la frontera
Las deducciones anteriores demuestran que las fronteras móviles pueden ser simuladas aplicando un simple Bounce-back más un término adicional de cantidad de movimiento que agrega la pared móvil a las partículas fluido.
1.3. Implementación computacional
El desarrollo del presente trabajo comenzó con una amplia revisión en las investigaciones realizadas por otros autores sobre la implementación de LBM en el flujo de Poiseuille [3][10][11][12], en la cavidad con pared móvil [6][13][7][14], y en investigaciones que incluyen ambos problemas al tiempo [15][16][18]. Con el conocimiento de la teoría existente y la naturaleza del método se manipularon los principales parámetros que influyen en la simulación, evitando violar los principios físicos existentes.
Se desarrolló un código bidimensional, utilizando el modelo LBM-BGK D2Q9. El código fue escrito en el lenguaje de programación FORTRAN 90. La compilación, ejecución y gráficos se realizaron utilizando los recursos existentes en el departamento de ingeniería mecánica de la Universidad Rovira i Virgili:
- Recurso de hardware: Beowulf de 24 procesadores AMD opteron 248 dual core (64 bits) y 7 procesadores Intel de 3 Ghz, con 3 terabytes de disco, unido con un gigaethernet en un ambiente Linux.
Los resultados son visualizados utilizando:
- Recurso software: Gnuplot y OpenDX
En la figura 3 se muestra el diagrama de flujo utilizado para el desarrollo del código computacional con el cual se obtuvieron los resultados. Existen diferentes plataformas tales como C, C++, Fortran, este último utilizado en la presente investigación para la implementación de estos tipos de modelos dentro de un código computacional [12]. Finalmente, se puede decir que los códigos se desarrollan en función de las necesidades del investigador y de la aplicación final del mismo.
2. RESULTADOS NUMÉRICOS
Para la simulación se utilizaron diferentes valores de tiempo de relajación τ, sin violar lo establecido en [9]. En este artículo solo se presentan los resultados que obtuvieron mayor exactitud. El número de iteraciones realizadas en cada problema es función de la condición de estado estable de la simulación, la cual fue establecida en
donde i en la sumatoria indica todos los nodos de la malla. Como se muestra en la figura 4, para que el flujo de Poiseuille alcance el estado estacionario toma unas 20.000 iteraciones aproximadamente, mientras que en la cavidad con tapa móvil, unas 150.000 iteraciones, aproximadamente.
Flujo de Poiseuille
Son pocas las soluciones analíticas que existen para las ecuaciones de Navier-Stokes. Una de ellas es el flujo plano de Poiseuille en un canal con ancho 2H donde el flujo es estable (∂ /∂ t = 0), en la dirección x, u = (u,υ) = (u(y),0), con presión constante y sin variación en la dirección x (∂ /∂ x=0). El flujo es inducido por una fuerza constante K = Kex. Resolviendo las ecuaciones de N-S con las condiciones antes mencionadas, se llega a la ecuación diferencial para u(y)
En las paredes del canal se imponen condiciones de frontera de no-deslizamiento (u(y) = 0 en y = -H e y=H). La solución analítica es una parábola de la forma
Para la solución LBM, se induce un desbalance en la primera columna de nodos del canal, que permite simular el cambio de presión entre la entrada y la salida del flujo dentro del canal. Se aplican condiciones periódicas a la entrada y la salida del canal. Las paredes inferior y superior del canal son ubicadas en los nodos y=1 e y=Ny', (para todos los nodos x), respectivamente. El tamaño del dominio es Nx × Ny' = 100 × 20, incluidos los nodos frontera. La condición inicial para la simulación es u(x,0)=0 para todos los nodos internos o nodos fluido. El valor utilizado para el tiempo de relajación es τ=0.54, lo cual corresponde a una viscosidad cinemática υ=0.0135 en unidades lattice.
La gráfica del lado izquierdo de la figura 5 muestra los resultados obtenidos para la velocidad u(x), con el LBM y la respectiva comparación con la solución analítica dada por la ecuación 19. También se muestra en el lado derecho de la figura 5 que la velocidad en la dirección vertical se puede considerar cero, ya que la magnitud numérica obtenida es de un orden menor a 10-6, pues cumple con la condición teórica del flujo de Poiseuille para la velocidad en esta dirección.
Cavidad con tapa móvil
Para validar los resultados obtenidos en la simulación de la cavidad con pared móvil se utiliza el trabajo desarrollado por Ghia [17]. La cavidad cuadrada con pared móvil es una de las más comunes y buenas dentro de los estándares existentes para validar resultados de un modelo o método numérico en la dinámica de fluidos computacional.
En el presente trabajo se comenzó la simulación con un tamaño mínimo de malla 64×64 y uno máximo de 500×500, con tres diferentes números de Reynolds 100, 1000 y 3200, como lo muestra la tabla 1. Se investigó sobre la influencia del tamaño de la malla N×Ny', la velocidad de la pared upared', y el parámetro de relajación τ en la exactitud de los resultados. La tabla 1 también muestra el valor de la viscosidad cinemática v del fluido simulado, calculada a partir de la ecuación 5.
En las figuras 6,7 y 8 se muestran los resultados obtenidos, utilizando diferente tamaño de malla y número de Reynolds comparados con el tradicional bencjmark para este tipo de configuración desarrollado por Guía, et al.El número de Reynolds es calculado por
La ecuación anterior permite determinar que para no violar los principios físicos, donde la velocidad de la pared ha de ser máximo un 10% de la velocidad del sonido [2, 3, 4] y el parámetro de relajación condiciona la viscosidad cinemática por medio de la ecuación (5), se deduce que en la utilización del método (LBM) se requiere aumentar el tamaño de la malla si se requiere aumentar el número de Reynolds en busca de flujos en régimen de transición y/o turbulentos.
Los resultados obtenidos validan el buen funcionamiento del código desarrollado por los autores en el lenguaje Fortran y dan luz a desarrollar nuevas aplicaciones en la dinámica de fluidos utilizando los LBM.
CONCLUSIONES
Se revisó el comportamiento de dos tipos de flujos laminares, simulados por medio del modelo LBM-BGK con simple tiempo de relajación. Los resultados obtenidos permiten observar el buen comportamiento del código desarrollado.
Para el flujo de Poiseuille, el flujo es inducido mediante un gradiente de presión entre la entrada y la salida del canal. El perfil de velocidad obtenido a lo largo del canal es descrito por una parábola con máxima velocidad en su línea central y con velocidad en la dirección vertical aproximadamente igual a cero, que se valida con los resultados teóricos.
En la cavidad con pared móvil se observó que para un mismo valor en el parámetro de relajación (τ), la exactitud de la simulación recae en el tamaño de la malla. A mayor numero de nodos que sean tomados en el dominio, el error disminuye. Lo anterior cae en la cuenta que la malla no puede ser tan grande porque un tamaño tal viola el máximo valor que debe tener el parámetro de relajación (0.5).
Las condiciones de frontera utilizadas, el bounce-back y el bounce-back con inclusión de cantidad de movimiento, demostraron ser un buen modelo a la hora de simular flujos laminares con pared quieta y pared móvil, respectivamente.
La comparación de los resultados en la cavidad de tapa móvil, respecto al benchmark (Ghia, et al.) es relacionado con la exactitud que se genera con los resultados finales y no se han comparado ni el tamaño de malla ni la resolución de la misma, ya que los principios físicos utilizados son diferentes y solo en un artículo posterior se comentará la eficiencia del método con respecto a los CFD tradicionales.
REFERENCIAS
[1] M. BOUZIDI, D. D´HUMIERES, P. LALLEMAND, L. Luo, "Lattice Boltzmann Equation on a Two-Dimensional Rectangular Grid", Journal of Computational Physics, 172, pp.704-717, 2001.
[2] X. HE, L-S LUO, "Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation", Phys. Rev. E 56 (6), pp. 6811-17, 1997.
[3] DIETER A., Wolf-Gladrow, "Lattice-Gas Cellular Automata and Lattice Boltzmann Models", Springer, 2000, pp. 40-65.
[4] SAURO SUCCI, The lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford, 2001, pp. 64- 93.
[5] P.L. BHATNAGAR, E.P. GROSs, and M. KROOK. "A model for collision processes in gases. I Small amplitude processes in charged and neutral one-component system", Phys. Rev., 94 pp. 511-25. 1954
[6] HOU, S., ZOU, Q., CHEN, S., DOOLEN, G., and COGLEY, AC., "Simulation ofcavity flow by the lattice Boltzmann method". Journal of Computational Physics, 118, pp. 329:347, 1995.
[7] J.S. WU and Y. L. SHAO, "Simulation of lid-driven cavity flows by parallel lattice Boltzmann method using multi-relaxation-time scheme", Int. J. Numer Meth. Fluids, 46, pp. 921-937, 2004.
[8] DAZHI, Yu, RENWEI, Mei, LI-SHI, Luo and WEI SHYY. "Viscous flow computations with the method of lattice Boltzmann equation". Progress in Aerospace Sciences 39, pp. 329-367, 2003.
[9] QUIAN, Y. H., D'HUMIÈRES D., LALLEMAND, P., "Lattice BGK models for Navier- Stokes equation", Europhysics Lett 17 (6), pp. 479-84, 1992.
[10] X. HE, Q. ZOU, L-S LUO and M. DEMBO, "Analytical solution of dimple flows and analysis of nonslip boundary conditions for the lattice Boltzmann BGK model", Journal of Statistical Physics, vol 87, Nos. 1/2, 1997.
[11] XIAOYI, H., and LI-SHI, LUO, "Lattice Boltzmann Model for the incompressible Navier- Stokes equation", Journal of Statistical Physics vol 88, pp. 927-935, 1997
[12] S. T. ENGLER, "Benchmarking the 2D lattice Boltzmann BGK model", Short communication. Amsterdam Center for Computacional Science, Amsterdam, The Netherlands.
[13] W. MILLER, "Flow in the driven cavity calculated by the lattice Boltzmann method", Phys Rev E; 51(4), pp. 3659-3669, 1995.
[14] D. V. PATIL, K. N. LAKSHMISHA and B. ROGG, "Lattice Boltzmann simulation of lid-driven flow in deep cavities", Computer & fluids [35], pp. 1116-1125, 2006
[15] S. CHEN, D. MARTÍNEZ and R. MEI, "On boundary conditions in lattice Boltzmann methods", Phys. Fluids 8 (9), pp. 2527-2536, 1996
[16] Y. T. CHEW, C. SHU and Y. PENG, "On implementation of finite volume lattice Boltzmann method". Journal of Statistical Physics, vol. 107, n.os 1/2, 2002.
[17] GHIA, U., GHIA, K. N. and SHIN, C. T., "High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method", J Comput Phys, 48, pp. 387-411, 1982.
[18] http://www.ccrl-nece.de/lba
|
Revista Ingeniería y Desarrollo |