​

En este artículo se presenta un método de alto orden de Galerkin discontinuo para problemas de flujo compresible, en los cuales es muy frecuente la aparición de choques. La estabilización se introduce mediante una nueva base de funciones. Esta base tiene la flexibilidad de variar localmente (en cada elemento) entre un espacio de funciones polinómicas continuas o un espacio de funciones polinómicas a trozos. Así, el método propuesto proporciona un puente entre los métodos estándar de alto orden de Galerkin discontinuo y los clásicos métodos de volúmenes finitos, manteniendo la localidad y compacidad del esquema. La variación de las funciones de la base se define automáticamente en función de la regularidad de la solución y la estabilización se introduce mediante el operador salto, estándar en los métodos Galerkin discontinuo. A diferencia de los clásicos métodos de limitadores de pendiente, la estrategia que se presenta es muy local y robusta, y es aplicable a cualquier orden de aproximación. Además, el método propuesto no requiere refinamiento adaptativo de la malla ni restricción del esquema de integración temporal. Se consideran varias aplicaciones de las ecuaciones de Euler que demuestran la validez y efectividad del método, especialmente para altos órdenes de aproximación.

​

==Abstract==

​

This article presents a high-order Discontinuous Galerkin method for compressible flow problems, in which is very frequent the formation of shocks. The stabilization is introduced by a new basis functions. This base has the flexibility to vary locally (within each element) between continuous polynomial functions space and a space of piecewise polynomial functions. Thus, the proposed method provides a bridge between the standard methods of high-order Discontinuous Galerkin and classical Finite Volume methods, maintaining the locality and compactness of the scheme. The variation of basis functions is automatically set according to the regularity of the solution and the stabilization is introduced by the jump operator, standard in Discontinuous Galerkin methods. Unlike the classical methods of slope limiting, the strategy here presented is very local, robust, and applies to any order of approximation. Moreover, the proposed method does not require adaptive mesh refinement techniques and it can be used with any temporal integration scheme. Several applications of the Euler equations are shown, demonstrating the validity and effectiveness of the method, especially for high orders of approximation.

​

==Palabras clave==

​

Galerkin discontinuo de alto orden ; Ecuaciones de Euler ; Captura de choques ; Funciones de forma ; Sensor

​

==Keywords==

​

High-order discontinous Galerkin ; Euler equations ; Shock-capturing ; Shape functions ; Sensor

​

==1. Introducción==

​

