RESUMEN

Este artículo describe una estrategia de discretización de la ecuación de difusión en mallas no estructuradas aplicando el método de los volúmenes finitos con las variables calculadas en el centroide de cada volumen. Esta aproximación está basada en el trabajo de Date [1] que usa una técnica iterativa conocida como corrección diferida para solucionar el cálculo del flujo difusivo en mallas no ortogonales. Se comprobó que para ángulos internos del elemento menores a 50°, el método propuesto por Date no converge y entonces se propone una nueva forma de calcular el gradiente que favorece la convergencia del problema. Se muestra un estudio de la convergencia donde se demuestra la alta efectividad del método propuesto. A partir de la solución de un problema típico, basado en la solución de la ecuación de Poisson, se realizó la comparación de los resultados obtenidos con los de la solución analítica, donde se observó una alta correspondencia de los resultados sin comprometer el tiempo de cálculo. Finalmente, se demostró la flexibilidad de la aproximación implementada al realizar simulaciones sobre mallas estructuradas y no estructuradas, usando elementos en forma de cuadriláteros y triángulos en 2D, y cubos curvilíneos y tetraedros en 3D.

Palabras clave:CFD, Volúmenes Finitos, Mallas no estructuradas, flujo difusivo, dominio complejo.

ABSTRACT

This paper describes a strategy to discretize the Poisson equation on unstructured meshes using the finite volume method based on the center of the cells. This approach is based on Date’s work[1] that uses an iterative technique known as deferred correction to obtain the right diffusive flow field in no orthogonal meshes. It was found that the method proposed by Date does not converge when the internal angles are less than 40°, then we proposes a new way to calculate the gradient in order to ensure the convergence. It shows a convergence study that demonstrates the high effectiveness of the proposed method. After solving a typical problem, based on the solution of the Poisson equation, we compared the results obtained with the analytical solution, where there was a high correspondence of results without compromising the computational time. Finally, we have demonstrate the flexibility of the approach implemented by performing simulations on structured and unstructured meshes, using elements in the form of quadrilaterals and triangles in 2D, and curvilinear cubes and tetrahedrons in 3D.

Keywords:CFD, Finite Volume Method, Unstructured Meshes, diffusive flux, complex domain

UIS Ingenierías, Volumen 11, No. 1, pags. 45 - 54, Junio 2012; Facultad de Ingenierías Fisicomecánicas, UIS

1. INTRODUCCIÓN

En prácticamente todas las aplicaciones reales de CFD se trabaja con dominios complejos. Un dominio se denomina complejo cuando no se puede describir o representar convenientemente en una malla cartesiana. Un dominio complejo, por lo general, está construido de elementos o volúmenes de control en forma de triángulos (en 2D) o tetraedros (en 3D) que pueden provenir de un mallado no estructurado; o en cuadriláteros curvilíneos (en 2D) o cubos curvilíneos (en 3D) si el mallado es estructurado. Estos tipos de elementos se usan debido a que permiten representar de forma más detallada el dominio y en particular la frontera. La desventaja que representa este tipo de discretización es que la estructura que se genera conduce a un consumo de tiempo de cálculo mucho mayor que para el caso de dominios cartesianos. El método general de cálculo basado en la aplicación de volúmenes de control comienza con la solución de una ecuación de conservación genérica del tipo:

Cuando se trabaja con geometrías complejas, los términos de la izquierda tienen tratamiento relativamente sencillo comparado con los pasos involucrados en la evaluación del componente difusivo. Si se define el flujo difusivo a través de una superficie S para un volumen de control como:

Usando el teorema de Gauss en la Ec (2), el flujo difusivo se puede evaluar como:

donde S^p_k es el área de flujo a través de la frontera del volumen P (∇ Ø ) ^k_p, es el gradiente de la variable Ø en la frontera k, es el vector normal a la superficie k (dirigido hacia afuera). El principal problema de la discretización del flujo difusivo en los dominios complejos se reduce al cálculo de una forma consistente del gradiente de Ø en cada una de las fronteras del volumen. El cálculo del gradiente está directamente relacionado con la estructura de la malla generada. Si el mallado es estructurado, el vector de área normal será paralelo a la línea que une los centroides del volumen P y cada uno de sus vecinos. Adicionalmente, la línea que une los dos centroides pasará por el punto medio de la respectiva frontera del volumen, punto de la Figura 1. Estas dos características hacen, como se verá más adelante, que la evaluación del gradiente sea trivial. Si el dominio es complejo, estas dos características no se cumplen y hacen que el cálculo del gradiente sea también complejo. Una forma de facilitar la evaluación del gradiente en dominios complejos es realizar el mallado de tal manera que se produzca una malla tal que la línea de unión de los dos centroides sea paralela al vector normal al área [2]. La desventaja de este método está en que no siempre es posible construir este tipo de mallas, y primordialmente, la mayoría del software para mallado no posee esta característica. Por otro lado, muchos otros autores han propuesto aproximaciones para el cálculo del término en mallas generales [3], [4], [5], [6], [7], [8].

