Método de Galerkin Discontinuo (DGM): Discretización de alta precisión para el análisis de flujos compresibles en mallas no estructuradas

Dentro de la celda, grado $K$. En la cara, un flujo numérico. En una malla no estructurada ese par es DGM.

Al mover p, el número de elementos o smooth/step cambian la curva interior y el salto en la cara.

En smooth, p=3 se pega a la curva. En step, el mismo p sigue oscilando en el borde.

Ecuación de Gobierno y Separación de Flujos#

Ecuación de conservación para flujo compresible:

$\frac{\partial \mathbf{U}}{\partial t} + \nabla \cdot \mathbf{F}^c(\mathbf{U}) - \nabla \cdot \mathbf{F}^v(\mathbf{U}, \nabla \mathbf{U}) = \mathbf{S}(\mathbf{U})$

Donde:

• $\mathbf{U} = [\rho,\ \rho u,\ \rho v,\ \rho w,\ \rho E]^T$: Vector de variables conservativas
• $\mathbf{F}^c$: Flujo convectivo (convective)
• $\mathbf{F}^v$: Flujo viscoso (viscous)
• $\mathbf{S}$: Término fuente

Expansión de Funciones de Base y Derivación de la Forma Débil#

Dentro de la celda $\Omega_e$, la solución se aproxima como una combinación lineal de funciones de base:

$\mathbf{U}_h(\mathbf{x}, t) = \sum_{l=1}^{N_p} \hat{\mathbf{U}}_l(t)\, \phi_l(\mathbf{x})$

$N_p$ es el número de grados de libertad dentro de la celda, determinado por el grado del polinomio de aproximación $K$ ($N_p = (K+1)(K+2)(K+3)/6$ en 3D).

Multiplicando por una función de prueba $\phi_k$ e integrando sobre $\Omega_e$:

$\int_{\Omega_e} \phi_k \frac{\partial \mathbf{U}_h}{\partial t}\, d\Omega + \int_{\Omega_e} \phi_k \nabla \cdot \mathbf{F}\, d\Omega = \int_{\Omega_e} \phi_k \mathbf{S}\, d\Omega$

Al aplicar la fórmula de Green (integración por partes), se divide en integral de volumen e integral de superficie (contorno):

$\hat{\mathbf{F}}$ en el contorno es el flujo numérico calculado utilizando información de las celdas adyacentes (por ejemplo, Roe, HLLC, etc.).

Matriz de Masa e Integración Temporal#

Organizando como una Ecuación Diferencial Ordinaria (ODE) para el $l$-ésimo coeficiente de la función de base:

$\mathbf{M} \frac{d\hat{\mathbf{U}}}{dt} = \mathbf{R}(\hat{\mathbf{U}})$

Donde la matriz de masa (mass matrix) es:

$M_{kl} = \int_{\Omega_e} \phi_k \phi_l\, d\Omega$

Si se utilizan funciones de base ortogonales (como la base de Taylor o los polinomios de Legendre), $\mathbf{M}$ se convierte en una matriz diagonal, eliminando el costo de calcular la matriz inversa. El vector residual $\mathbf{R}$ es la suma del término de integral de volumen y el término de flujo de superficie:

$R_k = \int_{\Omega_e} \nabla \phi_k \cdot \mathbf{F}\, d\Omega - \oint_{\partial \Omega_e} \phi_k \hat{\mathbf{F}} \cdot \mathbf{n}\, d\Gamma + \int_{\Omega_e} \phi_k \mathbf{S}\, d\Omega$

Aplicación de la Cuadratura de Gauss#

Las integrales de volumen y superficie se calculan numéricamente mediante la cuadratura de Gauss-Legendre. Para integrar con precisión un polinomio de grado $K$, se requiere una integración de al menos orden $2K-1$.

Número de puntos de Gauss por tipo de celda (por grados representativos):

Para los tipos de cara (face), se utiliza el triple de los puntos de la celda (proyección 3D).

$K=0,1,2,3$ en $N_p=(K+1)(K+2)(K+3)/6$ da 1, 4, 10 y 20. Cinco variables conservadas dan 5, 20, 50 y 100 incógnitas por celda. Los puntos 1, 4, 10 del tetraedro coinciden con $N_p$ en $K=0,1,2$, un grado bajo 4, 10, 20. Los 1, 8, 27 del hexaedro son $1^3,2^3,3^3$. La regla de triplicar asigna 3, 12 y 30 puntos de cara. El orden $2K-1$ vale 1, 3 y 5. Con menos puntos que grados de libertad, el flujo no lineal pierde ese orden primero. Quitar un punto hasta divergir está en el piso de cuadratura de la base de Taylor.