El desarrollo de métodos numéricos cada vez más precisos para la resolución de problemas de dinámica de fluidos es una imperante necesidad. Dos de los campos de la dinámica de fluidos computacional que requieren esquemas de máxima precisión son la turbulencia y la aeroacústica. Las ecuaciones involucradas son las de Navier-Stokes y las de Euler. En la última década, el desarrollo de métodos numéricos de alta precisión ha ido ligado, en general, al aumento del orden de los esquemas numéricos en mallas no estructuradas. Sin embargo, los esquemas de alto orden para el cálculo de flujo compresible son computacionalmente inestables debido a la aparición de oscilaciones espurias físicamente inaceptables [[#bib0005|[1]]]  and [[#bib0010|[2]]] .      

​

Los métodos de Galerkin discontinuo [[#bib0015|[3]]] , [[#bib0020|[4]]] , [[#bib0025|[5]]] , [[#bib0030|[6]]]  and [[#bib0035|[7]]]  se presentan como una buena alternativa para la solución de problemas de propagación de ondas y, en general, de leyes de conservación hiperbólicas. Estos métodos proporcionan un puente entre la aproximación mediante funciones de forma de elementos finitos continuos y la tecnología de los métodos de volúmenes finitos [[#bib0040|[8]]]  and [[#bib0045|[9]]] . Para problemas puramente hiperbólicos los métodos Galerkin discontinuo proporcionan un esquema muy compacto para cualquier orden de aproximación y cualquier tipo de malla no estructurada e incluso con nodos colgantes. Lamentablemente, en presencia de choques y cuando se utilizan altos órdenes de aproximación (''p''  ≥ 1) el mecanismo estabilizador introducido por los métodos Galerkin discontinuo mediante el operador salto no es suficiente. En estos casos es necesario utilizar otras técnicas más robustas de estabilización que, en general, van acompañadas del refinamiento adaptativo de la malla.      

​

Entre las técnicas posibles de resolución de choques destacan las técnicas de limitación [[#bib0025|[5]]] , [[#bib0050|[10]]] , [[#bib0055|[11]]] , [[#bib0060|[12]]]  and [[#bib0065|[13]]]  y las técnicas de difusión artificial [[#bib0070|[14]]] , [[#bib0075|[15]]]  and [[#bib0080|[16]]] . Los limitadores de pendiente se introducen inicialmente para los métodos de volúmenes finitos [[#bib0010|[2]]]  y posteriormente se extienden a métodos de Galerkin discontinuo [[#bib0020|[4]]]  and [[#bib0085|[17]]] . Estos esquemas buscan la monotonicidad de la solución en los valores promedios y, por lo tanto, reducen la aproximación a orden constante o, en el mejor de los casos, lineal, no solo en la región donde se localiza el choque, sino también en su entorno. Además, para poder garantizar estabilidad no lineal en una seminorma arbitraria solo pueden emplearse en combinación con una clase particular de esquemas de integración temporal Runge-Kutta explícitos que hasta el momento, según el conocimiento de los autores, como mucho son de cuarto orden [[#bib0090|[18]]] , [[#bib0095|[19]]]  and [[#bib0100|[20]]] . En resumen, las técnicas de limitadores de pendiente añaden difusión artificial de orden ''h''  (siendo ''h''  tamaño de la malla) proporcionando una solución estable, pero a costa de perder el alto orden de precisión. Esto hace necesario el uso de técnicas de refinamiento adaptativo, que además deben tener en cuenta el carácter direccional del choque.      

​

Otras técnicas de reconstrucción más sofisticadas que han centrado la atención de varios autores son los métodos de reconstrucción no oscilatorios de alto orden ENO y WENO [[#bib0095|[19]]] , [[#bib0105|[21]]]  and [[#bib0110|[22]]] . Estos métodos preservan la estabilidad no lineal del esquema así como el alto orden de aproximación proporcionando soluciones libres de oscilaciones. No obstante, añaden grados de libertad adicionales y pierden la compacidad del esquema, debido al uso de stencils amplios. Además el uso de estos métodos queda esencialmente restringido a mallas uniformes y estructuradas, raramente empleadas en el tratamiendo de choques y problemas complejos.      

​

La segunda alternativa consiste en la introducción de difusión artificial para obtener soluciones estables. Esta técnica ya fue introducida en los años cincuenta por Von Neumann y Richtmyer [[#bib0070|[14]]] . En su extensión a los métodos Galerkin discontinuo es habitual definir la constante en cada elemento, de esta forma se obtiene una solución precisa que mantiene el grado de la aproximación. En general, la solución es de orden superior al de los limitadores de pendiente. Recientemente se han presentado nuevos enfoques donde el orden obtenido es ''h'' /''p''  o ''h''<sup>''k''</sup>  con ''k''  ≤ ''p'' , siendo ''p''  el grado de aproximación  [[#bib0075|[15]]]  and [[#bib0080|[16]]] . Sin embargo, estos métodos presentan serias dificultades en su extensión al caso multidimensional. Por un lado, ajustar la magnitud de la difusión no es sencillo, así como tampoco el efecto de propagación multidireccional del choque. Por otro lado, añadir difusión constante en todo elemento puede derivar en una pérdida de resolución significativa cuando se utilizan mallas groseras en presencia de capas límite o estructuras no dimensionales, como los choques o discontinuidades. Además, la introducción de difusión artificial implica la resolución de una ecuación diferencial parcial de segundo orden, lo cual conlleva un error asociado a la ecuación original, así como un incremento del coste computacional  [[#bib0115|[23]]]  and [[#bib0120|[24]]] .      

​

En este trabajo se explota la flexibilidad que ofrecen los métodos de Galerkin discontinuo, variando el espacio de aproximación en cada elemento, así como también para distintos instantes de la simulación. Un espacio de funciones continuas debe emplearse para aproximar soluciones suaves, y un espacio de funciones discontinuas a trozos es apropiado para soluciones con discontinuidades o gradientes pronunciados. Basándose en esta idea, se introduce una nueva base de funciones polinómicas que dependen linealmente de un parámetro. Este parámetro permite definir funciones de la base localmente según la regularidad de la solución. Así, para un valor del parámetro la base permite una representación continua de alto orden de la solución. Variando el valor del parámetro se introduce una representación local, mediante una expansión polinómica a trozos en el interior de cada elemento. El objetivo reside en aprovechar la potente tecnología que proporcionan los métodos de volúmenes finitos en el contexto de leyes hiperbólicas de conservación, y en general, para problemas de convección dominante. Una de las principales ventajas del método propuesto respecto a los métodos de difusión artificial dependientes de uno o más parámetros es que en este caso el parámetro que define la base de funciones no es arbitrario y está controlado por un indicador de regularidad.

​

El método propuesto utiliza el principio estabilizador de los métodos de Galerkin discontinuo en la interfaz de los elementos mediante el operador salto y lo extiende al interior de ellos, tomando como idea un argumento similar a las técnicas de volúmenes finitos. A diferencia de estos últimos, el método propuesto tiene la propiedad de calibrar la magnitud de las discontinuidades introducidas y, por lo tanto, de controlar la difusión en el interior del elemento. Adicionalmente, la nueva base de funciones se define únicamente en el dominio espacial, y por lo tanto no interfiere en el esquema de integración temporal empleado, permitiendo el uso de esquemas explícitos o implícitos de cualquier orden.

​

El artículo se divide según sigue: en la sección [[#sec0010|2]]  se resume brevemente el método de Galerkin discontinuo estándar para funciones de forma hiperbólicas, en particular para las ecuaciones de Euler. En la sección [[#sec0015|3]]  se define la nueva base de funciones de forma y el criterio de regularidad impuesto. También se detalla el sensor de discontinuidades utilizado. La sección [[#sec0035|4]]  recoge varios tests numéricos que demuestran la aplicabilidad del método para diversas situaciones frente a otros métodos clásicos, así como su eficiencia para altos órdenes de aproximación. Finalmente en la sección [[#sec0055|5]]  se resumen las conclusiones obtenidas.      

​

<span id='sec0010'></span>

==2. El método de Galerkin discontinuo para las ecuaciones de Euler==

​

En esta sección se resume brevemente el método de Galerkin discontinuo para las ecuaciones de Euler. Sin pérdida de generalidad se asume el caso bidimensional <math display="inline">(x_1,x_2)\in {\mathbb{R}}^2</math> .      

​

Las ecuaciones de Euler de la dinámica de gases computacional [[#bib0005|[1]]]  expresan la conservación de la masa, el momento y la energía para un fluido ideal compresible y no viscoso. En ausencia de fuerzas externas, las ecuaciones de Euler se expresan de forma conservativa como el siguiente sistema

​

<span id='eq0005'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\frac{\partial }{\partial t}\left(\begin{array}{c}

\rho \\

\rho u\\

\rho v\\

\rho E\\

​

\end{array}\right)+\frac{\partial }{\partial x_1}\left(\begin{array}{c}

\rho u\\

\rho u^2+p\\

\rho uv\\

u(\rho E+p)

\end{array}\right)+\frac{\partial }{\partial x_2}\left(\begin{array}{c}

\rho v\\

\rho uv\\

\rho v^2+p\\

v(\rho E+p)

\end{array}\right)=\boldsymbol{0}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 1)

|}

​

definido en un dominio acotado <math display="inline">\Omega \subset {\mathbb{R}}^2</math>  con condiciones de contorno en ∂Ω. En esta expresión ''ρ''  representa la densidad del fluido, ''u  ''  y <math display="inline">v</math>  son las componentes del vector velocidad <math display="inline">\boldsymbol{v}</math> , ''p''  es la presión y ''E''  es la energía total. Para completar el sistema se define una ecuación de estado, que relaciona la energía interna con la presión y la densidad. Para un gas perfecto politrópico la ecuación de estado es

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>p=(\gamma -1)(E-\frac{\rho \vert \vert \boldsymbol{v}\vert \vert ^2}{2})</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | 

|}

​

donde ''γ''  es el exponente adiabático, con valor ''γ''  = 1,4 para el aire.      

​

De forma compacta el sistema [[#eq0005|(1)]]  se expresa como

​

<span id='eq0015'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\frac{\partial \boldsymbol{\mbox{U}}}{\partial t}+</math><math>\frac{\partial \boldsymbol{\mbox{F}}_k(\boldsymbol{\mbox{U}})}{\partial x_k}=</math><math>\boldsymbol{0}\quad para\quad k=1\mbox{,}2</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 2)

|}

​

donde se utiliza el criterio de Einstein para la notación, con

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{U}}=\left(\begin{array}{c}

\rho \\

\rho u\\

\rho v\\

\rho E\\

​

\end{array}\right)\quad \boldsymbol{\mbox{F}}_\boldsymbol{1}=</math><math>\left(\begin{array}{c}

\rho u\\

\rho u^2+p\\

\rho uv\\

u(\rho E+p)

\end{array}\right)\quad y\quad \boldsymbol{\mbox{F}}_\boldsymbol{2}=</math><math>\left(\begin{array}{c}

\rho v\\

\rho uv\\

\rho v^2+p\\

v(\rho E+p)

\end{array}\right)</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | 

|}

​

En el desarrollo posterior se emplea el número de Mach, definido como <math display="inline">M=\frac{\parallel \boldsymbol{v}\parallel }{c}</math>  donde <math display="inline">c=\sqrt{\gamma p/\rho }</math>  es la velocidad del sonido.      

​

El método de Galerkin discontinuo para un sistema hiperbólico de leyes de conservación se describe como sigue. Se divide el dominio computacional Ω en una partición de elementos no superpuestos entre sí:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\overline{\Omega }={\bigcup }_{e=1}^{n_{el}}{\overline{\Omega }}_e\quad tal\quad que\quad {\Omega }_i\bigcap {\Omega }_j=</math><math>\varnothing \quad para\quad i\not =j</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 3)

|}

​

donde n<sub>el</sub>  es el número total de elementos de la partición.      

​

Sea <math display="inline">\boldsymbol{\mbox{W}}\in [P^p({\Omega }_e)]^{n_{sd}+2}</math>  una función de test, donde n<sub>sd</sub>  respresenta la dimensión del espacio vectorial de trabajo (en este artículo se trabaja con n<sub>sd</sub>  = 2). Además, <math display="inline">P^p({\Omega }_e)</math>  es el conjunto de polinomios de grado menor o igual a ''p''  en Ω<sub>''e''</sub> . Multiplicando [[#eq0015|(2)]]  por la función de test, integrando sobre el elemento Ω<sub>''e''</sub>  y usando el teorema de la divergencia se obtiene la forma débil

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\int }_{{\Omega }_e}\boldsymbol{\mbox{W}}\cdot \frac{\partial \boldsymbol{\mbox{U}}_e}{\partial t}\quad d\Omega -</math><math>{\int }_{{\Omega }_e}\frac{\partial \boldsymbol{\mbox{W}}}{\partial x_k}\cdot \boldsymbol{\mbox{F}}_k(\boldsymbol{\mbox{U}}_e)\quad d\Omega +</math><math>{\int }_{\partial {\Omega }_e}\boldsymbol{\mbox{W}}\cdot \boldsymbol{\mbox{F}}_{\boldsymbol{\mbox{n}}_e}(\boldsymbol{\mbox{U}}_e)\quad d\Gamma =</math><math>0\quad \forall \boldsymbol{\mbox{W}}\in [P^p({\Omega }_e)]^{n_{sd}+2}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 4)

|}

​

donde '''U'''<sub>''e''</sub>  es la restricción de '''U'''  en el elemento Ω<sub>''e''</sub>  y <math display="inline">\boldsymbol{\mbox{F}}_{\boldsymbol{\mbox{n}}_e}(\boldsymbol{\mbox{U}}_e)</math>  es la función de flujo normal. Como consecuencia de la discontinuidad, la función de flujo normal <math display="inline">\boldsymbol{\mbox{F}}_{\boldsymbol{\mbox{n}}_e}</math>  tiene 2 valores en la interfaz de los elementos. Por ello se introduce un flujo numérico <math display="inline">{\overset{\mbox{ˆ}}{\boldsymbol{\mbox{F}}}}_{\boldsymbol{\mbox{n}}_e}(\boldsymbol{\mbox{U}}_e,\boldsymbol{\mbox{U}}_e^{out})</math>  que depende de los estados de la variable '''U'''  a ambos lados de cada interfaz entre elementos, y del vector unitario normal a esta, '''n'''<sub>''e''</sub> :

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{U}}_e^{out}=\underset{\varepsilon \rightarrow 0^+}{l\acute{\imath }m}\boldsymbol{\mbox{U}}(\boldsymbol{x}+</math><math>\varepsilon \boldsymbol{\mbox{n}}_e)\quad para\quad \boldsymbol{x}\in {\Omega }_e</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | 

|}

​

Para el flujo numérico <math display="inline">{\boldsymbol{\overset{\mbox{ˆ}}{\mbox{F}}}}_{\boldsymbol{\mbox{n}}_e}</math>  es posible emplear cualquiera de los Riemann ''solvers''  estándar de los métodos de volúmenes finitos  [[#bib0045|[9]]] . En los ejemplos numéricos que se presentan en este artículo se ha elegido el flujo de Roe [[#bib0125|[25]]] . Así, la forma débil del método Galerkin discontinuo para cada elemento Ω<sub>''e''</sub>  se escribe como

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\int }_{{\Omega }_e}\boldsymbol{\mbox{W}}\cdot \frac{\partial \boldsymbol{\mbox{U}}_e}{\partial t}\quad d\Omega -</math><math>{\int }_{{\Omega }_e}\frac{\partial \boldsymbol{\mbox{W}}}{\partial x_k}\cdot \boldsymbol{\mbox{F}}_k(\boldsymbol{\mbox{U}}_e)\quad d\Omega +</math><math>{\int }_{\partial {\Omega }_e}\boldsymbol{\mbox{W}}\cdot {\overset{\mbox{ˆ}}{\boldsymbol{\mbox{F}}}}_{\boldsymbol{\mbox{n}}_e}(\boldsymbol{\mbox{U}}_e,\boldsymbol{\mbox{U}}_e^{out})\quad d\Gamma =</math><math>0\quad \forall \boldsymbol{\mbox{W}}\in [P^p({\Omega }_e)]^{n_{sd}+2}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 5)

|}

​

Discretizando la forma débil en cada elemento se obtiene el sistema de ecuaciones diferenciales ordinarias

​

<span id='eq0045'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{M}}\frac{d\boldsymbol{U}}{dt}+\boldsymbol{\mbox{R}}(\boldsymbol{U})=</math><math>\boldsymbol{0}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 6)

|}

