Computational simulation of fluid flow by using element based finite volume method on unstructured grids

Carlos D. Araujo, José A. Rincón, Gilberto I. Materano y Alejandro A. Colman 

Laboratorio de Simulación Computacional, Escuela de Mecánica, Facultad de Ingeniería, Universidad del Zulia, Apartado 526. Maracaibo, Venezuela

Abstract 

The results obtained in the development of the computational program “MEFVOC”, to simulate laminar incompressible and steady flows, using Element based Finite Volume Methods on unstructured and bidimensional grids, is presented in this paper. There are references to problems for testing and evaluating numerical techniques related with this method, but it is not common to find results on unstructured grids like the ones reached in this work. The implementacion makes easier the extention to complex geometries. The problems resolved demonstrate concordance with other authors, using structured grid. In fact, it can be concluded that the computational program developed, succesfully resolves the problems proposed under the described conditions. 

Key words: Control volume finite element, unstructured grid, MEFVOC.

Desarrollo de un simulador de flujo de fluidos usando el método volúmenes finitos basado en elementos para mallas no estructuradas

Resumen 

En este trabajo, se presentan los resultados obtenidos luego del desarrollo de un simulador de problemas de flujo de fluidos en régimen laminar, estado estable, usando el método de volúmenes finitos basado en elementos, para dominios bidimensionales, llamado MEFVOC. De la revisión de referencias se observa con reiteración el uso de problemas típicos que permiten comprobar el buen funcionamiento del programa de cálculo, en dominios discretizados de forma estructurada. Por lo tanto, el aporte de este trabajo consiste en implementar computacionalmente el método, con la capacidad de reproducir la misma solución, pero al emplear una malla no estructurada, lo que facilita la extensión hacia geometrías complejas. Por el hecho de que los resultados obtenidos son congruentes con trabajos de otros autores en el área, al usar programas semejantes, está claro que el programa resuelve exitosamente los problemas planteados bajo las condiciones descritas. 

Palabras clave: Volumen finito basado en elementos, mallas no estructuradas, MEFVOC.

Recibido: 23-05-06  Revisada: 19-11-07

Introducción 

Una circunstancia muy especial en la solución de problemas de fenómenos de transporte, es el hecho de poder expresar todas las ecuaciones gobernantes de una forma general, lo que facilita el desarrollo de procedimientos comunes de solución [1, 2]. Para flujo permanente ésta puede ser escrita como:

donde F es una variable escalar dependiente, n representa el vector velocidad, r la densidad y G es el coeficiente difusivo correspondiente. El valor de es un término fuente por unidad de volumen.

Cuando el transporte de la variable escalar estudiada es debido a un fenómeno difusivo, entonces, la ecuación (1) se simplifica debido a la ausencia del término rnF, que representa el transporte convectivo [3]. 

Por otra parte, la solución de problemas de flujo de fluidos requiere resolver la ecuación de conservación de cantidad de movimiento. En este caso, el término fuente agrega al cálculo una variable adicional, la presión, que no es contenida de forma explícita en ninguna ecuación de conservación. Una alternativa para obtener el conjunto completo de ecuaciones es combinar la ecuación de cantidad de movimiento y el principio de conservación de masa mostrado en la ecuación (2), a través de procedimientos especiales de solución [4, 5].

Referente al método numérico, se utiliza volúmenes finitos basado en elementos, aplicado a problemas de difusión y flujo de fluidos en régimen laminar y en estado estable [1, 6]. En esta técnica, el dominio es dividido es elementos triangulares con un nodo en cada vértice, a su vez, volúmenes de control se forman aplicando el método de la mediatriz, que consiste en unir el centroide de cada triangulo con los puntos de mitad de lado de las aristas [7]. Las expresiones algebraicas son obtenidas mediante la integración de la ecuación (1) a lo largo de las fronteras de cada volumen de control, mostradas como líneas segmentadas en la Figura 1.

Figura 1. Triángulo interno típico [1].

En este procedimiento, la densidad y el coeficiente difusivo se suponen constantes en cada elemento, además, el término fuente solo admite una dependencia lineal con la variable dependiente. Por otra parte, el flujo convectivo difusivo es calculado suponiendo que en cualquier elemento, éste varía exponencialmente en dirección de la velocidad promedio y linealmente en la dirección normal. Si la situación es de difusión pura, la aproximación es lineal. El sistema de ecuaciones resultantes se obtiene mediante el ensamble de las matrices elementales locales y la especificación de las condiciones de borde apropiadas, que pueden ser del tipo Dirichlet o Neuman. Finalmente, el sistema es resuelto usando un método directo, que calcula la inversa de la matriz mediante una descomposición LU. 

