Método Numérico Basado en Cálculo Exterior Discreto Para la Discretización de la Ecuación de Convección - Difusión

Método Numérico Basado en Cálculo Exterior Discreto Para la Discretización de la Ecuación de Convección - Difusión

Marco A. Noguez^a, Salvador Botello^a, Rafael Herrera^1a

Resumen

Mientras que la discretización con Cálculo Exterior Discreto (DEC) del término elíptico de la ecuación de Transporte, también conocida como la Ecuación de Convección - Difusión, está bien documentada y establecida, la discretización del término convectivo, así como su estabilización siguen siendo un tema de investigación. En este trabajo se propone una discretización local de la Ecuación de Transporte homogénea, isótropa y con flujo incompresible, basada en DEC y fundamentada en argumentos geométricos. Debido a que el método numérico obtenido para el termino convectivo es similar al de Elemento Finito con funciones de interpolación lineales (FEML), se pueden aprovechar algunas técnicas conocidas para la estabilización de la ecuación. Empleando esta característica, las pruebas numéricas se llevan a cabo en mallas que varían de gruesas a finas para comparar los resultados obtenidos con DEC y FEML, y para mostrar convergencia numérica en problemas estacionarios. La estabilización de la ecuación es llevada a cabo mediante Difusión Artificial.

(¹) Apoyado por proyecto CONACYT CB-2015/256126

1 Introducción

Cálculo Exterior Discreto es un método numérico relativamente nuevo para la integración numérica de ecuaciones diferenciales parciales, basado en la discretización de la teoría del Cálculo Diferencial Exterior [7,4,3,5]. Fue propuesto originalmente por A. Hirani en su tesis doctoral en CALTECH [7] y, desde entonces, ha sido empleado exitosamente para la resolución de las ecuaciones de Navier-Stokes y Poisson [9], la ecuación de Darcy [8] y la ecuación de transporte para problemas con convección dominante con flujo incompresible [6]. En este último trabajo, los autores demostraron que DEC coincide con algunos esquemas numéricos como Elemento Finito, Volumen Finito y Diferencias Finitas al considerar problemas simples, lo cual conduce a una método estable tipo upwind, involucrando una discretización para la derivada de Lie.

La ecuación de transporte, conocida también como la ecuación de Convección - Difusión, es una ecuación diferencial parcial que modela el transporte de un campo escalar (e.g. temperatura, masa o momento) debido a procesos de transporte conocidos como convección (para transporte de calor) o advección (para transporte de cualquier sustancia), la cual ocurre siempre que un fluido se encuentra en movimiento, y difusión, debida al movimiento Browniano, principalmente. Esta ecuación surge en diversos problemas en ingeniería como en contaminación de acuíferos, donde el transporte del soluto se debe a la advección producida por la velocidad de partícula, la cual depende, usualmente, del espacio y del tiempo. Estos problemas suelen ser advectivos dominantes, y por lo tanto, producen inestabilidades en la mayoría de los métodos numéricos basados en discretización de dominios. En estos casos, se requiere introducir técnicas de estabilización al modelo numérico tales como "Stream-Upwind/Petrov-Galerkin" (SUPG) [1] y Galerkin/Least-Squares Method (GLS) [10]. Sin embargo, algunas técnicas de estabilización puede que no sean adecuadas para algunos problemas, mientras que refinar el mallado del dominio incrementa los costos computacionales.

En este trabajo, se propone una discretización local basada en DEC de la ecuación de Convección-Difusión para flujo incompresible, a partir de argumentos geométricos. Debido a que la formulación es similar a la de Elemento Finito, la estabilización del termino convectivo se lleva a cabo mediante Difusión Artificial.

El resto de este trabajo se organiza como sigue: En la sección 2 se revisan algunos conceptos básicos de DEC que son esenciales para el desarrollo de este metodo. En la sección 3 se describe el procedimiento para la discretización de la Ecuación de Convección - Difusión. Las pruebas numéricas son llevadas a cabo en mallas que varían de gruesas a finas para comparar DEC y Elemento Finito con funciones de interpolación lineales (FEML), así como para mostrar convergencia numérica. Estas pruebas se describen en la sección 4.