​

donde en este caso '''''U'''''  representa el vector de coeficientes de la aproximación, '''M'''  es la matriz de masa, que resulta diagonal por bloques, y '''R''' ('''''U''''' ) es el vector residuo. El sistema [[#eq0045|(6)]]  puede resolverse con diversos esquemas de integración temporal. En el presente trabajo se utiliza un método estándar Runge-Kutta explícito [[#bib0005|[1]]] . La condición de estabilidad viene determinada por una acotación del número de Courant:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>C=\vert {\lambda }_{m\acute{a}x}\vert \frac{\Delta t}{h}\leq \frac{1}{2p+1}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | 

|}

​

donde ''h''  es el tamaño elemental y ''p  ''  es el grado de la aproximación. El valor escalar <math display="inline">{\lambda }_{m\acute{a}x}</math>  es el valor propio máximo en valor absoluto de la matriz de Jacobi del flujo, <math display="inline">\frac{\partial \boldsymbol{\mbox{F}}}{\partial \boldsymbol{\mbox{U}}}</math> , es decir, <math display="inline">{\lambda }_{m\acute{a}x}=</math><math>{m\acute{a}x}_j\vert {\lambda }_j\vert </math>  para ''j''  = 1, …, n<sub>sd</sub>  + 2.

​

==Observación 2.1.                     ==

​

En el caso de un sistema hiperbólico de ecuaciones, como las ecuaciones de Euler, los valores propios ''λ''<sub>''j''</sub>  no son más que las velocidades de propagación de las cantidades características. Para más detalles, consultar  [[#bib0005|[1]]]  and [[#bib0130|[26]]] .

​

<span id='sec0015'></span>

==3. Base de funciones==

​

El método propuesto se caracteriza por la introducción de la nueva base de funciones polinómicas <math display="inline">{\overline{N}}_i</math> . Como es habitual en los métodos de Galerkin discontinuo, se aproxima la solución '''U'''  elemento a elemento con una base de funciones de grado p:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{U}}_e=\sum_{i=1}^{n_{en}(p)}{\overline{N}}_i(\boldsymbol{x};\alpha )\boldsymbol{U}_i\quad para\quad e=</math><math>1\ldots n_{el}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 7)

|}

​

donde n<sub>en</sub> (''p'' ) es el número de nodos en cada elemento (es decir, el número de grados de libertad por elemento).      

​

La base de funciones <math display="inline">\overline{N}(\boldsymbol{x};\alpha )</math>  tiene un papel crucial en la precisión de la aproximación. Las funciones <math display="inline">{\overline{N}}_i</math>  se definen en cada elemento como una combinación convexa de las funciones de forma continuas estándar de los métodos Galerkin discontinuo ''N''<sub>''i''</sub> , y unas funciones de forma constantes a trozos en el elemento, ''ϕ''<sub>''i''</sub> . Es decir,

​

<span id='eq0060'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\overline{N}}_i(\boldsymbol{x};\alpha ):=\alpha N_i(\boldsymbol{x})+</math><math>(1-\alpha ){\phi }_i(\boldsymbol{x})\quad para\quad i=</math><math>1,\ldots ,n_{en}(p)</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 8)

|}

​

donde ''α''  es un parámetro que toma valores entre 0 y 1 en función de la suavidad de la aproximación.      

​

