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

En el método de volúmenes finitos (FVM), la precisión y la estabilidad dependen del flujo numérico (numerical flux) en la cara de la celda. La clave es un problema de Riemann: dos estados constantes y un solo salto.

Al llevar el control de la razón de presiones desde cerca de 1 hasta 1000:1, las dos ondas exteriores pasan de rarefacción a choque.

El contacto del medio sigue ahí aunque la razón crezca. Las dos ondas exteriores son choques a la vez solo cuando la presión de la estrella $p^*$ queda por encima de la presión de cada lado.

¿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$$

El corte es una sola comparación entre $p^*$ y $p_K$. Si $p^*$ queda por encima de $p_K$, ese lado es un choque. $f_K$ es el salto de presión por $[A_K/(p^*+B_K)]^{1/2}$. Como $A_K=2/((\gamma+1)\rho_K)$, un fluido más denso cambia menos de velocidad ante el mismo salto. $B_K=(\gamma-1)p_K/(\gamma+1)$ ensancha un poco el denominador por encima de $p^*$. Si $p^*$ no supera a $p_K$, ese lado es una rarefacción. El factor $2c_K/(\gamma-1)$ crece con la velocidad del sonido y otra vez cuando $\gamma$ se acerca a 1. El exponente de la razón de presiones es $(\gamma-1)/(2\gamma)$. Si $\gamma$ baja, el exponente baja y $f_K$ depende menos de esa razón.

La suma de las dos $f$ más $u_R-u_L$ tiene que dar cero. Si $u_R$ es menor que $u_L$, el fluido se acumula en el centro. Entonces $f_L+f_R$ es positiva y $p^*$ queda por encima de ambas presiones laterales. Las dos ondas exteriores son choques. El extremo 1000:1 del control es ese límite. El $p_K$ del lado bajo queda bajo la presión de la estrella, así que esa $f_K$ se queda en la rama de choque. El lado alto puede seguir en rarefacción, según la razón. Los dos lados son choques solo cuando la diferencia de velocidad que cierra sobre el contacto es bastante grande. Tratar una celda como su propio tubo de choque es el cálculo de el tubo de choque de Godunov dentro de una celda.

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.

El artículo siguiente compara las técnicas de captura de interfaz (VOF, Level Set, Diffuse Interface).

Separar la advección de la presión en un flujo multifásico sigue en la separación TV del flujo de Baer–Nunziato.