Por otra parte, para el cálculo del campo de velocidades se sigue un procedimiento similar al empleado en la ecuación de conservación general para una cantidad F cualquiera, sólo que, el algoritmo “SIMPLER” se usa para relacionar la ecuación de cantidad de movimiento y conservación de masa formando un así un conjunto completo [5, 8, 9]. 

El desarrollo del trabajo se usan algunos parámetros adimensionales necesarios, tales como: el número de Reynolds (Re), el número de Prandtl (Pr), y el número de Grashof (Gr), los cuales pueden ser obtenidos mediante las siguientes expresiones [10]:

donde: L, g, b, Cp, k, s, Tc, y Tf son respectivamente, la longitud característica, la aceleración gravitacional, el coeficiente de expansión volumétrica, el calor específico, la conductividad térmica del fluido, la viscosidad dinámica y niveles de temperatura de referencia, donde Tc es mayor Tf.

Parte Experimental 

A fin de comprobar el correcto funcionamiento del programa desarrollado MEFVOC, se utilizan casos de estudios típicos documentados en la literatura, los cuales han servido a los investigadores del área para probar y evaluar, una y otra vez, técnicas numéricas. La estrategia empleada en este trabajo, no escapa a estas fases típicas de desarrollo. Por lo tanto, en primer lugar se presenta una situación simple de difusión pura, para luego incorporar otros aspectos que le den generalidad al programa. El segundo y tercer caso, se refiere a problemas de campo de flujo laminar en una cavidad cuadrada, motorizado por convección forzada y natural. En la Tabla 1 se puede observar un resumen de la data requerida en la solución de los problemas propuestos.

Tabla 1

Valores de coeficiente de difusión y fuente para cada variable resuelta en cada problema

Es importante mencionar que, aunque las referencias [1, 4, 6, 9-15] presentan formalmente el método empleado en este trabajo, este no es aplicado a casos de prueba en mallas no estructuradas. Por tal motivo, en todos los problemas propuestos, los valores de referencia son para mallas estructuradas. 

Resultados

Problema de difusión 

El problema consiste en obtener la distribución del campo de velocidad (w) cuando el flujo se encuentra completamente desarrollado. Esta situación es resuelta mediante el balance entre las fuerzas de presión y los esfuerzos viscosos causados por los gradientes existentes en la dirección perpendicular al flujo, lo que matemáticamente conduce a un caso particular de la ecuación (1), que es cuando el término convectivo es nulo [3]. El dominio de cálculo es un ducto cuadrado de lado L, con las condiciones de borde señaladas de valor prescrito mostrados en la Figura 2. El flujo es motorizado por un término fuente constante cuyo valor asignado es 1.

Figura 2. Ducto de sección cuadrada.

La discretización se realiza usando elementos triangulares dispuestos de forma no estructurada, como se muestra en la Figura 3.

Figura 3. Discretización del ducto cuadrado.

La malla fue creada utilizando un programa auxiliar, para luego, ser exportada a un archivo que almacena la información de los 283 nodos (16 en cada pared) y 504 elementos empleados. En la Figura 4 se puede observar que la distribución de velocidades es simétrica tanto al eje vertical como al eje horizontal y alcanza un valor máximo hacia el centro del ducto debido a la presencia del término fuente.

Figura 4. Resultados gráficos para el caso 1.

A fin de garantizar el correcto funcionamiento del programa desarrollado, se extrajeron los valores de velocidad sobre la línea de simetría horizontal A-A señalada en la Figura 3, y se compararon con los valores obtenidos con el programa CONDUCT [3], que resuelve las ecuaciones usando el método de volumen finito. Esto puede observarse en la Figura 5.

Figura 5. Comparación entre el programa MEFVOC y CONDUCT.