Una de las ventajas de los métodos de Galerkin discontinuo reside en su inherente mecanismo estabilizador mediante la función salto <math display="inline">\mbox{〚 } \boldsymbol{\mbox{U}}\cdot \boldsymbol{\mbox{n}}_e\mbox{〛} =</math><math>\boldsymbol{\mbox{U}}_e\cdot \boldsymbol{\mbox{n}}_e+</math><math>\boldsymbol{\mbox{U}}_e^{out}\cdot \boldsymbol{\mbox{n}}_{out}</math>  para aproximaciones constantes y lineales, véase por ejemplo [[#bib0030|[6]]] , [[#bib0060|[12]]]  and [[#bib0085|[17]]] . La base de funciones propuesta [[#eq0060|(8)]]  explota esta propiedad y la extiende al interior de los elementos, permitiendo así discontinuidades no solo en la interfaz, sino también en el interior del elemento. Asimismo, según la definición [[#eq0060|(8)]]  para ''α''  = 0 se obtiene la base de funciones discontinuas a trozos (<math display="inline">{\overline{N}}_i={\phi }_i</math> ) y para ''α''  = 1 se obtiene la base de funciones continuas (<math display="inline">{\overline{N}}_i=N_i</math> ), recuperando en este caso la forma estándar del método de Galerkin discontinuo.      

​

La [[#fig0005|figura 1]]  muestra una función <math display="inline">{\overline{N}}_i(\boldsymbol{x};\alpha )</math>  de grado ''p''  = 3 para distintos valores de ''α  '' . Los casos particulares <math display="inline">{\overline{N}}_i(\boldsymbol{x};\alpha )=</math><math>{\phi }_i</math>  y <math display="inline">{\overline{N}}_i(\boldsymbol{x};\alpha )=</math><math>N_i</math>  corresponden a las [[#fig0005|figura 1]] (a) y (b) respectivamente.

​

<span id='fig0005'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_703683201-1-s2.0-S0213131512000491-gr1.jpg|center|368px|Función de forma asociada a uno de los nodos arista de la base de funciones de ...]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 1.

​

Función de forma asociada a uno de los nodos arista de la base de funciones de grado ''p''  = 3 para distintos valores de ''α'' .                  

​

</span>

|}

​

La nueva base de funciones debe ser capaz de representar cualquier perfil de la solución, ya sea suave o con tendencia a formar frentes o discontinuidades. Por lo tanto, es de gran importancia elegir de manera adecuada, tanto el parámetro ''α''  como la base de funciones discontinuas ''ϕ'' . En los siguientes apartados se detallan ambas elecciones.      

​

===3.1. Definición de la base de funciones discontinuas===

​

La definición de la base de funciones discontinuas ''ϕ''<sub>''i''</sub>  tiene una estrecha relación con los métodos de volúmenes finitos  [[#bib0045|[9]]] , [[#bib0135|[27]]]  and [[#bib0140|[28]]] . En estos métodos se discretiza el dominio computacional en un conjunto de volúmenes de control (o celdas) no superpuestos entre sí y se impone la conservación de las ecuaciones en cada una de ellas. Siguiendo esta idea, se define la base ''ϕ''  de la siguiente manera.      

​

Se considera la partición arbitraria del elemento Ω<sub>''e''</sub>  en un conjunto de n<sub>en</sub> (''p'' ) volúmenes de control que no se solapan entre sí,

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\overline{\Omega }}_e={\bigcup }_k^{n_{en}(p)}{\overline{\Omega }}_e^k\quad tal\quad que\quad {\Omega }_e^l\bigcap {\Omega }_e^m=</math><math>\varnothing \quad para\quad l\not =m</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | 

|}

​

donde cada celda <math display="inline">{\Omega }_e^k</math>  contiene un nodo del elemento. La [[#fig0010|figura 2]]  muestra la partición para elementos de primer y segundo orden, donde se puede ver que las interfases quedan de forma arbitraria en una malla no estructurada.

​

<span id='fig0010'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_703683201-1-s2.0-S0213131512000491-gr2.jpg|center|366px|Partición en volúmenes de control de un elemento de primer y segundo orden. Los ...]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 2.

​

Partición en volúmenes de control de un elemento de primer y segundo orden. Los nodos del elemento están representados con puntos y los vértices de las celdas de control con cuadrados.

​

</span>

|}

​

==Observación 3.1.                     ==

​

La partición del elemento Ω<sub>''e''</sub>  en volúmenes de control no es única. Por simplicidad, en este trabajo se ha construido a partir de una subtriangulación del elemento que se obtiene uniendo los nodos. Se considera entonces el centroide de cada subtriángulo y los puntos medios de sus aristas. Uniendo el conjunto de puntos que forman de manera que cada polígono contenga uno y solo uno de los nodos elementales se obtiene una partición del elemento en n<sub>en</sub> (''p'' ) celdas.

​

Dada esta partición, cada función ''ϕ''<sub>''i''</sub>  se define como una función real <math display="inline">{\phi }_i(\boldsymbol{x}):{\Omega }_e\longrightarrow \mathbb{R}</math>  constante a trozos en cada volumen de control <math display="inline">{\Omega }_e^k</math> . Es decir, ∀'''''x'''''  ∈ Ω<sub>''e''</sub>

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\phi }_i(\boldsymbol{x})=\begin{array}{ll}

{\phi }_i^k=cst & si\quad \boldsymbol{x}\in {\Omega }_e^k\\

0 & si\quad \boldsymbol{x}\not\in {\Omega }_e^k\\

 & 

\end{array}\quad \forall \boldsymbol{x}\quad para\quad k=</math><math>1,\ldots ,n_{en}(p)</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 9)

|}

​

Por construcción, las funciones de forma continuas ''N''<sub>''i''</sub>  verifican la condición de reproducibilidad constante, o sea:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\sum_{i=1}^{n_{en}(p)}N_i(\boldsymbol{x})=1</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 10)

|}

​

Para que la nueva base de funciones <math display="inline">{\overline{N}}_i</math>  verifique también esta propiedad debe cumplirse

​

<span id='eq0080'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\sum_{i=1}^{n_{en}(p)}{\overline{N}}_i(\boldsymbol{x};\alpha )=</math><math>\alpha +(1-\alpha )\sum_{i=1}^{n_{en}(p)}{\phi }_i(\boldsymbol{x})=</math><math>1</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 11)

|}

​

o equivalentemente, <math display="inline">{\sum }_{i=1}^{n_{en}(p)}{\phi }_i^k=</math><math>1</math> .      

​

Existen varias opciones para definir <math display="inline">{\phi }_i^k</math>  de manera que se satisfaga [[#eq0080|(11)]] . La primera y más obvia opción es considerar <math display="inline">{\phi }_i^k=\frac{1}{n_{en}(p)}\quad \forall \quad k=</math><math>1,\ldots ,n_{en}(p)</math> . Esta opción lleva a una distribución constante y continua en cada elemento, sin ninguna aportación relevante a las funciones <math display="inline">{\overline{N}}_i(\boldsymbol{x};\alpha )</math> .      

​

Una segunda opción consiste en definir <math display="inline">{\phi }_i^k={\delta }_{ik}</math> . Esta opción proporciona una aproximación simple y eficaz, pero se descarta debido a que no tiene en cuenta el tamaño de la celda <math display="inline">{\Omega }_e^k</math>  correspodiente.      

​

Así pues, con el fin de establecer una relación entre el tamaño de las celdas y el valor de la función de forma asociada, se define <math display="inline">{\phi }_i^k</math>  como

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>{\phi }_i^k=\frac{1}{meas({\Omega }_e^k)}{\int }_{{\Omega }_e^k}N_i(\boldsymbol{x})\quad d\Omega </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 12)

|}

​

