El problema de Riemann: El corazón de CFD#

En el Método de Volúmenes Finitos (FVM), la precisión y estabilidad de la solución dependen de cómo se calcule el flujo numérico (numerical flux) en los límites de las celdas. La clave para determinar este flujo numérico es el problema de Riemann.

¿Qué es el problema de Riemann?#

Es un problema de valor inicial donde las condiciones iniciales consisten en dos estados constantes separados por una única discontinuidad:

$$\mathbf{U}(x, 0) = \begin{cases} \mathbf{U}_L & \text{if } x < 0 \\ \mathbf{U}_R & \text{if } x > 0 \end{cases}$$

Para las ecuaciones de Euler en 1D, la solución analítica de este problema consta de tres ondas (waves):

1. Onda hacia la izquierda (rarefacción o choque)
2. Discontinuidad de contacto (contact discontinuity)
3. Onda hacia la derecha (rarefacción o choque)

La idea de Godunov (1959)#

La visión fundamental de Godunov es simple pero poderosa:

“

Para obtener el flujo en el límite de una celda, se deben tomar los valores promedio de las dos celdas adyacentes como los estados izquierdo y derecho, y resolver el problema de Riemann.

Fórmula de actualización del Método de Volúmenes Finitos:

$$\mathbf{U}_i^{n+1} = \mathbf{U}_i^n - \frac{\Delta t}{\Delta x}\left(\hat{\mathbf{F}}_{i+1/2} - \hat{\mathbf{F}}_{i-1/2}\right)$$

Donde $\hat{\mathbf{F}}_{i+1/2}$ es precisamente el flujo numérico determinado a partir de la solución del problema de Riemann.

Exact Riemann Solver (Solucionador exacto)#

El solucionador exacto de Riemann encuentra la presión en la región central (star region), $p^*$, utilizando el método iterativo de Newton-Raphson. Se aplican las relaciones de Rankine-Hugoniot o las relaciones isentrópicas dependiendo del tipo de onda (choque/rarefacción).

Relación para la onda de choque (Shock wave) ($p^* > p_K$, $K = L$ o $R$):

$$f_K(p^*) = (p^* - p_K)\left[\frac{A_K}{p^* + B_K}\right]^{1/2}$$

Donde $A_K = \frac{2}{(\gamma+1)\rho_K}$, $B_K = \frac{\gamma-1}{\gamma+1}p_K$.

Relación para la onda de rarefacción (Rarefaction wave) ($p^* \leq p_K$):

$$f_K(p^*) = \frac{2c_K}{\gamma - 1}\left[\left(\frac{p^*}{p_K}\right)^{\frac{\gamma-1}{2\gamma}} - 1\right]$$

$p^*$ satisface la siguiente condición:

$$f_L(p^*) + f_R(p^*) + (u_R - u_L) = 0$$

Approximate Riemann Solver: ¿Por qué es necesario?#

Los solucionadores exactos son costosos porque requieren cálculos iterativos. En la mayoría de los códigos CFD prácticos, se utilizan solucionadores de Riemann aproximados (approximate Riemann solvers).

Roe Solver (1981)#

Roe linealiza el problema de Riemann no lineal, convirtiéndolo en un problema de Riemann lineal para la matriz $\hat{\mathbf{A}}$:

$$\hat{\mathbf{F}}_{i+1/2} = \frac{1}{2}(\mathbf{F}_L + \mathbf{F}_R) - \frac{1}{2}\sum_{k=1}^{3}|\hat{\lambda}_k|\hat{\alpha}_k\hat{\mathbf{r}}_k$$

Se construye $\hat{\mathbf{A}}$ utilizando el promedio de Roe (Roe average):

$$\hat{u} = \frac{\sqrt{\rho_L}\,u_L + \sqrt{\rho_R}\,u_R}{\sqrt{\rho_L} + \sqrt{\rho_R}}, \quad \hat{H} = \frac{\sqrt{\rho_L}\,H_L + \sqrt{\rho_R}\,H_R}{\sqrt{\rho_L} + \sqrt{\rho_R}}$$

Ventajas: Captura con precisión las discontinuidades de contacto. Desventajas: Puede violar las condiciones de entropía (requiere una corrección de entropía o "entropy fix").

HLLC Solver (Toro, 1994)#

HLLC es una mejora del solucionador HLL que restaura la discontinuidad de contacto (Contact). Estima tres velocidades de onda $S_L$, $S^*$, $S_R$:

$$\hat{\mathbf{F}}_{i+1/2} = \begin{cases} \mathbf{F}_L & \text{if } S_L > 0 \\ \mathbf{F}^*_L & \text{if } S_L \leq 0 < S^* \\ \mathbf{F}^*_R & \text{if } S^* \leq 0 \leq S_R \\ \mathbf{F}_R & \text{if } S_R < 0 \end{cases}$$

Flujo en la región central:

$$\mathbf{F}^*_K = \mathbf{F}_K + S_K(\mathbf{U}^*_K - \mathbf{U}_K)$$

HLLC es sencillo de implementar y robusto, por lo que es el más utilizado en códigos de flujo multifásico compresible.

AUSM+ (Liou, 1996)#

La familia AUSM separa el flujo en términos convectivos (convective) y de presión (pressure):

$$\hat{\mathbf{F}}_{i+1/2} = \dot{m}_{1/2}\boldsymbol{\Psi}_{1/2} + \hat{p}_{1/2}\mathbf{D}$$

Donde $\dot{m}_{1/2}$ es el flujo másico en el límite de la celda y $\boldsymbol{\Psi}$ es la cantidad convectada. Es estable incluso en flujos de bajo Mach, lo que lo hace ideal para extenderlo a esquemas de todas las velocidades (all-speed schemes).

Extensión a flujos multifásicos#

En flujos multifásicos, dado que las EOS en ambos lados de la interfaz son diferentes, los solucionadores de Riemann deben generalizarse a versiones multi-material.

Enfoques representativos:

• Ghost Fluid Method: Construye celdas fantasma (ghost cells) del fluido del lado opuesto cerca de la interfaz y resuelve problemas de Riemann de un solo material independientes en cada fase.
• HLLC para multi-material: Determina $p^*$, $u^*$ aplicando cada EOS a ambos lados de la región central.

En el próximo artículo, compararemos las técnicas de captura de interfaz (VOF, Level Set, Diffuse Interface).

Lleva la relación de presiones a 1000:1 para ver cómo ambas ondas se vuelven choques.