|
|
| Line 1: |
Line 1: |
| − | ==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<sup>a</sup>, Salvador Botello<sup>a</sup>, Rafael Herrera<span id="fnc-1"></span>[[#fn-1|<sup>1</sup>]]<sup>a</sup>'''
| + | == Abstract == |
| | | | |
| − | ==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.
| |
| | | | |
| − | <span id="fn-1"></span>
| + | == Full document == |
| − | <span style="text-align: center; font-size: 75%;">([[#fnc-1|<sup>1</sup>]]) Apoyado por proyecto CONACYT CB-2015/256126</span>
| + | <pdf>Media:Draft_Herrera_100377898-7751-document.pdf</pdf> |
| − | | + | |
| − | ==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 <span id='citeF-7'></span><span id='citeF-4'></span><span id='citeF-3'></span><span id='citeF-5'></span>[[#cite-7|[7,4,3,5]]]. Fue propuesto originalmente por A. Hirani en su tesis doctoral en CALTECH <span id='citeF-7'></span>[[#cite-7|[7]]] y, desde entonces, ha sido empleado exitosamente para la resolución de las ecuaciones de Navier-Stokes y Poisson <span id='citeF-9'></span>[[#cite-9|[9]]], la ecuación de Darcy <span id='citeF-8'></span>[[#cite-8|[8]]] y la ecuación de transporte para problemas con convección dominante con flujo incompresible <span id='citeF-6'></span>[[#cite-6|[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) <span id='citeF-1'></span>[[#cite-1|[1]]] y Galerkin/Least-Squares Method (GLS) <span id='citeF-10'></span>[[#cite-10|[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:
| + | |
| − | | + | |
| − | <ol>
| + | |
| − | | + | |
| − | <li>Expresar los operadores diferenciales de la ecuación en lenguaje de Cálculo Diferencial Exterior. </li>
| + | |
| − | <li>Discretizar los operadores involucrados por sus homólogos continuos. </li>
| + | |
| − | | + | |
| − | </ol>
| + | |
| − | | + | |
| − | Los operadores diferenciales involucrados en la ecuación diferencial parcial son sustituidos por las siguientes identidades para el gradiente, rotacional, divergencia y laplaciano (ver <span id='citeF-2'></span>[[#cite-2|[2]]] y <span id='citeF-7'></span>[[#cite-7|[7]]] para la demostración de estas identidades y familiarización con formas diferenciales)
| + | |
| − | | + | |
| − | <span id="eq-1.a"></span>
| + | |
| − | <span id="eq-1.b"></span>
| + | |
| − | <span id="eq-1.c"></span>
| + | |
| − | <span id="eq-1.d"></span>
| + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>\nabla \phi = d \phi , </math>
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;" | (1.a)
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math> \nabla \times \underline{u}u = (\star d \underline{u}u^{\flat })^{\sharp }, </math>
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;" | (1.b)
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math> \nabla \cdot \underline{u}u = \star d \star \underline{u}u^{\flat }, </math>
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;" | (1.c)
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math> \nabla ^{2}\phi = \star d \star d \phi , </math>
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;" | (1.d)
| + | |
| − | |}
| + | |
| − | |}
| + | |
| − | | + | |
| − | donde <math display="inline">\phi </math> es una función escalar (o una <math display="inline">0-</math>forma), <math display="inline">\underline{u}u</math> es un campo vectorial, <math display="inline">d</math> es el operador Derivada Exterior y <math display="inline">\star </math> es el operador Estrella de Hodge. Los operadores <math display="inline">(\flat )</math> y <math display="inline">(\sharp )</math> dependen de la métrica inducida por el espacio y transforman vectores a <math display="inline">1-</math>formas diferenciales y viceversa, respectivamente.
| + | |
| − | | + | |
| − | Para discretizar los operadores <math display="inline">d</math> y <math display="inline">\star </math> en 2D, consideremos una malla formada por un solo elemento <math display="inline">[v_{1},v_{2},v_{3}]</math> orientado en sentido contrario a las manecillas del reloj, como el mostrado en la Figura [[#img-1|1]], formado por los vértices <math display="inline">[v_{1}], [v_{2}], [v_{3}]</math> y aristas <math display="inline">[v_{1},v_{2}]</math>, <math display="inline">[v_{2},v_{3}]</math>, <math display="inline">[v_{3},v_{1}]</math>. El operador frontera <math display="inline">\partial </math> que actúa sobre triángulos orientados (<math display="inline">[v_{1},v_{2}, v_{3}]</math>) y aristas (<math display="inline">[v_{1},v_{2}], [v_{2},v_{3}]</math> y <math display="inline">[v_{3},v_{1}]</math>), como se menciona en <span id='citeF-5'></span>[[#cite-5|[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 ([[#img-1|1]]), el operador frontera que actúa sobre <math display="inline">[v_{1},v_{2},v_{3}]</math> y produce una suma orientada de aristas es
| + | |
| − | | + | |
| − | <span id="eq-2"></span>
| + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>\partial _{2,1} = \begin{pmatrix}1 \\ 1 \\ -1 \end{pmatrix}, </math>
| + | |
| − | |}
| + | |
| − | |}
| + | |
| − | | + | |
| − | donde los subíndices <math display="inline">\partial _{2,1}</math> indican que estamos enviando la frontera de un elemento <math display="inline">2-</math>dimensional y obtenemos una suma orientada de elementos <math display="inline">1-</math>dimensionales. De manera similar, el operador frontera que actúa sobre aristas y produce una suma orientada de vértices es
| + | |
| − | | + | |
| − | <span id="eq-2"></span>
| + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>\partial _{1,0} = \begin{pmatrix}-1 & 0 & 1 \\ 1 & -1 & 0 \\ 0 & 1 & -1 \end{pmatrix}. </math>
| + | |
| − | |}
| + | |
| − | |}
| + | |
| − | | + | |
| − | Por la dualidad (entre los operadores diferenciación y frontera) presente en el Teorema de Green (vea <span id='citeF-5'></span>[[#cite-5|[5]]]), la discretización de la derivada exterior <math display="inline">d</math> está dada por:
| + | |
| − | | + | |
| − | <span id="eq-2"></span>
| + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>d = (\partial )^{T}. </math>
| + | |
| − | |}
| + | |
| − | |}
| + | |
| − | | + | |
| − | Por ejemplo, la discretización del operador <math display="inline">d</math> que actúa sobre <math display="inline">0-</math>formas diferenciales y produce <math display="inline">1-</math>formas es
| + | |
| − | | + | |
| − | <span id="eq-2"></span>
| + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>d_{0,1} = \begin{pmatrix}-1 & 1 & 0 \\ 0 & -1 & 1 \\ 0 & 1 & -1 \end{pmatrix}. </math>
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;" | (2)
| + | |
| − | |}
| + | |
| − | | + | |
| − | <div id='img-1'></div>
| + | |
| − | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
| + | |
| − | |-
| + | |
| − | |[[Image:Draft_Herrera_100377898-orientedTriangle1.png|300px|Elemento triangular [v₁,v₂,v₃] orientado según indican las flechas.]]
| + | |
| − | |- style="text-align: center; font-size: 75%;"
| + | |
| − | | colspan="1" | '''Figura 1:''' Elemento triangular <math>[v_{1},v_{2},v_{3}]</math> orientado según indican las flechas.
| + | |
| − | |}
| + | |
| − | | + | |
| − | Para discretizar el operador Estrella de Hodge <math display="inline">\star </math>, es necesario definir una malla dual que capture la idea de direcciones ortogonales (Ver Figura [[#img-3|3]]):
| + | |
| − | | + | |
| − | * El dual del triangulo <math display="inline">2-</math>dimensional es el circuncentro <math display="inline">0-</math>dimensional del triangulo (punto <math display="inline">C</math> en la Figura [[#img-2a|2a]]).
| + | |
| − | * El dual de una arista <math display="inline">1-</math>dimensional es el segmento <math display="inline">1-</math>dimensional que une el punto medio de la arista y el circuncentro C (segmento <math display="inline">l_{1}</math>, de la Figura [[#img-2|2]]).
| + | |
| − | * El dual de un vértice <math display="inline">0-</math>dimensional es el cuadrilátero <math display="inline">2-</math>dimensional formado por los puntos medios de las aristas adyacentes, el circuncentro y el vértice (''e.g.'' el dual del vértice <math display="inline">[v_{1}]</math> en la Figura [[#img-3a|3a]] es el cuadrilátero <math display="inline">[v_{1},p_{1},C,p_{3}]</math>).
| + | |
| − | | + | |
| − | 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 [[#img-3|3]]):
| + | |
| − | | + | |
| − | <span id="eq-3"></span>
| + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>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}, </math>
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;" | (3)
| + | |
| − | |}
| + | |
| − | | + | |
| − | donde los subíndices <math display="inline">M_{1,1}</math> indican que estamos enviando elementos <math display="inline">1-</math>dimensionales de la malla primal a la dual. De manera similar, la segunda matriz Estrella de Hodge está dada por (ver Figura [[#img-3|3]]):
| + | |
| − | | + | |
| − | <span id="eq-4"></span>
| + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>M_{0,2} = \begin{pmatrix}A_{1} & 0 & 0 \\ 0 & A_{2} & 0 \\ 0 & 0 & A_{3} \\ \end{pmatrix}, </math>
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;" | (4)
| + | |
| − | |}
| + | |
| − | | + | |
| − | donde <math display="inline">A_{1}</math>, <math display="inline">A_{2}</math> y <math display="inline">A_{3}</math> son las áreas de las celdas duales. La matriz inversa de <math display="inline">M_{0,2}</math>, la matriz <math display="inline">M_{2,0}</math> envía elementos <math display="inline">2-</math>dimensionales de la malla dual a elementos <math display="inline">0-</math>dimensionales de la malla original.
| + | |
| − | | + | |
| − | <div id='img-2a'></div>
| + | |
| − | <div id='img-2b'></div>
| + | |
| − | <div id='img-2'></div>
| + | |
| − | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
| + | |
| − | |-
| + | |
| − | |[[Image:Draft_Herrera_100377898-dualVertex_circumcenter.png|600px|]]
| + | |
| − | |[[Image:Draft_Herrera_100377898-dualEdge.png|600px|]]
| + | |
| − | |- style="text-align: center; font-size: 75%;"
| + | |
| − | | (a)
| + | |
| − | | (b)
| + | |
| − | |- style="text-align: center; font-size: 75%;"
| + | |
| − | | colspan="2" | '''Figura 2'''
| + | |
| − | |}
| + | |
| − | <div id='img-3a'></div>
| + | |
| − | <div id='img-3'></div>
| + | |
| − | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
| + | |
| − | |-
| + | |
| − | |[[Image:Draft_Herrera_100377898-dualCell.png|300px|]]
| + | |
| − | |- style="text-align: center; font-size: 75%;"
| + | |
| − | | (a)
| + | |
| − | |- style="text-align: center; font-size: 75%;"
| + | |
| − | | colspan="1" | '''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:
| + | |
| − | | + | |
| − | <span id="eq-5"></span>
| + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>\underline{u}u\cdot (\nabla \phi ) = k \nabla ^{2}\phi + q, </math>
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;" | (5)
| + | |
| − | |}
| + | |
| − | | + | |
| − | donde <math display="inline">\underline{u}u</math> es la velocidad de partícula, <math display="inline">k</math> es el coeficiente de difusión y <math display="inline">q</math> 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 ([[#eq-5|5]]), y por difusión, modelado, a su vez, por el lado derecho de la Ecuación ([[#eq-5|5]]).
| + | |
| − | | + | |
| − | Para expresar esta ecuación en notación de Cálculo Diferencial Exterior, se sustituyen las Ecuaciones ([[#eq-1.a|1.a]]) y ([[#eq-1.d|1.d]]) en ([[#eq-5|5]])
| + | |
| − | | + | |
| − | <span id="eq-6"></span>
| + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>\underline{u}u \cdot (d\phi ) = k\star d \star d \phi + q, </math>
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;" | (6)
| + | |
| − | |}
| + | |
| − | | + | |
| − | y, sustituyendo los operadores <math display="inline">\star </math> y <math display="inline">d</math> por sus homólogos discretos ([[#eq-2|2]]), ([[#eq-3|3]]) y ([[#eq-4|4]]) vistos anteriormente, se obtiene
| + | |
| − | | + | |
| − | <span id="eq-7"></span>
| + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>\underline{u}u \cdot (d\phi ) = k M_{2,0}D_{1,2}^{dual}M_{1,1}D_{0,1} \phi + q, </math>
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;" | (7)
| + | |
| − | |}
| + | |
| − | | + | |
| − | donde <math display="inline">D_{1,2}^{dual} = -D_{0,1}^{T}</math>.
| + | |
| − | | + | |
| − | Para discretizar el termino <math display="inline">\underline{u}u \cdot (d \phi )</math> de la Ecuación ([[#eq-6|6]]), consideremos una función <math display="inline">\phi </math> por sus valores discretos en los vértices de un triángulo <math display="inline">[v_{1},v_{2},v_{3}]</math>:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>\phi \simeq \begin{pmatrix}\phi _{1} \\ \phi _{2} \\ \phi _{3} \end{pmatrix} </math>
| + | |
| − | |}
| + | |
| − | |}
| + | |
| − | | + | |
| − | Podemos aproximar la derivada direccional de <math display="inline">\phi </math> en el vértice <math display="inline">[v_{1}]</math> en dirección del vector <math display="inline">[v_{2}] - [v_{1}]</math> como:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>d(\phi _{[v_{1}]})_{[v_{1},v_{2}]} = \phi _{2} - \phi _{1}, </math>
| + | |
| − | |}
| + | |
| − | |}
| + | |
| − | | + | |
| − | de manera similar, la derivada direccional de <math display="inline">\phi </math> en dirección del vector <math display="inline">[v_{3}] - [v_{1}]</math> está dada por:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>d(\phi _{[v_{1}]})_{[v_{1},v_{3}]} = \phi _{3} - \phi _{1}. </math>
| + | |
| − | |}
| + | |
| − | |}
| + | |
| − | | + | |
| − | Por lo tanto, para encontrar el vector gradiente <math display="inline">\underline{W}W_{1}</math> de <math display="inline">\phi </math> en el vértice <math display="inline">[v_{1}]</math>, podemos plantear el siguiente sistema de ecuaciones:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>\underline{W}W_{1} \cdot ([v_{2}] - [v_{1}]) = \phi _{2} - \phi _{1}, </math>
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math> \underline{W}W_{1} \cdot ([v_{3}] - [v_{1}]) = \phi _{3} - \phi _{1}, </math>
| + | |
| − | |}
| + | |
| − | |}
| + | |
| − | | + | |
| − | donde:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>[v_{1}] = (x_{1},y_{1}), </math>
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math> [v_{2}] = (x_{2},y_{2}), </math>
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math> [v_{3}] = (x_{3},y_{3}). </math>
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" |
| + | |
| − | |}
| + | |
| − | |}
| + | |
| − | | + | |
| − | Resolviendo el sistema anterior para <math display="inline">\underline{W}W_{1}</math> se obtiene:
| + | |
| − | | + | |
| − | <span id="eq-8"></span>
| + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>\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}, </math>
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;" | (8)
| + | |
| − | |}
| + | |
| − | | + | |
| − | donde <math display="inline">A</math> es el área del triangulo. Siguiendo el mismo procedimiento para encontrar <math display="inline">\underline{W}W_{2}</math> y <math display="inline">\underline{W}W_{3}</math>, se observa que:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>\underline{W}W_{1} = \underline{W}W_{2} = \underline{W}W_{3}. </math>
| + | |
| − | |}
| + | |
| − | |}
| + | |
| − | | + | |
| − | Más aún, el vector gradiente obtenido en la Ecuación ([[#eq-8|8]]) coincide con el que se obtiene con FEML. Consideremos ahora la velocidad de partícula <math display="inline">\underline{v}v</math> discretizada sobre los vértices del triángulo:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>\underline{u}u_{1} = (u_{1,1}, u_{1,2}), </math>
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math> \underline{u}u_{2} = (u_{2,1}, u_{2,2}), </math>
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math> \underline{u}u_{3} = (u_{3,1}, u_{3,2}). </math>
| + | |
| − | |}
| + | |
| − | |}
| + | |
| − | | + | |
| − | 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 ([[#eq-5|5]]):
| + | |
| − | | + | |
| − | <span id="eq-9"></span>
| + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>\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} , </math>
| + | |
| − | |}
| + | |
| − | | style="width: 5px;text-align: right;white-space: nowrap;" | (9)
| + | |
| − | |}
| + | |
| − | | + | |
| − | donde:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>V_{1,1} = u_{1,1}(y_{2} - y_{3}) + u_{1,2}(x_{3} - x_{2}), </math>
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math> V_{2,1} = u_{2,1}(y_{2} - y_{3}) + u_{2,2}(x_{3} - x_{2}), </math>
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math> V_{3,1} = u_{3,1}(y_{2} - y_{3}) + u_{3,2}(x_{3} - x_{2}),</math>
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math> V_{1,2} = u_{1,1}(y_{3} - y_{1}) + u_{1,2}(x_{1} - x_{3}), </math>
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math> V_{2,2} = u_{2,1}(y_{3} - y_{1}) + u_{2,2}(x_{1} - x_{3}), </math>
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math> V_{3,2} = u_{3,1}(y_{3} - y_{1}) + u_{3,2}(x_{1} - x_{3}), </math>
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math> V_{1,3} = u_{1,1}(y_{1} - y_{2}) + u_{1,2}(x_{2} - x_{1}), </math>
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math> V_{2,3} = u_{2,1}(y_{1} - y_{2}) + u_{2,2}(x_{2} - x_{1}), </math>
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math> V_{3,3} = u_{3,1}(y_{1} - y_{2}) + u_{3,2}(x_{2} - x_{1}). </math>
| + | |
| − | |}
| + | |
| − | |}
| + | |
| − | | + | |
| − | Finalmente, sustituyendo la Ecuación ([[#eq-9|9]]) en ([[#eq-7|7]]), obtenemos la discretización basada en DEC de la Ecuación de Transporte para flujo incompresible:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>[M_{0,2}V + kD_{0,1}^{T}M_{1,1}D_{0,1}]\phi = M_{0,2}q. </math>
| + | |
| − | |}
| + | |
| − | |}
| + | |
| − | | + | |
| − | ==4 Pruebas Numéricas==
| + | |
| − | | + | |
| − | Para las pruebas numéricas se consideró como dominio una elipse cuya ecuación está dada por:
| + | |
| − | | + | |
| − | {| class="formulaSCP" style="width: 100%; text-align: left;"
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | {| style="text-align: left; margin:auto;width: 100%;"
| + | |
| − | |-
| + | |
| − | | style="text-align: center;" | <math>x^{2} + 4y^{2}-100=0, </math>
| + | |
| − | |}
| + | |
| − | |}
| + | |
| − | | + | |
| − | con las siguientes características (ver Figura [[#img-4|4]]):
| + | |
| − | | + | |
| − | * Condiciones tipo Dirichlet <math display="inline">\phi = 0</math> en la frontera de la elipse.
| + | |
| − | * Fuentes <math display="inline">q=1</math> sobre dos discos de radio <math display="inline">0.5</math> y centro en <math display="inline">(-7.5;0)</math> y <math display="inline">(7.5;0)</math>, respectivamente.
| + | |
| − | * Coeficiente de difusión <math display="inline">k = 0.5</math>.
| + | |
| − | * Velocidad de partícula <math display="inline">\underline{v}v = (2y;-x)</math>.
| + | |
| − | | + | |
| − | <div id='img-4'></div>
| + | |
| − | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
| + | |
| − | |-
| + | |
| − | |[[Image:Draft_Herrera_100377898-EllipticDomain.png|480px|Elipse considerado para las pruebas numéricas]]
| + | |
| − | |- style="text-align: center; font-size: 75%;"
| + | |
| − | | colspan="1" | '''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 [[#img-6|6]].
| + | |
| − | | + | |
| − | <div id='img-5a'></div>
| + | |
| − | <div id='img-5b'></div>
| + | |
| − | <div id='img-5'></div>
| + | |
| − | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
| + | |
| − | |-
| + | |
| − | |[[Image:Draft_Herrera_100377898-mesh1C.png|600px|]]
| + | |
| − | |[[Image:Draft_Herrera_100377898-mesh2C.png|600px|]]
| + | |
| − | |- style="text-align: center; font-size: 75%;"
| + | |
| − | | (a)
| + | |
| − | | (b)
| + | |
| − | |- style="text-align: center; font-size: 75%;"
| + | |
| − | | colspan="2" | '''Figura 5'''
| + | |
| − | |}
| + | |
| − | <div id='img-6a'></div>
| + | |
| − | <div id='img-6b'></div>
| + | |
| − | <div id='img-6'></div>
| + | |
| − | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
| + | |
| − | |-
| + | |
| − | |[[Image:Draft_Herrera_100377898-mesh3C.png|600px|]]
| + | |
| − | |[[Image:Draft_Herrera_100377898-mesh4C.png|600px|]]
| + | |
| − | |- style="text-align: center; font-size: 75%;"
| + | |
| − | | (a)
| + | |
| − | | (b)
| + | |
| − | |- style="text-align: center; font-size: 75%;"
| + | |
| − | | colspan="2" | '''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 [[#img-6b|6b]]) se muestra en la Figura [[#img-7|7]]. Para fines de comparación, el resultado obtenido con FEML empleando la malla más fina (Figura [[#img-6b|6b]]) es tomado como la mejor aproximación de la solución analítica. El Cuadro [[#table-1|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.
| + | |
| − | | + | |
| − | <div id='img-7'></div>
| + | |
| − | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
| + | |
| − | |-
| + | |
| − | |[[Image:Draft_Herrera_100377898-analiticsC.png|480px|Resultado obtenido con '''DEC-based''' considerando la malla más fina (19679 nodos y 38388 elementos.)]]
| + | |
| − | |- style="text-align: center; font-size: 75%;"
| + | |
| − | | colspan="1" | '''Figura 7:''' Resultado obtenido con '''DEC-based''' considerando la malla más fina (<math>19679</math> nodos y <math>38388</math> elementos.)
| + | |
| − | |}
| + | |
| − | | + | |
| − | | + | |
| − | {| class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
| + | |
| − | |+ style="font-size: 75%;" |<span id='table-1'></span>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.
| + | |
| − | |-
| + | |
| − | |
| + | |
| − | |
| + | |
| − | |
| + | |
| − | | colspan='2' style="border-left: 2px solid;border-right: 2px solid;" | Valores Máximos
| + | |
| − | |- style="border-top: 2px solid;"
| + | |
| − | | colspan='1' style="border-left: 2px solid;" | Malla
| + | |
| − | | Nodos
| + | |
| − | | Elementos
| + | |
| − | | FEML
| + | |
| − | | DEC-b
| + | |
| − | | colspan='1' style="border-right: 2px solid;" | DSP
| + | |
| − | |- style="border-top: 2px solid;"
| + | |
| − | | colspan='1' style="border-left: 2px solid;" | Figura [[#img-5a|5a]]
| + | |
| − | | <math>491</math>
| + | |
| − | | <math>860</math>
| + | |
| − | | <math>0.10852</math>
| + | |
| − | | colspan='1' | <math>0.10714</math>
| + | |
| − | | colspan='1' style="border-right: 2px solid;" | <math>0.01571</math>
| + | |
| − | |- style="border-top: 2px solid;"
| + | |
| − | | colspan='1' style="border-left: 2px solid;" | Figura [[#img-5|5]]
| + | |
| − | | <math>983</math>
| + | |
| − | | <math>1772</math>
| + | |
| − | | <math>0.11105</math>
| + | |
| − | | colspan='1' | <math>0.10983</math>
| + | |
| − | | colspan='1' style="border-right: 2px solid;" | <math>0.01302</math>
| + | |
| − | |- style="border-top: 2px solid;"
| + | |
| − | | colspan='1' style="border-left: 2px solid;" | Figura [[#img-6a|6a]]
| + | |
| − | | <math>5246</math>
| + | |
| − | | <math>10010</math>
| + | |
| − | | <math>0.13583</math>
| + | |
| − | | colspan='1' | <math>0.13653</math>
| + | |
| − | | colspan='1' style="border-right: 2px solid;" | <math>0.01368</math>
| + | |
| − | |- style="border-top: 2px solid;border-bottom: 2px solid;"
| + | |
| − | | colspan='1' style="border-left: 2px solid;" | Figura [[#img-6b|6b]]
| + | |
| − | | <math>19679</math>
| + | |
| − | | <math>38388</math>
| + | |
| − | | <math>0.12285</math>
| + | |
| − | | colspan='1' | <math>0.12295</math>
| + | |
| − | | colspan='1' style="border-right: 2px solid;" | <math>1 \times 10^{-4}</math>
| + | |
| − | | + | |
| − | |}
| + | |
| − | | + | |
| − | La Figura [[#img-9|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 [[#img-6b|6b]].
| + | |
| − | | + | |
| − | El problema presenta inestabilidades para ambos métodos cuando se emplean las mallas mostradas en las Figuras [[#img-5a|5a]] y [[#img-5|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.
| + | |
| − | | + | |
| − | <div id='img-8a'></div>
| + | |
| − | <div id='img-8b'></div>
| + | |
| − | <div id='img-8'></div>
| + | |
| − | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
| + | |
| − | |-
| + | |
| − | |[[Image:Draft_Herrera_100377898-mesh1C.png|600px|]]
| + | |
| − | |[[Image:Draft_Herrera_100377898-mesh2C.png|600px|]]
| + | |
| − | |- style="text-align: center; font-size: 75%;"
| + | |
| − | | (a)
| + | |
| − | | (b)
| + | |
| − | |- style="text-align: center; font-size: 75%;"
| + | |
| − | | colspan="2" | '''Figura 8'''
| + | |
| − | |}
| + | |
| − | <div id='img-9a'></div>
| + | |
| − | <div id='img-9b'></div>
| + | |
| − | <div id='img-9'></div>
| + | |
| − | {| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
| + | |
| − | |-
| + | |
| − | |[[Image:Draft_Herrera_100377898-mesh3C.png|600px|]]
| + | |
| − | |[[Image:Draft_Herrera_100377898-mesh4C.png|600px|]]
| + | |
| − | |- style="text-align: center; font-size: 75%;"
| + | |
| − | | (a)
| + | |
| − | | (b)
| + | |
| − | |- style="text-align: center; font-size: 75%;"
| + | |
| − | | colspan="2" | '''Figura 9:''' Valores de temperatura a lo largo de una sección sobre el eje mayor de la elipse para diferentes mallas. <math>(a)</math> Gráfica para la malla mostrada en la Figura [[#img-5a|5a]]; <math>(b)</math> Gráfica para la malla mostrada en la Figura [[#img-5|5]]; <math>(c)</math> Gráfica para la malla mostrada en la Figura [[#img-6a|6a]]; <math>(d)</math> 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===
| + | |
| − | | + | |
| − | <div id="cite-1"></div>
| + | |
| − | '''[[#citeF-1|[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.
| + | |
| − | | + | |
| − | <div id="cite-2"></div>
| + | |
| − | '''[[#citeF-2|[2]]]''' Keenan Crane y col. ''Digital Geometry Processing with Discrete Exterior Calculus''. En: SIG- GRAPH ’13 (2013).
| + | |
| − | | + | |
| − | <div id="cite-3"></div>
| + | |
| − | '''[[#citeF-3|[3]]]''' Keenan Crane. ''Discrete Differential Geometry'' (2019) https://www.cs.cmu.edu/ kmcrane/Projects/DDG/paper.pdf
| + | |
| − | | + | |
| − | <div id="cite-4"></div>
| + | |
| − | '''[[#citeF-4|[4]]]''' Andrew Gillette. ''Notes on Discrete Exterior Calculus'' (2009) https://www.math.arizona.edu/ agillette/research/decNotes.pdf
| + | |
| − | | + | |
| − | <div id="cite-5"></div>
| + | |
| − | '''[[#citeF-5|[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
| + | |
| − | | + | |
| − | <div id="cite-6"></div>
| + | |
| − | '''[[#citeF-6|[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.
| + | |
| − | | + | |
| − | <div id="cite-7"></div>
| + | |
| − | '''[[#citeF-7|[7]]]''' Anil Nirmal Hirani. ''Discrete Exterior Calculus''. Tesis doct. California Institute of Technology, 2003.
| + | |
| − | | + | |
| − | <div id="cite-8"></div>
| + | |
| − | '''[[#citeF-8|[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.
| + | |
| − | | + | |
| − | <div id="cite-9"></div>
| + | |
| − | '''[[#citeF-9|[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.
| + | |
| − | | + | |
| − | <div id="cite-10"></div>
| + | |
| − | '''[[#citeF-10|[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.
| + | |