donde <math display="inline">{\Omega }_e^k</math>  es el subelemento que contiene el nodo '''''x'''''<sub>''i''</sub> . De este modo, se prioriza la contribución de ''ϕ''<sub>''i''</sub>  al elemento para cada celda, de manera que esta sea equivalente a la contribución de la función de forma ''N''<sub>''i''</sub>  del nodo <math display="inline">\boldsymbol{x}_i\in {\Omega }_e^k</math> .

​

==Observación 3.2.                     ==

​

La condición de reproducibilidad constante [[#eq0080|(11)]]  es independiente del valor del parámetro ''α'' . Es decir,

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\frac{\partial }{\partial \alpha }{\int }_{{\Omega }_e^k}{\overline{N}}_i(\boldsymbol{x};\alpha )\quad d\boldsymbol{x}=</math><math>0\quad \forall k=1,\ldots ,n_{en}(p)</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | 

|}

​

===3.2. Parámetro ''α''===

​

Una de las principales ventajas de los métodos Galerkin discontinuo es su flexibilidad respecto al espacio de aproximación. Esencialmente, cualquier espacio lineal puede usarse como espacio de aproximación local, y este puede variar elemento a elemento, así como para cada instante de tiempo. La introducción del parámetro ''α''  permite explotar esta propiedad adaptando la base localmente en función de la regularidad de la solución.      

​

Por su definición, ver [[#eq0060|(8)]] , las funciones de forma <math display="inline">{\overline{N}}_i(\boldsymbol{x};\alpha )</math>  dependen linealmente de ''α'' , que actúa como función de peso entre las bases ''N''<sub>''i''</sub>  y ''ϕ''<sub>''i''</sub> . El espacio de aproximación queda definido localmente según el valor de ''α'' . Tal y como se ha mencionado en la sección anterior, para ''α''  = 1 se tiene la base de funciones estándar, continua en el elemento, mientras que para ''α''  = 0 se introducen saltos en el interior de los elementos (en particular, en la interfaz de los subelementos <math display="inline">{\Omega }_e^k</math> ), obteniendo una base constante a trozos en Ω<sub>''e''</sub> . La magnitud de los saltos interiores está controlada por ''α''  según la expresión (1− ''α'' ) 〚 ''ϕ''<sub>''i''</sub> ('''x'''<sub>''j''</sub> )'''n'''<sub>''j''</sub>  〛 para los puntos <math display="inline">\boldsymbol{\mbox{x}}_j\in \partial {\Omega }_e^k</math> , siendo '''n'''<sub>''j''</sub>  el vector normal unitario a '''x'''<sub>''j''</sub> . El control, vía el parámetro ''α'' , de la potencia del salto también permite un control sobre la estabilización introducida por los Riemann solvers.      

​

La elección de ''α''  permite una cierta flexibilidad. En primer lugar se podría recurrir a una función ''switch''  (''α''  = 0 o 1) entre una discretización con volúmenes finitos constante a trozos o una discretización clásica de Galerkin discontinuo. Esta opción, además de ser muy restrictiva, no resulta apropiada para problemas estacionarios, ya que el uso de funciones no diferenciables dificulta alcanzar el estado estacionario. Este hecho ya ha sido señalado por otros autores haciendo referencia a técnicas de limitación o estabilización basadas en operadores discontinuos  [[#bib0145|[29]]]  and [[#bib0150|[30]]] . Otras opciones razonables consisten en modelar el parámetro ''α  ''  como una función <math display="inline">C^k</math>  con ''k''  ≥ 0. La introducción de funciones continuas permite una variación suave elemento a elemento de ''α'' , así como un rango más amplio de valores. Por simplicidad, en este artículo se ha considerado ''α  ''  como una función <math display="inline">C^0</math> , pero otras alternativas pueden ser también válidas. La [[#fig0015|figura 3]]  muestra el perfil de una posible función ''α''  y una función tipo switch.

​

<span id='fig0015'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_703683201-1-s2.0-S0213131512000491-gr3.jpg|center|359px|Función-tipo α continua y switch.]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 3.

​

Función-tipo ''α''  continua y switch.                  

​

</span>

|}

​

Con el fin de adaptar la base de funciones a la regularidad de la solución, ''α''  debe estar necesariamente controlada por un detector de choques que permita discernir gradientes pronunciados frente a perfiles suaves. En este trabajo se utiliza el sensor propuesto recientemente por Persson y Peraire  [[#bib0075|[15]]]  and [[#bib0155|[31]]] , ''S''<sub>''e''</sub> (''s'' ), definido a nivel elemental Ω<sub>''e''</sub> . Este sensor depende únicamente de la cantidad física ''s'' , obtenida a partir de la variable de estado '''U''' . Esta cantidad se conoce como ''variable de detección'' . Para el caso particular de las ecuaciones de Euler, otras referencias indican que el número de Mach es una variable de detección más precisa que otras cantidadas físicas como la entropía o la densidad, propuestas en otros trabajos  [[#bib0065|[13]]]  and [[#bib0075|[15]]] .      

​

En la [[#fig0015|figura 3]]  se distinguen 3 regiones en el perfil de la función ''α'' (''s'' ): una primera región tal que si ''S'' (''s'' ) < ''s''<sub>1</sub>  entonces ''α''  = 1, una segunda región donde ''α''  toma valores entre 0 y 1, y una tercera región ''S'' (''s'' ) > ''s''<sub>0</sub>  para la cual ''α''  = 0. La aproximación de la solución en el interior del elemento es continua para valores del sensor de la primera región (i.e., ''α''  = 1). En este sentido ''s''<sub>1</sub>  representa un umbral para discernir entre funciones suaves o funciones con gradientes pronunciados o discontinuidades.

​

==Observación 3.3.                     ==

​

La experiencia con diferentes cáculos indica que la función ''α''  se comporta de manera robusta respecto a estos 2 valores, es decir, pequeñas variaciones de ''s''<sub>0</sub>  y ''s''<sub>1</sub>  no afectan al comportamiento global de la aproximación.

​

Seguidamente se define el sensor de discontinuidades ''S'' (''s'' ), así como los valores ''s''<sub>1</sub>  y ''s''<sub>0</sub> .      

​

====3.2.1. El sensor de discontinuidades====

​

El sensor de discontinuidades se define a nivel elemental por un operador no lineal <math display="inline">S_e(s):{\Omega }_e\longrightarrow \mathbb{R}</math> , donde la variable de detección ''s''  = ''s'' ('''U''' ) es una componente o función del vector de estado '''U''' . El sensor proporciona una medida escalar de la regularidad de la aproximación.      

​

A continuación se resume brevemente el sensor para el caso unidimensional, detallado en trabajos previos de los autores [[#bib0075|[15]]]  and [[#bib0080|[16]]] . Posteriormente se extiende al caso bidimensional.      

​

Sea la base de polinomios ortonormal de Legendre, ''P''<sub>''i''</sub>  para ''i''  = 1, …, n<sub>en</sub> (''p'' ). Las aproximaciones de grado ''p''  y ''p''  − 1 de la variable de detección ''s''  en esta base se expresan como

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>s(x)=\sum_{i=1}^{n_{en}(p)}s_iP_i(x)\quad \overset{\mbox{ˆ}}{s}(x)=</math><math>\sum_{i=1}^{n_{en}(p-1)}s_iP_i(x)</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 13)

|}

​

donde, debido a las propiedades jerárquicas de la base, <math display="inline">\overset{\mbox{ˆ}}{s}(x)</math>  se ha obtenido truncando la serie polinómica ''s'' (''x'' ).      

​

El sensor de discontinuidades para el elemento unidimensional ''I''<sub>''e''</sub>  queda definido como

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>S_e(s)={log}_{10}(\frac{\parallel s-\overset{\mbox{ˆ}}{s}{\parallel }_2^2}{\parallel s{\parallel }_2^2})</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | 

|}

​

donde ∥ • ∥ <sub>2</sub>  es la norma ''L''<sub>2</sub>  de la variable. El cálculo del sensor con estas normas se simplifica gracias a la ortonormalidad de la base de polinomios de Legendre. En efecto, sea <math display="inline">\boldsymbol{\mbox{s}}^T=(s_1,\ldots ,s_{n_{en}(p)})</math>  el vector de valores nodales de la aproximación ''s'' (''x'' ), se tiene

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\parallel s{\parallel }_2^2={\int }_{I_e}s^2\quad dx=</math><math>{\int }_{I_e}(\sum_i^{n_{en}(p)}s_iP_i(x))^2\quad dx=</math><math>\boldsymbol{\mbox{s}}^T\boldsymbol{\mbox{M}}\boldsymbol{\mbox{s}}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | 

|}

​

donde '''M'''  es la matriz de masa diagonal obtenida con los polinomios ortonormales de Legendre:

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{M}}_{ij}={\int }_{I_e}P_i(x)P_j(x)\quad dx=</math><math>{\delta }_{ij}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | 

|}

​

Con esta definición y teniendo en cuenta la jerarquía de la base, el sensor de discontinuidades resulta

​

<span id='eq0115'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>S_e(s)={log}_{10}(\frac{\left(0,\ldots ,0,s_{n_{en}(p)})^T\quad \boldsymbol{\mbox{M}}\quad (0,\ldots ,0,s_{n_{en}(p)}\right)}{\boldsymbol{\mbox{s}}^T\boldsymbol{\mbox{M}}\boldsymbol{\mbox{s}}})</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 14)

|}

​

El sensor de discontinuidades en 2D se construye de manera análoga pero utilizando la base de polinomios ortonormales de Koornwinder [[#bib0160|[32]]]  para cada elemento. La expresión <math display="inline">\parallel s-\overset{\mbox{ˆ}}{s}{\parallel }_2^2</math>  en [[#eq0115|(14)]]  representa una medida de las frecuencias altas de la aproximación y vive en el espacio de polinomios <math display="inline">([P^p]\[P^{p-1}])^2</math> . En el caso bidimensional la matriz de masa se construye utilizando la base de polinomios de Koornwinder. La norma asociada a las frecuencias altas de la aproximación se obtiene de nuevo mediante una proyección de la solución ''s'' ('''''x''''' ) en el espacio de polinomios <math display="inline">([P^p]\[P^{p-1}])^2</math> . Así pues, las proyecciones en norma ''L''<sub>2</sub>  del sensor [[#eq0115|(14)]]  para el caso multidimensional se expresan como

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>E=\parallel s{\parallel }_2^2=\boldsymbol{\mbox{s}}^T\boldsymbol{\mbox{M}}\boldsymbol{\mbox{s}},\quad E_H=</math><math>\parallel s-\overset{\mbox{ˆ}}{s}{\parallel }_2^2=</math><math>\boldsymbol{\mbox{s}}^T\boldsymbol{\mbox{M}}_H\boldsymbol{\mbox{s}}\quad </math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | 

|}

​

donde '''M'''<sub>''H''</sub>  no es más que la matriz de masa calculada con los monomios de Koornwinder de grado únicamente ''p'' . Para ello se considera '''V'''  la matriz de Vandermonde en la base de Koornwinder, con las columnas ordenadas por bloques de polinomios de menor a mayor grado. Sea n<sub>L</sub>  el número de modos hasta grado ''p''  − 1 y n<sub>H</sub>  los restantes (n<sub>L</sub>  + n<sub>H</sub>  = n<sub>en</sub> (''p'' )), la matriz '''M'''<sub>''H''</sub>  es la proyección <math display="inline">\boldsymbol{\mbox{M}}_H=\boldsymbol{\mbox{F}}_H^T\boldsymbol{\mbox{M}}\boldsymbol{\mbox{F}}_H</math> , donde '''F'''<sub>''H''</sub>  es la ''matriz filtro''  siguiente

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>\boldsymbol{\mbox{F}}_H=\boldsymbol{\mbox{V}}\boldsymbol{\mbox{P}}_H\boldsymbol{\mbox{V}}^{-1}\quad con\quad \boldsymbol{\mbox{P}}_H=</math><math>diag(\underset{n_L}{\underbrace{0\ldots 0}},\underset{n_H}{\underbrace{1\ldots 1}})</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | 

|}

​

Con estas definiciones el sensor [[#eq0115|(14)]]  para cualquier dimensión es el logaritmo del cociente entre la energía asociada a las altas frecuencias y la energía total

​

<span id='eq0130'></span>

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>S_e(s)={log}_{10}(E_H/E)</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | ( 15)

|}

​

El sensor, [[#eq0115|(14)]]  y [[#eq0130|(15)]] , considera la solución de alto orden como si estuviera comprendida por una secuencia de modos de Fourier. Para una función continua los coeficientes de la serie de Fourier tienden a 0 rápidamente, a una velocidad aproximada de 1/''p''<sup>2</sup> . Teniendo en cuenta que la definición del sensor involucra cantidades al cuadrado y logaritmos, se espera que para funciones suaves <math display="inline">C^k</math> , con ''k''  ≥ 0, el sensor se comporte como un valor de orden ≲−4 log <sub>10</sub> (''p'' ). Sin embargo, en una discontinuidad todos los modos de Fourier tienen un valor significativo en la aproximación y el sensor ''S''<sub>''e''</sub>  toma valores mayores al orden esperado. Esta idea guarda una estrecha relación con los indicadores del error para métodos espectrales [[#bib0165|[33]]] .

​

==Observación 3.4.                     ==

​

El siguiente teorema relativo a la expansión en series de Fourier justifica la hipótesis relativa al comportamiento de ''S''<sub>''e''</sub> .

​

==Teorema 3.1.                     ==

​

Sea la función ''f'' (''x  '' ). Si se considera que la aproximación en términos de serie de Fourier periódica <math display="inline">SF(f)={\sum }_{n=-\infty }^{\infty }g_ne^{inx}</math> . Si <math display="inline">f(x)\in C^{k-1}</math>  y ∂<sup>''k''</sup>''f'' /∂''x''<sup>''k''</sup>  es continua, entonces |''g''<sub>''n''</sub> | ∼ ''n''<sup>−''k''</sup>  para ''n'' ⟶ ∞.

​

En base a este comportamiento la función ''α'' (''S'' ) queda unívocamente determinada por los valores ''s''<sub>1</sub>  = ''C'' (− 4 log <sub>10</sub> (''p'' )) y ''s''<sub>0</sub>  = −4 log <sub>10</sub> (''p'' ), siendo C una constante estrictamente positiva. Por experiencia numérica de los autores, se ha comprobado que el valor ''C''  = 2 proporciona resultados óptimos.      

​

La [[#fig0020|figura 4]]  muestra la tasa de decaimiento para los coeficientes de la expansión de Fourier de una función continua y una función escalón, comparadas con los límites ''s''<sub>1</sub>  y ''s''<sub>0</sub>  determinados. Se puede comprobar que en efecto, la suavidad de la función dicta la velocidad de decaimiento de los coeficientes de su serie de Fourier; veáse por ejemplo [[#bib0170|[34]]] .

​

<span id='fig0020'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_703683201-1-s2.0-S0213131512000491-gr4.jpg|center|352px|Tasa de decaimiento de los coeficientes de Fourier para curvas de diferente ...]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 4.

​

Tasa de decaimiento de los coeficientes de Fourier para curvas de diferente regularidad para el caso bidimensional. La curva log <sub>10</sub> (1/''n''<sup>8</sup> ) corresponde al umbral ''s''<sub>1</sub>  establecido.                  

​

</span>

|}

​

==Observación 3.5.                     ==

​

El sensor de discontinuidades ''S''<sub>''e''</sub> (''s'' ) es poco fiable para aproximaciones de bajo orden (''p''  ≤ 2). La hipótesis de caracterizar el comportamiento de una expansión polinómica con la serie discreta de Fourier y, por lo tanto, la validez del teorema, se verifica para grado ''p''  suficientemente alto. En consecuencia se tiene que el sensor de discontinuidades es más preciso para aproximaciones de alto orden. Los tests numéricos ratifican esta propiedad.

​

<span id='sec0035'></span>

==4. Ejemplos numéricos==

​

En esta sección se presentan diversos ejemplos numéricos de flujo compresible estacionario, en regímenes transónicos y supersónicos. Para todos los tests se han utlizado mallas triangulares y un esquema de integración temporal de tipo Runge-Kutta explícito. Para problemas estacionarios, el criterio de parada se impone mediante el residuo '''R''' ('''''U''''' ), ver [[#eq0045|(6)]] . Esta cantidad se evalúa a nivel local nodo a nodo, debido a que la solución en estado estacionario no tiene un comportamiento uniforme en todo el dominio. Los resultados obtenidos muestran el buen comportamiento de la metodología propuesta en problemas de aplicación práctica. Se demuestra el carácter local del método, así como la eficacia para altos órdenes de aproximación.      

​

===4.1. Tubo de choque o problema de Riemann===

​

El primer ejemplo es el problema clásico del tubo de choque [[#bib0175|[35]]] , ampliamente usado para validar código y metodologías en dinámica de fluidos computacional. El dominio espacial es el rectángulo de dimensiones [0, 1] × [0, 0,4]. La condición inicial utilizada es la siguiente

​

{| class="formulaSCP" style="width: 100%; text-align: center;" 

|-

| 

{| style="text-align: center; margin:auto;" 

|-

| <math>(\rho ,\boldsymbol{v},p)=\begin{array}{ll}

 & ({\rho }_l,\boldsymbol{v}_l,p_l)=(1,\boldsymbol{0},1)\quad para\quad x<1/2\\

 & ({\rho }_r,\boldsymbol{v}_r,p_r)=(0.1,\boldsymbol{0},0.125)\quad para\quad x>1/2\\

 & 

\end{array}</math>

|}

| style="width: 5px;text-align: right;white-space: nowrap;" | 

|}

​

La solución en el estado estacionario se caracteriza por una onda de choque, una discontinuidad de contacto y un abanico de expansión.      

​

El objetivo de este ejemplo es mostrar la importancia de requerir continuidad a esta función, proporcionando un rango más amplio de valores de ''α'' . Para ello se consideran 2 funciones ''α'' (''S  '' ): una función <math display="inline">C^0</math>  y una función tipo switch con umbral ''s''<sub>1</sub> ; veáse la [[#fig0025|figura 5]] (a). Todos los cálculos se han efectuado en una malla de 20 elementos longitudinales, representada en la [[#fig0025|figura 5]] (b) y con aproximaciones de grado ''p''  = 6.

​

<span id='fig0025'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_703683201-1-s2.0-S0213131512000491-gr5.jpg|center|298px|Tubo de choque: (a) funciones α(S) y (b) malla computacional.]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 5.

​

Tubo de choque: (a) funciones ''α'' (''S'' ) y (b) malla computacional.                  

​

</span>

|}

​

Los resultados obtenidos con ambas aproximaciones se muestran en la [[#fig0030|figura 6]] . Para poder hacer una mejor comparación con la solución unidimensional se ha representado un corte a lo largo de la sección horizontal para las variables densidad, velocidad y presión. Las soluciones obtenidas con la función ''α'' (''S'' ) continua proporcionan perfiles más angulosos que se ajustan mejor a la solución exacta. En cambio, la función switch introduce más difusión, llegando incluso a inhibir el perfil escalón característico de la solución (veáse por ejemplo la variable densidad en la  [[#fig0030|figura 6]] d). Además, con ''α'' -switch se obtienen perfiles que no son localmente suaves. Aunque los resultados obtenidos son satisfactorios, el método propuesto permite utilizar refinamiento adaptativo para obtenter una mejor resolución en las discontinuidades de contacto, por ejemplo, en la variable velocidad.

​

<span id='fig0030'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_703683201-1-s2.0-S0213131512000491-gr6.jpg|center|380px|Tubo de choque: perfiles para las variables densidad, velocidad y presión ...]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 6.

​

Tubo de choque: perfiles para las variables densidad, velocidad y presión calculadas con aproximaciones de grado ''p''  = 6. Arriba: con una función ''α'' (''S'' ) continua. Abajo: con una función ''α'' (''S'' ) discontinua tipo switch.                  

​

</span>

|}

​

A continuación, para confirmar el buen comportamiento del método, se ha representado en la [[#fig0035|figura 7]]  el error absoluto entre la aproximación y la solución aproximada para las variables densidad, velocidad y presión. El análisis del error confirma los resultados observados en los perfiles de las variables: los errores obtenidos con ''α'' (''S'' ) continua son ligeramente menores que los errores correspondientes a la solución con ''α'' -switch. Sin embargo, destaca en la variable velocidad el pico alrededor del punto ''x''  = 0,9, debido al suavizado en la discontinuidad de contacto. Para obtener mayor precisión se deben aplicar técnicas de refinamiento local, ya que una reducción de la cantidad de difusión (i.e., valores de ''α''  más cercanos a 1) produciría la aparición de oscilaciones en las variables densidad y presión.

​

<span id='fig0035'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_703683201-1-s2.0-S0213131512000491-gr7.jpg|center|310px|Tubo de choque: error absoluto para las funciones α(S) testeadas.]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 7.

​

Tubo de choque: error absoluto para las funciones ''α'' (''S'' ) testeadas.                  

​

</span>

|}

​

===4.2. Flujo supersónico alrededor de un bache circular===

​

El siguiente ejemplo es un test clásico en simulaciones de flujo compresible estacionario, tanto en régimen transónico como supersónico. Veánse, por ejemplo, los trabajos [[#bib0180|[36]]] , [[#bib0185|[37]]]  and [[#bib0190|[38]]] . El problema consiste en un canal de dimensiones rectangulares con un bache circular de altura 4% en la pared inferior. El bache está centrado en la mitad de la pared inferior y tiene amplitud igual a 1. El número de Mach de entrada es ''M''  = 1,4. En las paredes superior e inferior se han utilizado condicones de contorno de pared sólida no penetrante.      

​

El objetivo de este test es mostrar, por una parte, el buen comportamiento del método para mallas estructuradas y alto orden de aproximación, y, por otra parte, comprobar la robustez del método, mostrando como a medida que se refina la malla, la amplitud del choque también disminuye. Se ha resuelto el problema en 2 mallas uniformes, la primera una malla grosera de 588 elementos ([[#fig0040|fig. 8]] a) y en segundo lugar se ha utilizado una malla refinada de 1.452 elementos ([[#fig0040|fig. 8]] b). En ambos casos se ha utiliizado una aproximación de alto orden de grado ''p''  = 5.

​

<span id='fig0040'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_703683201-1-s2.0-S0213131512000491-gr8.jpg|center|340px|Flujo supersónico alrededor de un bache circular: mallas de cálculo.]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 8.

​

Flujo supersónico alrededor de un bache circular: mallas de cálculo.

​

</span>

|}

​

En las [[#fig0045|figura 9]] (a) y (b) se muestra el número de Mach obtenido con ambas aproximaciones. Cualitativamente ambas aproximaciones son comparables: el choque se captura con precisión sin propagarse en los elementos vecinos. Es remarcable, especialmente para la malla grosera, la precisión del método incluso con elementos grandes: en efecto, el choque queda contenido en un elemento sin necesidad de adaptar la malla, mostrando la capacidad del método de obtener choques de espesor del orden ''h'' /''p''  mucho menor que el tamaño del elemento. Con el fin de hacer una comparación más precisa entre ambas aproximaciones se considera una sección a lo largo del dominio correspondiente a ''y''  = 0,4. En la [[#fig0050|figura 10]]  se muestran los resultados obtenidos. Por una parte, se aprecia como en ambos casos el método permite capturar choques con alta precisión evitando oscilaciones espurias, típicamente presentes en los métodos de alto orden (''p''  ≥ 1). Por otra parte, se puede comprobar que con la malla refinada se reduce el espesor del choque y se obtienen gradientes más fuertes.

​

<span id='fig0045'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_703683201-1-s2.0-S0213131512000491-gr9.jpg|center|357px|Flujo supersónico alrededor de un bache circular: número de Mach.]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 9.

​

Flujo supersónico alrededor de un bache circular: número de Mach.

​

</span>

|}

​

<span id='fig0050'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_703683201-1-s2.0-S0213131512000491-gr10.jpg|center|356px|Flujo supersónico alrededor de un bache circular: número de Mach.]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 10.

​

Flujo supersónico alrededor de un bache circular: número de Mach.

​

</span>

|}

​

Finalmente, en las [[#fig0055|figura 11]] (a) y (b) se muestra el mapa de elementos donde se ha detectado el choque. En efecto, puesto que el sensor de discontinuidades actúa elemento a elemento, el espesor de la zona detectada para la malla fina es mucho menor que en el caso de las malla grosera.

​

<span id='fig0055'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_703683201-1-s2.0-S0213131512000491-gr11.jpg|center|342px|Flujo supersónico alrededor de un bache circular: valor del parámetro α de la ...]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 11.

​

Flujo supersónico alrededor de un bache circular: valor del parámetro ''α''  de la base de funciones. Azul (0 ≤ ''α''  ≤ 1), amarillo (''α''  = 1).                  

​

</span>

|}

​

===4.3. Supersonic NACA 0012===

​

El último ejemplo corresponde al flujo aldedor de un perfil NACA0012. El flujo de la corriente libre viene dado por un número de Mach igual a 1,2 y un ángulo de incidencia ''β''  = 0°. La malla utilizada es de 450 elementos triangulares, y 14 sobre la cuerda del perfil (en cada cara). Se ha utilizado una aproximación de grado ''p''  = 4, en contraposición a los trabajos realizados hasta el momento por otros autores  [[#bib0195|[39]]] , [[#bib0200|[40]]]  and [[#bib0205|[41]]] , donde se utilizan métodos de volúmenes finitos (''p''  = 0) o aproximaciones lineales con mallas no estructuradas y refinadas para resolver el problema. Persson y Peraire [[#bib0075|[15]]]  resuelven también este problema con alto orden y Galerkin discontinuo utilizando un método clásico de difusión artificial, pero la especificación del término de difusión supone un problema añadido que no es fácil de de solventar ([[#fig0060|fig. 12]] ).

​

<span id='fig0060'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_703683201-1-s2.0-S0213131512000491-gr12.jpg|center|357px|Flujo supersónico alrededor de un perfil NACA0012: malla de cálculo.]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 12.

​

Flujo supersónico alrededor de un perfil NACA0012: malla de cálculo.

​

</span>

|}

​

Pese a que la malla es muy grosera (compárese con la malla utilizada en otros trabajos [[#bib0210|[42]]]  de más de 10.000 elementos), el choque queda resuelto con precisión en el interior de un solo elemento, tal y como se puede apreciar en la [[#fig0065|figura 13]] (b), que muestra el detalle del número de Mach sobre la malla de fondo. Nótese que se ha utilizado una malla muy grosera (excepto en la punta del perfil aerodinámico) sin refinamiento adaptativo alrededor del choque. De nuevo en este ejemplo, la transición brusca entre la región supersónica y subsónica queda resuelta dentro de un único elemento. Cabe destacar que, en contraposición al trabajo [[#bib0075|[15]]]  donde la solución pierde la simetría en la cola del ala, pese a utilizar el mismo sensor de discontinuidades, los resultados obtenidos con el presente método mantienen la simetría.

​

<span id='fig0065'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_703683201-1-s2.0-S0213131512000491-gr13.jpg|center|357px|Flujo supersónico alrededor de un perfil NACA0012: número de Mach y detalle de ...]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 13.

​

Flujo supersónico alrededor de un perfil NACA0012: número de Mach y detalle de la malla. ''M''<sub>∞</sub>  = 1,25.                  

​

</span>

|}

​

En la [[#fig0070|figura 14]]  se representa el mapa de valores de ''α'' . El sensor permanece inactivo en prácticamente todo el dominio, en contraste con otros métodos donde se aplica estabilización no solo en la zona del choque, sino también en los elementos adyancentes al perfil  [[#bib0075|[15]]]  and [[#bib0200|[40]]] .

​

<span id='fig0070'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_703683201-1-s2.0-S0213131512000491-gr14.jpg|center|357px|Flujo supersónico alrededor de un perfil NACA0012: elementos detectados tal que ...]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 14.

​

Flujo supersónico alrededor de un perfil NACA0012: elementos detectados tal que ''α''  < 1. ''M''<sub>∞</sub>  = 1,25.                  

​

</span>

|}

​

Finalmente se muestra el coeficiente de presión, que constituye una medida común para validar la precisión del método. El perfil a ambos lados del contorno de la NACA es simétrico, y no se observan apenas oscilaciones. Los resultados además son comparables con los obtenidos en la literatura [[#bib0215|[43]]]  ([[#fig0075|fig. 15]] ).

​

<span id='fig0075'></span>

​

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 

|-

|

​

​

[[Image:draft_Content_703683201-1-s2.0-S0213131512000491-gr15.jpg|center|309px|Flujo supersónico alrededor de un perfil NACA0012: distribución de la presión en ...]]

​

​

|-

| <span style="text-align: center; font-size: 75%;">

​

Figura 15.

​

Flujo supersónico alrededor de un perfil NACA0012: distribución de la presión en la superficie superior e inferior del perfil del ala. ''M''<sub>∞</sub>  = 1,25.                  

​

</span>

|}

​

<span id='sec0055'></span>

==5. Conclusiones==

​

En este trabajo se presenta una nueva formulación para los métodos de Galerkin discontinuo que permite tratar soluciones discontinuas, no solo en la interfaz de los elementos, sino también dentro del propio elemento. El método combina la habilidad de los métodos Galerkin discontinuo de alto orden para obtener soluciones precisas utilizando mallas groseras con la tecnología estabilizadora de los métodos de volúmenes finitos. La estabilización es inherente en el método y se evita el uso de parámetros no arbitrarios. Pese a que el método permite el uso de elementos grandes y no requiere refinamiento para obtener soluciones precisas, el uso de técnicas adaptativas es compatible con el método. Con el fin de mejorar la precisión de los resultados obtenidos los autores creen que debería hacerse un estudio más preciso entre la opción de refinamiento h o p localmente.

​

==References==