¿Cómo manejar la interfaz?#

En un cálculo compresible de dos fluidos, lo que suele atascar no es el flujo numérico. Es cómo una celda anota la superficie donde los dos fluidos se tocan.

Hay tres ramas.

1. Volume of Fluid (VOF)
2. Level Set
3. Diffuse Interface (Interfaz difusa)

Al cambiar el modo entre VOF, Level Set y Diffuse, el mismo círculo cambia de espesor y de signo.

VOF encierra el color en un corte recto de la celda. Level set dibuja un contorno cero liso. Diffuse abre una banda de unas pocas celdas.

1. Volume of Fluid (VOF)#

Principio#

Se realiza el seguimiento de la fracción de volumen $\alpha$ que ocupa un fluido específico en cada celda.

$$\frac{\partial \alpha}{\partial t} + \mathbf{u} \cdot \nabla \alpha = 0$$

Si $\alpha = 1$, la celda está llena del fluido A; si $\alpha = 0$, está llena del fluido B; y si $0 < \alpha < 1$, la interfaz pasa por el interior de la celda.

Reconstrucción de la interfaz (Interface Reconstruction)#

El núcleo de VOF es reconstruir la posición de la interfaz dentro de la celda a partir de los valores de $\alpha$:

• SLIC (Simple Line Interface Calculation): Aproxima la interfaz de forma paralela a los ejes de coordenadas.
• PLIC (Piecewise Linear Interface Calculation): Aproxima la interfaz como una línea con una inclinación arbitraria.

En PLIC, el vector normal de la interfaz $\mathbf{n}$ se estima a partir del gradiente del campo $\alpha$:

$$\mathbf{n} = \frac{\nabla \alpha}{|\nabla \alpha|}$$

Ventajas y desventajas#

• Ventajas: La masa se conserva exactamente (forma conservativa).
• Desventajas: La reconstrucción de la interfaz es compleja y costosa. Difícil de extender a 3D. El cálculo de la curvatura es inexacto.

2. Level Set Method#

Principio#

Define la interfaz como el contorno cero de una función de distancia con signo (signed distance function) $\phi$.

$$\frac{\partial \phi}{\partial t} + \mathbf{u} \cdot \nabla \phi = 0$$

• $\phi > 0$: Fluido A
• $\phi < 0$: Fluido B
• $\phi = 0$: Interfaz

Reinicialización (Re-initialization)#

Durante el proceso de advección, $\phi$ pierde la propiedad de función de distancia con signo ($|\nabla\phi| = 1$). Para restaurar esto, se resuelve una ecuación de reinicialización:

$$\frac{\partial \phi}{\partial \tau} + \text{sign}(\phi_0)(|\nabla\phi| - 1) = 0$$

Donde $\tau$ es el tiempo ficticio (pseudo-time).

Curvatura y tensión superficial#

Una gran ventaja de Level Set es que el cálculo de las magnitudes geométricas es natural:

$$\kappa = -\nabla \cdot \left(\frac{\nabla\phi}{|\nabla\phi|}\right)$$

La tensión superficial se convierte en una fuerza volumétrica mediante el modelo CSF (Continuum Surface Force):

$$\mathbf{f}_s = \sigma \kappa \delta(\phi) \nabla\phi$$

Ventajas y desventajas#

• Ventajas: Fácil cálculo de curvatura/normales. Maneja cambios de topología de forma natural.
• Desventajas: Falla en la conservación de la masa. La posición de la interfaz puede desplazarse durante el proceso de reinicialización.

CLSVOF: La combinación de ambos#

Se utiliza mucho el método Coupled Level Set and VOF (CLSVOF), que combina la conservación de masa de VOF con las ventajas geométricas de Level Set.

3. Diffuse Interface Method (Interfaz difusa)#

Cambio de filosofía#

Mientras que los dos métodos anteriores ven la interfaz como algo afilado (sharp), el método de interfaz difusa considera que la interfaz tiene intrínsecamente un espesor finito.

En este enfoque, no se requiere una ecuación de seguimiento de interfaz separada, y la fracción de volumen se incluye como parte del sistema de leyes de conservación.

5-equation Model (Allaire et al., 2002; Kapila et al., 2001)#

Es el modelo reducido más utilizado:

$$\frac{\partial (\alpha_1 \rho_1)}{\partial t} + \nabla \cdot (\alpha_1 \rho_1 \mathbf{u}) = 0$$ $$\frac{\partial (\alpha_2 \rho_2)}{\partial t} + \nabla \cdot (\alpha_2 \rho_2 \mathbf{u}) = 0$$ $$\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u} \otimes \mathbf{u} + p\mathbf{I}) = 0$$ $$\frac{\partial E}{\partial t} + \nabla \cdot ((E + p)\mathbf{u}) = 0$$ $$\frac{\partial \alpha_1}{\partial t} + \mathbf{u} \cdot \nabla \alpha_1 = 0$$

La última línea es la que separa. La advección de $\alpha_1$ es no conservativa. Ahí es donde la presión de la interfaz deja de oscilar.

¿Por qué es necesaria la forma no conservativa?#

Si se resuelve $\alpha_1$ de forma conservativa:

$$\frac{\partial \alpha_1}{\partial t} + \nabla \cdot (\alpha_1 \mathbf{u}) = \alpha_1 \nabla \cdot \mathbf{u}$$

Si $\alpha_1 \nabla \cdot \mathbf{u}$ a la derecha no es cero en la divergencia discreta, la presión oscila.

Al desarrollar la advección no conservativa $\partial\alpha/\partial t + \mathbf{u}\cdot\nabla\alpha = 0$ queda $\partial\alpha/\partial t + \nabla\cdot(\alpha\mathbf{u}) - \alpha\nabla\cdot\mathbf{u} = 0$. Pasar el último término al otro lado da la ecuación de arriba. En un contacto, presión y velocidad son uniformes, así que $\nabla\cdot\mathbf{u}$ debería ser cero. Si la divergencia discreta no entrega ese cero, el lado derecho sobrevive. En una celda de interfaz con $\alpha=0.5$, si la divergencia residual es $\delta$, la fuente es $0.5\delta$. Alfa se mueve en esa cantidad, la densidad de mezcla cambia y la ecuación de estado no regresa a la misma presión. La semilla de la oscilación es ese término.

La forma no conservativa no construye ese lado derecho. Alfa sigue a $\mathbf{u}$, y la diferencia de flujos en las caras no incrementa alfa por sí sola. Por eso la última línea del modelo de cinco ecuaciones es no conservativa. La masa la guardan las dos líneas de arriba, $\partial(\alpha_k\rho_k)/\partial t + \nabla\cdot(\alpha_k\rho_k\mathbf{u})=0$. No hace falta escribir alfa en forma conservativa para que la masa de cada fase siga en su sitio.

La reinicialización ya trae su número en la fórmula. Cuando $|\nabla\phi|=1$ se rompe, $\mathrm{sign}(\phi_0)(|\nabla\phi|-1)$ empuja $\phi$ en el tiempo ficticio $\tau$. Si $|\nabla\phi|$ supera 1, el término recorta $\phi$. Si queda por debajo de 1, empuja al revés. El único alto es $|\nabla\phi|=1$. Si el cero $\phi=0$ resbala aunque sea una celda, la masa se escapa. La integral de alfa en VOF no tiene ese resbalón. A cambio, la normal $\mathbf{n}=\nabla\alpha/|\nabla\alpha|$ solo tiene pendiente en las celdas donde alfa está entre 0 y 1, así que la curvatura monta el ruido de esas celdas. $\kappa=-\nabla\cdot(\nabla\phi/|\nabla\phi|)$ devuelve la misma normal de forma lisa cuando $\phi$ es una distancia. CLSVOF pega las dos. Masa desde $\alpha$, curvatura desde $\phi$.

7-equation Model (tipo Baer-Nunziato)#

Cada fase puede llevar su propia presión y su propia velocidad. Ese camino ya está abierto en el modelo de siete ecuaciones de Baer–Nunziato.

$$\frac{\partial \alpha_k \rho_k}{\partial t} + \nabla \cdot (\alpha_k \rho_k \mathbf{u}_k) = 0, \quad k = 1, 2$$

Se converge al equilibrio de presión/velocidad mediante un proceso de relajación. Matemáticamente es el más completo, pero el coste computacional es elevado.

Resumen comparativo#

Una alfa conservativa hace oscilar la presión en el contacto#

En compresible, las filas que se quedan son las dos de abajo. El modelo de cinco ecuaciones guarda la masa en $\partial(\alpha_k \rho_k)/\partial t + \nabla\cdot(\alpha_k \rho_k \mathbf{u})=0$ y deja solo la advección de alfa en forma no conservativa. Si la divergencia discreta no anula $\alpha\nabla\cdot\mathbf{u}$, una celda de interfaz con $\alpha=0.5$ empuja la densidad de mezcla y la presión se mueve.

Volver a afilar la banda en cada paso es otra decisión. Esa operación sigue en el afilado de interfaz por antidifusión.