Nótese, que existe correspondencia entre el perfil de velocidad predicho por el CONDUCT y el alcanzado por MEFVOC, tanto en forma como en magnitud. Por lo tanto, este caso permite constatar el correcto funcionamiento de las rutinas de cálculo de los coeficientes de la matriz local y el ensamble de las ecuaciones en la matriz de conductividad y el vector de carga global. El valor máximo de velocidad adimensional obtenido analíticamente es de 0,0737 [3] y el arrojado por el programa es 0,0735.

Casos de flujo de fluidos 

Los dos problemas descritos a continuación, comparten la misma configuración geométrica, a saber, una cavidad cuadrada de lado unitario. Por tal motivo, se emplea la misma discretización en ambos casos, la cual puede ser observada en la Figura 6. La malla empleada consta de 245 nodos y 433 elementos triangulares dispuestos de forma no estructurada, con una densidad de nodos mayor hacia las esquinas inferiores de la cavidad, a fin de capturar el fenómeno que tienen lugar en esa zona, a medida que el fenómeno convectivo se hace importante.

Figura 6. Discretización del dominio de cálculo.

Adicionalmente, se muestran los volúmenes de control que lucen como polígonos en el dominio de cálculo, sobre los cuales, se resuelven las ecuaciones de conservación que rigen el fenómeno estudiado. 

Convección forzada en una cavidad cuadrada 

En este caso, se analiza un problema en donde el movimiento del fluido está motorizado por una pared con libertad de desplazarse. Un esquema de esta situación se muestra en la Figura 7, en donde tres paredes permanecen inmóviles y la cuarta (tope), se mueve con una velocidad horizontal U. El fluido se supone incompresible y con propiedades constantes.

Figura 7. Convección forzada en una cavidad cuadrada.

Las Figuras 8 y 9 muestran los valores de la velocidad adimensional (u/u_p) en la línea de simetría vertical B-B mostrada en la Figura 7, en donde se comparan con los resultados obtenidos por Prakash citados en [10, 12].

Figura 8. Velocidad en la línea de simetría vertical (Re = 100).

Figura 9. Velocidad en la línea de simetría vertical (Re = 400).

Entre los puntos que representan los valores de validación y la curva arrojada por el programa desarrollado, hay correspondencia, tanto en forma como en magnitud, observándose desviaciones menores al 8%, aun cuando, los valores de validación corresponden a una malla uniforme de 441 nodos (21×21), mientras que la usada en este trabajo es una malla no estructurada de 245 nodos.

Finalmente, se obtuvo convergencia en la conservación de masa admitiéndose como criterio un desbalance máximo adimensional de 1×10^–7. 

Convección natural en una cavidad cuadrada 

Este es otro problema típico que se muestra de forma esquemática en la Figura 10. Aquí, las paredes verticales son isotermas de temperatura cuyo valor es especificado a Tc y Tf, bajo la condición de que Tc es mayor a Tf. En cuanto a los bordes superior e inferior se suponen adiabáticos.

Figura 10. Convección natural en una cavidad cuadrada.

Respecto al flujo, su comportamiento se modela como incompresible y las propiedades constantes, excepto para el término fuente de la ecuación de cantidad de movimiento en la dirección vertical, donde se utiliza la ecuación de Bussinesq para relacionar las ecuaciones de energía y cantidad de movimiento [8]. Por ultimo, cualquier componente del vector velocidad es cero en las paredes. 

Las Figuras 11 y 12 muestran los resultados de distribución de velocidad en el plano horizontal medio C-C, los cuales se comparan con los obtenidos en una malla de 19×19 por Prakash citados en [10, 12], para dos números de Grashof diferentes.

Figura 11. Distribución de velocidad V en el plano horizontal medio (Gr = 1000).

Figura 12. Distribución de velocidad V en el plano horizontal medio (Gr = 10000).

Desde el punto de vista cualitativo, el comportamiento del campo de flujo sigue el patrón esperado. Hacia la pared caliente, la densidad es baja y, por lo tanto, se motorizan corrientes ascendentes. Por otra parte, hacia la pared fría la densidad es alta y, corrientes descendentes surgen de forma natural. 

Desde el punto de vista cuantitativo, resulta evidente en estas figuras que la solución alcanzada por el programa es prácticamente igual a la de la referencia, con desviaciones menores al 5%. Por lo tanto, se demuestra la correcta implementación de computacional del método. 

Conclusiones 