La idea principal de estas aproximaciones es el uso de variables adicionales, que incrementan la molécula computacional y de esta forma la complejidad de la matriz. La solución de este sistema de ecuaciones es mucho más costosa en términos de tiempo de cálculo. Esta dificultad se puede evitar usando una estrategia muy conocida en CFD propuesta por Muzaferija [9], divulgada por Peric [4], que ha sido actualizada por diversos autores tales como Traoré [7] y Date [3] conocida como corrección diferida. Bajo este contexto, la corrección diferida es un método iterativo para aproximar el flujo difusivo en la solución de las ecuaciones de Navier-Stokes, en problemas basados en la ecuación de Poisson o de Laplace en mallas no ortogonales. La aproximación del término se puede realizar a partir del proceso iterativo referenciado en [4], [7] y [1], entre otros, el cual debe converger hacia el valor correcto de cuando el número de iteraciones tiende a infinito. En el trabajo de Ferziger [4] se plantea la corrección diferida en una dirección perpendicular a la línea que une los centroides de los elementos vecinos, dirección en la Figura 1. De acuerdo con Traoré [7], esta estrategia presenta problemas de divergencia para ángulos muy obtusos, especialmente si los elementos son en forma de cuadriláteros. Desde esta perspectiva, Traoré [7] sugiere el uso de una corrección diferida mejorada donde la dirección de corrección se toma ahora en una dirección tangente a la frontera de análisis, dirección en la Figura 2. Sin embargo, como los autores comentan, este método presenta mayores errores de cálculo si existe una gran distorsión general en la malla. Este aumento del error puede estar relacionado con el hecho que en la formulación de Traoré no se evalúa el gradiente en el centro de la frontera. Date [3] presentó una metodología para calcular el gradiente usando una corrección diferida tomando la dirección de corrección paralela a la frontera, pero asegurando la evaluación del gradiente en el centro del área. Sin embargo, este último trabajo no garantiza la convergencia para ángulos muy cerrados. En este artículo se presenta una metodología para el cálculo del flujo difusivo, basado en el trabajo de Date [1], realizando un modificación en el cálculo del gradiente, que evita los problemas de convergencia, incluso para ángulos muy agudos.

2. METODOLOGÍA DE DISCRETIZACIÓN DEL FLUJO DIFUSIVO

El flujo difusivo se encuentra presente en la solución de las ecuaciones de Navier-Stokes, en problemas relacionados con la ecuación de Poisson o la ecuación de Laplace. Este trabajo se soluciona la ecuación de Poisson, que en un dominio n-dimensional se puede escribir como:

con Ø = Ø_o en S, es decir, con unos valores o condiciones conocidos en la frontera. f y Γ son funciones definidas en el dominio Ω, es la variable desconocida, y S es la frontera de Ω. Para resolver la Ec. (4) por volúmenes finitos se debe realizar una integración sobre un volumen diferencial, obteniéndose:

Aplicando el teorema de Green en el lado izquierdo de la Ec. (5) queda:

donde el lado izquierdo de la Ec. (6) se puede interpretar como el balance del flujo de la variable Ø a través de las fronteras del elemento, y el lado derecho como la velocidad de generación o destrucción de la variable en el elemento.

2.1 Discretización del flujo difusivo

En la Figura 3 se muestra un elemento típico para la discretización por volúmenes finitos. Cabe anotar que, para que el método implementado funcione, no es necesario que el elemento provenga de un mallado estructurado o no estructurado. La metodología presentada puede funcionar para cualquier geometría implementando el cálculo de ciertos valores geométricos que se enumerarán más adelante. La formulación se realiza considerando el problema en dos dimensiones, pero se deja al lector demostrar que se puede extender a los casos de una y tres dimensiones.

En la Figura 3 se muestra un volumen principal rodeado de un volumen vecino etiquetado con la letra E. Para la formulación basada en volúmenes finitos se desea realizar un balance de la variable Ø en el punto medio de las fronteras del volumen (punto c en la figura). Denotando D_k como el flujo difusivo de la variable Ø en las fronteras y como el flujo que proviene de los volúmenes vecinos y llega al volumen a través de la frontera k, se tiene entonces:

que como se mencionó anteriormente se debe evaluar en el punto medio de cada frontera. La evaluación del flujo difusivo de la variable Ø entre el volumen P y un vecino E en el centro de la frontera (punto de la Figura 3) se puede evaluar mediante:

Observe que la difusión se evalúa directamente en términos de un gradiente normal a la superficie. Esta característica permite que la evaluación del transporte difusivo planteado en este artículo se pueda realizar en poliedros en 3D.

Para la evaluación del transporte difusivo en la Ec. (8), el área normal a la frontera se conocerá a partir de la información geométrica del elemento.Γ_c representa la constante difusiva en la frontera, que se puede evaluar a partir del peso geométrico de cada volumen a cada lado en la frontera. Lo que se desea evaluar a continuación es el gradiente en dirección normal de la variable Ø en el punto c. Para realizar esta operación, se debe reconocer que el punto c, en general, no estará en el punto medio entre los puntos P₂ y E₂ de la Figura 3. Por lo tanto, para conservar un segundo orden de exactitud en la evaluación del gradiente, se emplea una serie de Taylor expandida a partir de los punto P₂ y E₂.

Eliminando la segunda derivada a partir de las dos ecuaciones anteriores, se puede demostrar según Date [1] que:

Combinando las Ec. (11) y (12), se puede expresar el transporte difusivo total como

La ecuación anterior no tiene ningún uso práctico debido a que usa valores de la variable en puntos ficticios como P₂, E₂ y en los vértices a y b que son desconocidos. Para poder usar esta formulación se debe evaluar los valores de a partir de Ø los valores en los puntos conocidos P y E.

Interpolación de Ø en P₂, E₂, a y b

Si se asume que la variación de Ø entre P y P₂ es lineal, se puede usar una expansión lineal,

Donde

y

Bajo argumentos similares, se pueden calcular las restantes variables desconocidas

y

onde NE es el número de elementos vecinos al nodo.

2.2 Ecuación discretizada final

Substituyendo las Ec. (15) a (20) en la Ec. (13) y realizando algunas operaciones algebraicas, se puede obtener la ecuación discretizada

donde

A partir de la deducción de la Ec. (22) se pueden realizar los siguientes comentarios:

2.3 Evaluación del gradiente de Ø

Para evaluar el término D_i en la Ec. (22) se debe evaluar el gradiente de Ø en coordenadas cartesianas, esta evaluación se puede realizar mediante la aplicación del teorema de Gauss, donde se aproxima la derivada en el centro del volumen de control a través del valor medio en la celda:

Se puede considerar entonces la derivada como la divergencia del vector Ø î y transformar la integral sobre el volumen en una integral de superficie[4]:

La ecuación anterior muestra que para calcular el gradiente de con respecto a x en el centro del volumen de control, se debe realizar la sumatoria de los productos Ø de con la componente x del vector área en todas las fronteras del volumen de control. Si se divide por el volumen de la celda:

La formulación presentada, se tiene como puntos de referencia P₂ y E₂, para la evaluación del gradiente se deben usar estos puntos como referencia. En este caso, Ø _c,k se calcula mediante:

La aparición de Δ Ø_p en la Ec. (26) sugiere que la Ec. (25) es implícita en la evaluación de . Para resolver este inconveniente se propone la solución mediante un proceso iterativo. Dado que el cálculo del gradiente se hace a partir de un esquema de segundo orden, y como se verá más adelante, para mallas muy deformadas no siempre se garantiza la convergencia de la solución. Una solución para esta problemática es disminuir el grado de exactitud para el cálculo del gradiente. En este caso también se aumenta el error de truncamiento correspondiente en la solución global. Otra solución puede ser relajar el cálculo del gradiente a partir de los valores disponibles del gradiente de la iteración anterior, mediante el uso de un factor de relajación ω_g. De esta forma, 2.4

2.4 Evaluación de las condiciones de frontera

Para describir la aplicación de las condiciones de frontera, considere un volumen de control cercano a la frontera, ver la Figura 4, con el lado ab coincidiendo con la frontera del dominio. Se define un punto B en la frontera tal que:

La evaluación del flujo difusivo de la variable en la frontera se puede determinar mediante:

Para este caso, el gradiente de Ø se representará mediante una fórmula en diferencias retrasada de primer orden

donde, una vez más,

Reemplazando las Ec. (30) y (31) en la Ec. (29), queda:

donde el coeficiente de difusión se puede evaluar mediante:

Se evaluarán tres condiciones de frontera comunes:

Valor de conocido en la frontera

Inicialmente se reescribirá la Ec. (21) tal que:

donde nb - B implica que la frontera que contribuye con su condición de frontera queda excluida de la sumatoria y será evaluada a partir del término -D_B. Tal como en la formulación tradicional de volúmenes finitos, la condición de frontera se incluye a través de los términos S_u y S_p como

De esta forma

Valor del flujo de Ø conocido en la frontera

En el caso en que se especifique un flujo en la frontera se conoce, y apelando nuevamente a la Ec. (35), se tiene que:

Valor del flujo de Ø asociado con un coeficiente convectivo en la frontera.

En el caso en que se especifique un valor asociado con un coeficiente convectivo en la frontera se conoce, y reemplazando en la Ec. (35), se tiene que:

De esta forma,

donde

Las condiciones de frontera aquí presentadas típicamente involucran condiciones físicas tales como flujos de entrada, paredes, líneas de simetría, flujos de salida.

2.5 Evaluación del gradiente de en la frontera

En el caso en que un lado del volumen pertenezca a la frontera del dominio, la evaluación del gradiente no se puede realizar tal como se muestra en la Ec. (25). Para cada condición de frontera especificada anteriormente se debe imponer el valor de Ø _c correspondiente. A continuación se presenta el valor Ø _c de para cada condición de frontera.

1. Valor de Ø conocido en la frontera. En este caso Ø _c = Ø _b.

2. Valor del flujo de Ø _c = Ø _b conocido en la frontera. A partir de la Ec. (9), y utilizando una aproximación lineal de primer orden , donde , siendo G_n el valor del flujo especificado en la frontera. El valor de se calcular a partir de los valores conocidos de la iteración anterior.

3. Valor de Ø_P2 asociado con un coeficiente convectivo en la frontera. En este caso se deja al lector demostrar que .

2.6 Procedimiento de cálculo general

La Ec. (21) se escribe para cada volumen de control definido en el dominio de cálculo, que junto con la aplicación de las condiciones de frontera, produce un sistema de ecuaciones de la forma Ax = b que se resuelven mediante un método conveniente para este tipo de ecuaciones. La solución del sistema de ecuaciones se realiza mediante dos lazos de control diferentes. El primer lazo, que se denominarán iteraciones mayores, están asociadas con la evaluación del término en la Ec. (21) a partir de valores de de la iteración anterior. En segundo lazo, el cual se denominará lazo interior, está asociado con el tiempo que toma el solucionador del sistema de ecuaciones resultante en resolver el sistema. De esta forma se establecen dos criterios de parada, uno para las iteraciones mayores definidos como para todos las celdas en el dominio. El segundo criterio de parada (o convergencia) está asociado con el solucionador de ecuaciones, si se utilizase un solucionador de ecuaciones iterativo.

2.7 Evaluación de los parámetros geométricos involucrados

En la Figura 5 se muestra un elemento típico al que se le deben definir las cantidades geométricas listadas más adelante. Dada la variedad de elementos que se pueden emplear, se especificarán en este artículo aquellas cantidades que son comunes a todos los posibles elementos a usar. Cabe mencionar que el desarrollo mostrado aquí se implementó para elementos cuadriláteros curvilíneos y triangulares en 2D y elementos cúbicos curvilíneos y tetraedros en 3D. Para estos tipos de elementos se proporcionará alguna referencia donde se puede encontrar el cálculo respectivo. En la Tabla 1 se listan los diferentes parámetros geométricos necesarios. En la tabla son las coordenadasde cada uno de los nodos pertenecientes al elemento .

3. VALIDACIÓN NUMÉRICA

En esta sección se muestran dos experimentos numéricos usados para comprobar la conveniencia del uso de la formulación implementada.

La primera prueba consiste en aplicar una ecuación de Poisson conocida sobre un dominio y comprobar los resultados numéricos obtenidos a partir de la simulación con los resultados analíticos. En este caso se utilizará el dominio mostrado en la Figura 6 variando el ángulo de inclinación.En este artículo la ecuación a resolver es tomada de [7]:

La solución exacta está dada por Ø (x,y)= x³ + y² + xc + C En este artículo se ha tomado la constante C=10 para evitar problemas con los errores de comparación de valores cercanos al cero cuando las coordenadas son cercanas al origen. El ángulo de inclinación se varió de acuerdo con los valores: 0, 20, 30, 38, 45, 50, 60, 70, 80, 85, mientras que se dejó el parámetro de relajación ωg inicialmente igual a cero, teniendo en este caso, el método original propuesto por Date en 2005.