Secuencia de Implementación del Algoritmo#

En la implementación real, se procede en el siguiente orden:

1. Clasificación de celdas y caras: Tetra/Hexa/Prisma/Pirámide/Poliedro.
2. Preparación de funciones de base: Calcular $\phi_l, \nabla\phi_l$ adecuadas para cada tipo de celda × grado de aproximación.
3. Caché de puntos de Gauss: Almacenar las coordenadas y pesos del elemento de referencia (reference element) en la memoria.
4. Inicialización de variables: Establecer $\hat{\mathbf{U}}_l^0$ (condiciones de flujo libre, etc.).
5. Bucle (avance temporal):
   • Integral de volumen: $\int_{\Omega_e} \nabla\phi_k \cdot \mathbf{F}\ d\Omega$ — Calcular y sumar flujos en cada punto de Gauss.
   • Flujo de superficie: $\oint_{\partial\Omega_e} \phi_k \hat{\mathbf{F}}\cdot\mathbf{n}\ d\Gamma$ — Calcular el flujo numérico.
   • Actualización RK: Integración temporal explícita como SSP-RK3.

import numpy as np
 
def compute_volume_residual(U_hat, phi, dphi, gauss_w, flux_func):
    """Cálculo del residual de volumen (celda individual)"""
    Np, Ng = phi.shape   # (número de funciones de base, número de puntos de Gauss)
    R = np.zeros_like(U_hat)
    for g in range(Ng):
        # Reconstrucción de magnitudes físicas en el punto de Gauss
        U_g = phi[:, g] @ U_hat          # shape: (nvar,)
        F_g = flux_func(U_g)             # (nvar, ndim)
        # Contribución al residual de volumen
        for k in range(Np):
            R[k] += gauss_w[g] * (dphi[k, :, g] @ F_g.T)
    return R

Selección del Flujo Numérico — Advertencias Prácticas#

En DGM, el flujo numérico está directamente relacionado con la estabilidad y la precisión.

• Flujo de Roe: Captura bien las discontinuidades de contacto, pero puede provocar choques de expansión (expansion shock) si no hay corrección de entropía.
• HLLC: Más estable que Roe, pero en casos raros puede ocurrir una disipación excesiva.
• LxF (Lax-Friedrichs): Implementación simple, pero la disipación es grande, lo que reduce la ganancia de alto orden.

Errores comunes:

1. Puntos de Gauss insuficientes: En flujos no lineales, se producen errores de aliasing que reducen el orden.
2. Ignorar la diagonalización de la matriz de masa: Si se utilizan funciones de base de FEM comunes, se requiere resolver un sistema lineal en cada paso, lo cual es ineficiente.
3. Inconsistencia en la dirección normal de la cara: Verificar siempre que el signo de la normal de la celda adyacente sea opuesto en la cara compartida.

Verificación de Precisión — Vórtice de Taylor-Green#

Al comparar la solución analítica incompresible con los resultados de DGM, el DGM de orden $K$ debería mostrar teóricamente una convergencia de $O(h^{K+1})$. En la práctica, si el error no disminuye $2^{K+1}$ veces cuando la malla se refina al doble, conviene revisar el orden de cuadratura o el tratamiento de contorno.

Partir el espaciado a la mitad da el cociente $2^{K+1}$: 4, 8 y 16 en $K=1,2,3$. Clavado en 2, el orden observado es 1. Con puntos un grado abajo, subir $K$ tampoco da ese cociente.

¿Fiarse de K=1 cuando el factor de refinamiento no es 4?#

En $K=1$, bajar por 2 y no por 4 no es una corrida $O(h^{K+1})$. Los puntos 1, 4, 10 quedan una columna bajo los grados de libertad 4, 10, 20. Con esa columna, el estudio de orden no separa el integrador de la cuadratura.

El contraste es $N_p=(K+1)(K+2)(K+3)/6$ y $2^{K+1}$. En $K=2$, 10 grados de libertad y cociente 8; en $K=3$, 20 y 16. LxF puede comerse la ganancia aunque el cociente cuadre. Donde esa disipación aplasta un vórtice pequeño: el esquema de baja disipación.

Referencias#

1. Luo, H. et al., "A discontinuous Galerkin method based on a Taylor basis for the compressible flows on arbitrary grids," J. Comput. Phys. 227 (2008).
2. Luo, H. et al., "Hybrid discontinuous Galerkin–finite volume techniques for compressible flows on unstructured meshes," AIAA Paper (2012).