Los resultados obtenidos al emplear el programa MEFVOC para la solución de problemas de campo escalar difusivo y de flujo de fluidos bidimensional, en régimen laminar, estado estable, con mallas no estructuradas, corresponden con los obtenidos por otros autores al usar programas semejantes. Por lo tanto, el programa de desarrollado ha sido implementado correctamente, y debido a su potencialidad de trabajar con mallas no estructuradas, puede ser aplicado a situaciones de geometría compleja. 

Por otra parte, el programa MEFVOC representa una herramienta útil para la comprensión de los fenómenos de transporte y del método de volumen finito basado en elementos aplicado a mallas no estructuradas, y constituye la primera fase de desarrollo hacia la obtención de un simulador de problemas de transporte de mayor alcance. 

Finalmente, en este trabajo se presentan resultados para casos típicos de prueba y evaluación de técnicas numéricas sobre mallas no estructuradas, usando el método de volumen finito mediante elementos, lo que constituye una información que no es común encontrar en la bibliografía. 

Agradecimientos 

Los autores manifiestan su agradecimiento al Laboratorio de Simulación Computacional, Departamento de Energía, Escuela de Mecánica, Universidad del Zulia, por hacer disponible los recursos computacionales necesarios para la ejecución de este trabajo.

Referencias

1. W.J. Minkowycz, E.M. Sparrow, G. E. Schneider, R. H. Pletcher, “Handbook of Numerical Heat Transfer”, Jhon Wiley & Sons Inc, pag. 422, USA, 1988.        [ Links ]

2. Ferziger Joel, Peric Milovan, “Computacional Methods for fluids dynamics”, Editorial: Springer, Heildenberg, Alemania, 1997.        [ Links ]

3. S. V. Patankar, “Computation of conduction and duct flow heat transfer”, Innovative Reserch, Inc, Minneapolis, USA, 1991.        [ Links ]

4. José Rincón, “Elemento finito basado en el método de volúmenes de control”, trabajo de Ascenso, Universidad del Zulia, Maracaibo, Venezuela, 1990.        [ Links ]

5. Versteeg, Malalasekera, “Computacional Fluid Dynamics”, Longman Group Ltd, Malaysia, 1995.        [ Links ]

6. B.R Baliga, S.V. Patankar: “A control volume finite-element method for two dimensional fluid flow and heat transfer”, Numerical Heat Transfer, vol 6, pp. 245-261, 1983.        [ Links ]

7. Clovis R. Maliska: “Transferencia de calor e macanica dos fluidos computacional”, Livros Técnicos e Científicos Editora, S. A. Rio de Janeiro, Brasil, 1995.        [ Links ]

8. S.V. Patankar, “Numerical Heat Transfer and fluid flow”, McGraw-Hill, pag. 85, Washington, 1980.        [ Links ]

9. Michael J. Row, “A new control volumen based finite element procedure for the numerical solution of the fluid flow and scalar transport equations”, PhD Thesis, Universidad de Waterloo, Ontario, 1985.        [ Links ]

10. Di Mateo, C. “Simulación Computacional de Campos de Flujo usando Volúmenes de Control”. Tesis de Grado, Universidad del Zulia, Facultad de Ingeniería, Venezuela, 1987.        [ Links ]

11. B.R Baliga, S.V. Patankar, T.T. Pham: “Solution of some two dimensional incompressible fluid flow and heat transfer problems, using a control volume finite-element method”, Numerical Heat Transfer, vol 6, pp. 263-282, 1983.        [ Links ]

12. C. Prakash: “An improved control volume finite-element method for heat and mass transfer and for fluid flow using equal order velocity pressure interpolation”, Numerical Heat Transfer, vol 9, pp. 253-271, 1986.        [ Links ]

13. Morillo Trina, “Desarrollo de un programa de elementos finitos de propósitos generales para problemas de convección difusión”, Tesis de grado, Universidad del Zulia, Facultad de Ingeniería, Venezuela, 1997.        [ Links ]

14. M. Reyes, J. Rincon, R. Health, “Simulación de flujo turbulento utilizando el modelo Kappa Epsilon y el método de elementos finitos basados en volúmenes de control”, Memorias del IV CIMENICS, Ciudad Guayana, Venezuela, 1998.        [ Links ]

15. Ivo Wenneker, “Computation of flows using unstructured staggered grids”, Thesis, Delf University of technology, Rotherdam, 2002.        [ Links ]