2 Preliminares de Cálculo Exterior Discreto

La integración numérica de ecuaciones diferenciales parciales mediante Cálculo Exterior Discreto consiste en los siguientes pasos:

1. Expresar los operadores diferenciales de la ecuación en lenguaje de Cálculo Diferencial Exterior.
2. Discretizar los operadores involucrados por sus homólogos continuos.

Los operadores diferenciales involucrados en la ecuación diferencial parcial son sustituidos por las siguientes identidades para el gradiente, rotacional, divergencia y laplaciano (ver [2] y [7] para la demostración de estas identidades y familiarización con formas diferenciales)

$\nabla \phi =d\phi ,$ (1.a) $\nabla \times {\underline {u}}u=(\star d{\underline {u}}u^{\flat })^{\sharp },$ (1.b) $\nabla \cdot {\underline {u}}u=\star d\star {\underline {u}}u^{\flat },$ (1.c) $\nabla ^{2}\phi =\star d\star d\phi ,$ (1.d)

donde ${\textstyle \phi }$ es una función escalar (o una ${\textstyle 0-}$forma), ${\textstyle {\underline {u}}u}$ es un campo vectorial, ${\textstyle d}$ es el operador Derivada Exterior y ${\textstyle \star }$ es el operador Estrella de Hodge. Los operadores ${\textstyle (\flat )}$ y ${\textstyle (\sharp )}$ dependen de la métrica inducida por el espacio y transforman vectores a ${\textstyle 1-}$formas diferenciales y viceversa, respectivamente.

Para discretizar los operadores ${\textstyle d}$ y ${\textstyle \star }$ en 2D, consideremos una malla formada por un solo elemento ${\textstyle [v_{1},v_{2},v_{3}]}$ orientado en sentido contrario a las manecillas del reloj, como el mostrado en la Figura 1, formado por los vértices ${\textstyle [v_{1}],[v_{2}],[v_{3}]}$ y aristas ${\textstyle [v_{1},v_{2}]}$, ${\textstyle [v_{2},v_{3}]}$, ${\textstyle [v_{3},v_{1}]}$. El operador frontera ${\textstyle \partial }$ que actúa sobre triángulos orientados (${\textstyle [v_{1},v_{2},v_{3}]}$) y aristas (${\textstyle [v_{1},v_{2}],[v_{2},v_{3}]}$ y ${\textstyle [v_{3},v_{1}]}$), como se menciona en [5], puede ser expresado como un operador matricial, al considerar cada triangulo, arista y vértice como un elemento de una base vectorial. Por lo tanto, para la Figura (1), el operador frontera que actúa sobre ${\textstyle [v_{1},v_{2},v_{3}]}$ y produce una suma orientada de aristas es

$\partial _{2,1}={\begin{pmatrix}1\\1\\-1\end{pmatrix}},$

donde los subíndices ${\textstyle \partial _{2,1}}$ indican que estamos enviando la frontera de un elemento ${\textstyle 2-}$dimensional y obtenemos una suma orientada de elementos ${\textstyle 1-}$dimensionales. De manera similar, el operador frontera que actúa sobre aristas y produce una suma orientada de vértices es

$\partial _{1,0}={\begin{pmatrix}-1&0&1\\1&-1&0\\0&1&-1\end{pmatrix}}.$

Por la dualidad (entre los operadores diferenciación y frontera) presente en el Teorema de Green (vea [5]), la discretización de la derivada exterior ${\textstyle d}$ está dada por:

$d=(\partial )^{T}.$

Por ejemplo, la discretización del operador ${\textstyle d}$ que actúa sobre ${\textstyle 0-}$formas diferenciales y produce ${\textstyle 1-}$formas es

$d_{0,1}={\begin{pmatrix}-1&1&0\\0&-1&1\\0&1&-1\end{pmatrix}}.$

(2)

Elemento triangular [v₁,v₂,v₃] orientado según indican las flechas.

Figura 1: Elemento triangular $[v_{1},v_{2},v_{3}]$ orientado según indican las flechas.