Se desea estudiar la sensibilidad del método implementado con respecto al ángulo de inclinación representado por θ. En la Tabla 2 se muestra el error máximo para cada uno de los ángulos. Los valores no presentados de los otros ángulos indican que la solución diverge para esos casos. Se pudo comprobar que para una malla con elementos cuadriláteros no es posible obtener una convergencia para valores de θ superiores a 40°. Se observó, que cuando el sistema converge, el error máximo es independiente del ángulo θ y depende exclusivamente del criterio de parada establecido. En la tabla también se reporta el número de iteraciones mayores requeridas para alcanzar la convergencia. La tolerancia impuesta como criterio de parada es impuesta a ε = 1x10-4. Como se observa, a medida que el ángulo de inclinación aumenta, el número de iteraciones también lo hace. Con este estudio se pudo determinar que cuando el sistema no convergía, no lo hacía debido a que se presentaban grandes oscilaciones en la evaluación del gradiente entre una iteración y otra. Dado que para valores de θ mayores a 40° se encontró que no era posible alcanzar la convergencia del sistema, se realizaron estudios de sensibilidad del parámetro ω_g.

En este estudio, se varió en el intervalo en incrementos de para el rango de ωg propuesto anteriormente, en las Tabla 3 y 4 se muestran los resultados obtenidos. Sólo se presentan los valores hasta ωg <0.4 debido a que se verifican las mismas tendencias para todo el rango de valores. Como se muestra en la tabla, es posible alcanzar la convergencia para todo el rango de ángulos propuestos al variar el factor de relajación. También se observa, que para el criterio de parada impuesto para la convergencia en las iteraciones mayores, en general, el error de la solución numérica con respecto a la solución analítica aumenta a medida que aumenta el factor de relajación. Una tendencia igual se observa en el número de iteraciones. Adicionalmente, se comprobó que si se desea que el error de cálculo tenga igual orden de magnitud en todos los casos, se requerirán mayor número iteraciones. Cabe anotar que aunque en este ejercicio se requirieron valores de relajación bajos, no es seguro que en todos los problemas se conserve esta tendencia. Se comprobó que para algunos problemas en 3D se necesitaron valores de relajación hasta de ω_g=0.8.

Se realizó un segundo estudio numérico donde se deseaba experimentar sobre la adaptabilidad del método a otras formas de volúmenes de control y otras estructuras de malla. Para este estudio, se tomaron dos tipos de malla para el rango de ángulos indicados anteriormente: (a) a partir de un mallado estructurado se generó una malla de elementos triangulares, ver Figura 7; y (b) a partir de un mallado no estructurado se generó una malla de elementos triangulares, ver Figura 8. Cuando la malla se obtiene a partir de un mallado estructurado, la distorsión de la malla puede ser mínima y, por lo tanto, la línea que une los centroides de los elementos pasaría en la gran mayoría de los casos por el centro de cada frontera. Como consecuencia, no podría ser posible verificar la importancia de asegurar el cálculo del gradiente en el punto medio de la frontera, lo cual es evidente desde el punto de vista físico, que debe ser el punto de evaluación del gradiente. En la Tabla 5 se muestran los resultados de las simulaciones para cada una de las mallas. En todos los casos el error respecto de la solución analítica fue menor a 1x10^-4.

Como se observa en la Tabla 5, la ubicación del punto de evaluación del gradiente influye en la convergencia del sistema, indicando que se debe tener en cuenta. Sin embargo, con el estudio presentado, no es evidente la influencia de este aspecto en la exactitud de la solución. Se pudo verificar mediante la ejecución del mismo estudio con el método propuesto por Traoré [7] que no tener en cuenta este aspecto genera un error en la solución de un orden de magnitud mayor respecto del obtenido con el método aquí desarrollado, inclusive para mallas poco distorsionadas.

4. OBSERVACIONES

En este artículo se describe una metodología para la discretización del flujo difusivo en una malla general basado en el trabajo de Date [1]. Este método es aplicable a un elemento de cualquier topología, siempre que se pueda disponer de la información geométrica correspondiente. Se comprobó que para ángulos internos del elemento menores a 50°, el método propuesto por Date no converge y entonces se propuso una nueva forma de calcular el gradiente que permitió obtener la convergencia del problema en cada uno de los casos planteados. Se demostró la eficacia y eficiencia del método en los casos donde el método original diverge. En los ejemplos considerados, se observa que se conserva el segundo orden de exactitud aún para mallas muy deformadas.

NOMENCLATURA