Para discretizar el operador Estrella de Hodge ${\textstyle \star }$, es necesario definir una malla dual que capture la idea de direcciones ortogonales (Ver Figura 3):

• El dual del triangulo ${\textstyle 2-}$dimensional es el circuncentro ${\textstyle 0-}$dimensional del triangulo (punto ${\textstyle C}$ en la Figura 2a).
• El dual de una arista ${\textstyle 1-}$dimensional es el segmento ${\textstyle 1-}$dimensional que une el punto medio de la arista y el circuncentro C (segmento ${\textstyle l_{1}}$, de la Figura 2).
• El dual de un vértice ${\textstyle 0-}$dimensional es el cuadrilátero ${\textstyle 2-}$dimensional formado por los puntos medios de las aristas adyacentes, el circuncentro y el vértice (e.g. el dual del vértice ${\textstyle [v_{1}]}$ en la Figura 3a es el cuadrilátero ${\textstyle [v_{1},p_{1},C,p_{3}]}$).

Por lo tanto, se requieren dos matrices de estrella de Hodge: una que relaciona aristas primales y duales, y otra que relaciona vertices y celdas duales. La primera de ellas está dada por (ver Figura 3):

$M_{1,1}={\begin{pmatrix}{\frac {l_{1}}{L_{1}}}&0&0\\0&{\frac {l_{2}}{L_{2}}}&0\\0&0&{\frac {l_{3}}{L_{3}}}\\\end{pmatrix}},$

(3)

donde los subíndices ${\textstyle M_{1,1}}$ indican que estamos enviando elementos ${\textstyle 1-}$dimensionales de la malla primal a la dual. De manera similar, la segunda matriz Estrella de Hodge está dada por (ver Figura 3):

$M_{0,2}={\begin{pmatrix}A_{1}&0&0\\0&A_{2}&0\\0&0&A_{3}\\\end{pmatrix}},$

(4)

donde ${\textstyle A_{1}}$, ${\textstyle A_{2}}$ y ${\textstyle A_{3}}$ son las áreas de las celdas duales. La matriz inversa de ${\textstyle M_{0,2}}$, la matriz ${\textstyle M_{2,0}}$ envía elementos ${\textstyle 2-}$dimensionales de la malla dual a elementos ${\textstyle 0-}$dimensionales de la malla original.

File:Draft Herrera 100377898-dualVertex circumcenter.png	File:Draft Herrera 100377898-dualEdge.png
(a)	(b)
Figura 2

File:Draft Herrera 100377898-dualCell.png

(a)

Figura 3: Celdas duales

3 Discretización de la Ecuación de Convección - Difusión

La ecuación de Convección - Difusión en 2D homogénea, isótropa y estacionaria con flujo incompresible está dada por:

${\underline {u}}u\cdot (\nabla \phi )=k\nabla ^{2}\phi +q,$

(5)

donde ${\textstyle {\underline {u}}u}$ es la velocidad de partícula, ${\textstyle k}$ es el coeficiente de difusión y ${\textstyle q}$ el término fuente. Se trata de una ecuación elíptica que modela el transporte de una cantidad escalar debido a los procesos de convección, modelado por el lado izquierdo de la Ecuación (5), y por difusión, modelado, a su vez, por el lado derecho de la Ecuación (5).

Para expresar esta ecuación en notación de Cálculo Diferencial Exterior, se sustituyen las Ecuaciones (1.a) y (1.d) en (5)

${\underline {u}}u\cdot (d\phi )=k\star d\star d\phi +q,$

(6)

y, sustituyendo los operadores ${\textstyle \star }$ y ${\textstyle d}$ por sus homólogos discretos (2), (3) y (4) vistos anteriormente, se obtiene

${\underline {u}}u\cdot (d\phi )=kM_{2,0}D_{1,2}^{dual}M_{1,1}D_{0,1}\phi +q,$

(7)

donde ${\textstyle D_{1,2}^{dual}=-D_{0,1}^{T}}$.

Para discretizar el termino ${\textstyle {\underline {u}}u\cdot (d\phi )}$ de la Ecuación (6), consideremos una función ${\textstyle \phi }$ por sus valores discretos en los vértices de un triángulo ${\textstyle [v_{1},v_{2},v_{3}]}$:

$\phi \simeq {\begin{pmatrix}\phi _{1}\\\phi _{2}\\\phi _{3}\end{pmatrix}}$

Podemos aproximar la derivada direccional de ${\textstyle \phi }$ en el vértice ${\textstyle [v_{1}]}$ en dirección del vector ${\textstyle [v_{2}]-[v_{1}]}$ como:

$d(\phi _{[v_{1}]})_{[v_{1},v_{2}]}=\phi _{2}-\phi _{1},$

de manera similar, la derivada direccional de ${\textstyle \phi }$ en dirección del vector ${\textstyle [v_{3}]-[v_{1}]}$ está dada por:

$d(\phi _{[v_{1}]})_{[v_{1},v_{3}]}=\phi _{3}-\phi _{1}.$

Por lo tanto, para encontrar el vector gradiente ${\textstyle {\underline {W}}W_{1}}$ de ${\textstyle \phi }$ en el vértice ${\textstyle [v_{1}]}$, podemos plantear el siguiente sistema de ecuaciones:

${\underline {W}}W_{1}\cdot ([v_{2}]-[v_{1}])=\phi _{2}-\phi _{1},$ ${\underline {W}}W_{1}\cdot ([v_{3}]-[v_{1}])=\phi _{3}-\phi _{1},$

donde:

$[v_{1}]=(x_{1},y_{1}),$ $[v_{2}]=(x_{2},y_{2}),$ $[v_{3}]=(x_{3},y_{3}).$

Resolviendo el sistema anterior para ${\textstyle {\underline {W}}W_{1}}$ se obtiene:

${\underline {W}}W_{1}={\frac {1}{2A}}{\begin{bmatrix}\phi _{1}(y_{2}-y_{3})+\phi _{2}(y_{3}-y_{1})+\phi _{3}(y_{1}-y_{2})\\-(\phi _{1}(x_{2}-x_{3})+\phi _{2}(x_{3}-x_{1})+\phi _{3}(x_{1}-x_{2}))\end{bmatrix}},$

(8)

donde ${\textstyle A}$ es el área del triangulo. Siguiendo el mismo procedimiento para encontrar ${\textstyle {\underline {W}}W_{2}}$ y ${\textstyle {\underline {W}}W_{3}}$, se observa que:

${\underline {W}}W_{1}={\underline {W}}W_{2}={\underline {W}}W_{3}.$

Más aún, el vector gradiente obtenido en la Ecuación (8) coincide con el que se obtiene con FEML. Consideremos ahora la velocidad de partícula ${\textstyle {\underline {v}}v}$ discretizada sobre los vértices del triángulo:

${\underline {u}}u_{1}=(u_{1,1},u_{1,2}),$ ${\underline {u}}u_{2}=(u_{2,1},u_{2,2}),$ ${\underline {u}}u_{3}=(u_{3,1},u_{3,2}).$

Tomando el producto interno del vector gradiente con la velocidad de partícula en cada vértice, obtenemos el operador matricial local para el término convectivo de la Ecuación (5):

${\underline {u}}u\cdot \nabla \phi \simeq V\phi ={\frac {1}{2A}}{\begin{bmatrix}V_{1,1}&V_{1,2}&V_{1,3}\\V_{2,1}&V_{2,2}&V_{2,3}\\V_{3,1}&V_{3,2}&V_{3,3}\end{bmatrix}}{\begin{bmatrix}\phi _{1}\\\phi _{2}\\\phi _{3}\end{bmatrix}},$

(9)

donde:

$V_{1,1}=u_{1,1}(y_{2}-y_{3})+u_{1,2}(x_{3}-x_{2}),$ $V_{2,1}=u_{2,1}(y_{2}-y_{3})+u_{2,2}(x_{3}-x_{2}),$ $V_{3,1}=u_{3,1}(y_{2}-y_{3})+u_{3,2}(x_{3}-x_{2}),$ $V_{1,2}=u_{1,1}(y_{3}-y_{1})+u_{1,2}(x_{1}-x_{3}),$ $V_{2,2}=u_{2,1}(y_{3}-y_{1})+u_{2,2}(x_{1}-x_{3}),$ $V_{3,2}=u_{3,1}(y_{3}-y_{1})+u_{3,2}(x_{1}-x_{3}),$ $V_{1,3}=u_{1,1}(y_{1}-y_{2})+u_{1,2}(x_{2}-x_{1}),$ $V_{2,3}=u_{2,1}(y_{1}-y_{2})+u_{2,2}(x_{2}-x_{1}),$ $V_{3,3}=u_{3,1}(y_{1}-y_{2})+u_{3,2}(x_{2}-x_{1}).$

Finalmente, sustituyendo la Ecuación (9) en (7), obtenemos la discretización basada en DEC de la Ecuación de Transporte para flujo incompresible:

$[M_{0,2}V+kD_{0,1}^{T}M_{1,1}D_{0,1}]\phi =M_{0,2}q.$

4 Pruebas Numéricas

Para las pruebas numéricas se consideró como dominio una elipse cuya ecuación está dada por:

$x^{2}+4y^{2}-100=0,$

con las siguientes características (ver Figura 4):

• Condiciones tipo Dirichlet ${\textstyle \phi =0}$ en la frontera de la elipse.
• Fuentes ${\textstyle q=1}$ sobre dos discos de radio ${\textstyle 0.5}$ y centro en ${\textstyle (-7.5;0)}$ y ${\textstyle (7.5;0)}$, respectivamente.
• Coeficiente de difusión ${\textstyle k=0.5}$.
• Velocidad de partícula ${\textstyle {\underline {v}}v=(2y;-x)}$.

Elipse considerado para las pruebas numéricas

Figura 4: Elipse considerado para las pruebas numéricas

La discretización del dominio se llevó a cabo mediante mallas simpliciales que varían de gruesas a finas y se muestran en la Figura 6.

File:Draft Herrera 100377898-mesh1C.png	File:Draft Herrera 100377898-mesh2C.png
(a)	(b)
Figura 5

File:Draft Herrera 100377898-mesh3C.png	File:Draft Herrera 100377898-mesh4C.png
(a)	(b)
Figura 6: Mallas empleadas para la simulación numérica.

El resultado obtenido con nuestra propuesta basada en DEC (DEC-b), empleando la malla más fina (Figura 6b) se muestra en la Figura 7. Para fines de comparación, el resultado obtenido con FEML empleando la malla más fina (Figura 6b) es tomado como la mejor aproximación de la solución analítica. El Cuadro 1 muestra las características de cada malla, una comparación entre los valores máximos de la función obtenidos con DEC-b y FEML y la discrepancia (DSP) entre los valores de DEC-b y la mejor aproximación de la solución analítica.

Resultado obtenido con DEC-based considerando la malla más fina (19679 nodos y 38388 elementos.)

Figura 7: Resultado obtenido con DEC-based considerando la malla más fina ($19679$ nodos y $38388$ elementos.)

Tabla. 1 Resultados obtenidos en las pruebas numéricas. Se muestra una comparación entre los valores máximos de la función obtenidos con DEC-b y FEML.

			Valores Máximos
Malla	Nodos	Elementos	FEML	DEC-b	DSP
Figura 5a	$491$	$860$	$0.10852$	$0.10714$	$0.01571$
Figura 5	$983$	$1772$	$0.11105$	$0.10983$	$0.01302$
Figura 6a	$5246$	$10010$	$0.13583$	$0.13653$	$0.01368$
Figura 6b	$19679$	$38388$	$0.12285$	$0.12295$	$1\times 10^{-4}$

La Figura 9 muestra gráficas de los valores de la función obtenida con DEC-b y FEML sobre el eje mayor de la elipse, así como una comparación con el resultado obtenido con FEML empleando la malla mostrada en la Figura 6b.

El problema presenta inestabilidades para ambos métodos cuando se emplean las mallas mostradas en las Figuras 5a y 5 y, como consecuencia, la estabilización con difusión artificial es necesaria para evitar oscilaciones en la solución; debido a esto, la solución obtenida con FEML y DEC-b está por debajo de la mejor aproximación de la solución analítica. Como se puede observar en los resultados mostrados, DEC-b y FEML se comportan de manera similar con mallas con pocos elementos. Para mallas finas, DEC-b y FEML convergen al mismo resultado.

File:Draft Herrera 100377898-mesh1C.png	File:Draft Herrera 100377898-mesh2C.png
(a)	(b)
Figura 8

File:Draft Herrera 100377898-mesh3C.png	File:Draft Herrera 100377898-mesh4C.png
(a)	(b)
Figura 9: Valores de temperatura a lo largo de una sección sobre el eje mayor de la elipse para diferentes mallas. $(a)$ Gráfica para la malla mostrada en la Figura 5a; $(b)$ Gráfica para la malla mostrada en la Figura 5; $(c)$ Gráfica para la malla mostrada en la Figura 6a; $(d)$ Gráfica para la malla mas fina.

5 Conclusiones

En este trabajo se propuso una discretización local para la Ecuación de Transporte basada en Cálculo Exterior Discreto. Dada la similitud de esta discretización con Elemento Finito con funciones de interpolación lineales, la estabilización del problema puede llevarse a cabo empleando técnicas conocidas para Elemento Finito, como difusión artificial.

Las pruebas numéricas muestran un comportamiento similar de ambos métodos numéricos para mallas gruesas; para mallas finas, DEC-b y FEML convergen al mismo resultado.

BIBLIOGRAFÍA

[1] Brooks Alexander y Hughes Thomas. Streamline Upwind/Petrov-Galerkin Formulations for Convection Dominated Flows with Particular Emphasis on the Incompressible Navier-Stokes Equations. En: Computer Methods in Applied Mechanics and Engineering (sep. de 1982), p ágs. 199-259.

[2] Keenan Crane y col. Digital Geometry Processing with Discrete Exterior Calculus. En: SIG- GRAPH ’13 (2013).

[3] Keenan Crane. Discrete Differential Geometry (2019) https://www.cs.cmu.edu/ kmcrane/Projects/DDG/paper.pdf

[4] Andrew Gillette. Notes on Discrete Exterior Calculus (2009) https://www.math.arizona.edu/ agillette/research/decNotes.pdf

[5] Humberto Esqueda, Rafael Herrera y Salvador Botello. A Geometric Description of Discrete Exterior Calculus for General Triangulations. En: Revista Internacional de M'etodos Num'ericos para Cálculo y Diseño en Ingeniería 35.1 (feb. de 2019). https://www.scipedia.com/public/Herrera_et_al_2018b

[6] Michael Griebel, Christian Rieger y Alexander Schier. Upwind Schemes for Scalar Advection- Dominated Problems in the Discrete Exterior Calculus. En: Bothe D., Reusken A. (eds) Trans- port Processes at Fluidic Interfaces. Advances in Mathematical Fluid Mechanics. Birkhäuser, Cham (jul. de 2017), págs. 145-175. https://link.springer.com/chapter/10.10072F978-3-319-56602-3_6.

[7] Anil Nirmal Hirani. Discrete Exterior Calculus. Tesis doct. California Institute of Technology, 2003.

[8] Anil Hirani, K Nakshatrala y Jehanzeb Chaudhry. Numerical Method for Darcy Flow Deri- ved Using Discrete Exterior Calculus. En: International Journal for Computational Methods in

Engineering Science and Mech 16 (oct. de 2008). doi: 10.1080/15502287.2014.977500.

[9] Mamdouh Mohamed, Anil Hirani y Ravi Samtaney. Discrete Exterior Calculus of Incompressible Navier - Stokes Equations Over Surface Simplicial Meshes. En: Journal of Computational Physics 312 (mayo de 2016), págs. 175-191. https://doi.org/10.1016/j.jcp.2016.02.028.

[10] Hughes Thomas, Franca Leopoldo y Hulbert G. A New Finite Element Formulation for Computational Fluid Dynamics: VIII. The Galerkin/Least-Squares Method for Advective-Diffusive equations. En: Computer Methods in Applied Mechanics and Engineering (1989), págs. 173-189.