[Reseña] Se soltó el grillete capilar y quedó CFL 0.05 — el techo real del transporte de interfaz en VOF

El valor con que un esquema compresivo levanta la interfaz viene dividido por el número de Courant. Al agrandar el paso temporal, ese cociente es lo primero que desaparece.

En las conclusiones de Janodet et al. (2025) hay una frase así: "el algoritmo propuesto permite simular flujos gas-líquido realistas con pasos temporales mayores que la restricción capilar — siempre que se satisfagan las demás restricciones de paso temporal". El énfasis es mío. Se soltó un grillete al volver implícita la tensión superficial, y el mayor número de CFL que el artículo pudo correr en la práctica fue 0.05. El esquema que transporta la función de color había tomado el paso temporal en su lugar. Este texto muestra de dónde sale ese techo, sobre el diagrama NVD, y luego pone un número al costo con un experimento de transporte en un vórtice.

La restricción capilar en sí y la forma implícita de romperla ya se trataron en el artículo sobre tensión superficial implícita. Aquí va solo la continuación.

El valor que levanta la interfaz viene de aguas abajo#

En VOF algebraico — transportar la función de color directamente en lugar de reconstruir la interfaz — hay una sola razón por la que las interfaces engordan: el valor de cara aguas arriba siempre las difumina. Por eso los esquemas compresivos tiran del valor de cara hacia la celda de aguas abajo. Si se toma el valor de aguas abajo tal cual, el escalón queda encerrado en una sola celda.

El problema es que ese valor no garantiza acotación (boundedness). En cuanto la función de color cae por debajo de 0 o sube por encima de 1, la densidad se vuelve negativa y el cálculo termina. Hace falta una regla que diga hasta dónde se puede tirar. Esa regla toma el número de Courant como argumento, y eso es todo lo que trata este texto.

Conviene manipularlo en la simulación siguiente.

Courant C0.05

blend gamma_f0.60transition cells0slab width0

Start at C = 0.05: the amber window fills almost the whole box, the blue HYPER-C curve pins itself to phi~_f = 1, and the slab keeps a two-cell edge forever. Drag C toward 0.9 and the ceiling min(1, phi~_D/C) folds down onto the dashed upwind diagonal — the window is a sliver, the pink face states have nowhere to sit, and the slab bleeds out over a dozen cells.

Con Courant C en 0.05, la región admisible ámbar de la izquierda llena casi toda la caja y la losa de la derecha conserva un borde de dos celdas. Al llevar C hasta 0.9, el techo se dobla sobre la diagonal punteada, que es aguas arriba. Los puntos rosados se quedan sin lugar donde posarse y la losa se desangra con el tiempo. El esquema es el mismo, la malla es la misma. Solo creció el paso temporal.

Sobre la caja NVD, el número de Courant baja el techo#

El diagrama de variable normalizada (NVD) toma una cara y reescala los valores de la celda aguas arriba $U$ , la celda donante $D$ y la celda aceptora $A$ de esta forma.

\tilde{\phi}_D = \frac{\phi_D - \phi_U}{\phi_A - \phi_U}, \qquad \tilde{\phi}_f = \frac{\phi_f - \phi_U}{\phi_A - \phi_U}

$\tilde{\phi}_D$ indica dónde queda la celda donante entre aguas arriba y la aceptora; $\tilde{\phi}_f$ , dónde aterriza el valor de cara. $\tilde{\phi}_f = \tilde{\phi}_D$ es aguas arriba y $\tilde{\phi}_f = 1$ es aguas abajo.

El criterio de acotación convectiva de Leonard (CBC, Convection Boundedness Criterion) fija dónde debe estar el valor de cara bajo avance temporal explícito.

\tilde{\phi}_D \le \tilde{\phi}_f \le \min\!\left(1,\ \frac{\tilde{\phi}_D}{C}\right), \qquad 0 \le \tilde{\phi}_D \le 1

Aquí $C = u\,\Delta t/\Delta x$ es el número de Courant de esa cara. El significado del techo $\tilde{\phi}_D/C$ es directo: un paso drena un volumen $C$ de la celda donante, y si la función de color contenida en ese volumen supera lo que la celda tenía, la celda se vuelve negativa. Esa condición es $C \cdot \tilde{\phi}_f \le \tilde{\phi}_D$ , que reordenada da la desigualdad anterior.

Con $C = 0.05$ el techo es $20\,\tilde{\phi}_D$ . Basta que $\tilde{\phi}_D$ pase apenas de 0.05 para poder subir el valor de cara hasta 1. Con $C = 0.8$ el techo es $1.25\,\tilde{\phi}_D$ : una franja delgada justo encima de la diagonal. La compresión disponible es proporcional a $1/C$ .

Por qué CICSAM se queda en 0.01 y THINC/QQ en 0.05#

CICSAM mezcla dos curvas dentro de esta caja. Una recorre el techo mismo, HYPER-C; la otra es la más suave ULTIMATE-QUICKEST. El peso de mezcla $\gamma_f$ sale del ángulo entre la normal a la interfaz y el vector de cara. Interfaz perpendicular a la cara empuja hacia HYPER-C; oblicua, hacia UQ. Comprimir una interfaz oblicua genera arrugas artificiales en escalera, y eso es lo que la mezcla evita.

El botón CICSAM blend y el deslizador blend gamma_f de arriba son esa mezcla. Al bajar $\gamma_f$ la curva se despega del techo y la losa engorda de inmediato. También se ve que en una interfaz unidimensional alineada $\gamma_f = 1$ , con lo cual CICSAM colapsa de hecho sobre HYPER-C.

De ahí viene el techo práctico de CFL cercano a 0.01 en CICSAM. En una interfaz tridimensional real, $\gamma_f$ oscila entre 0 y 1, y solo la componente HYPER-C satisface el CBC por sí sola: la parte UQ debe volver a limitarse aparte. Para que la mezcla no cruce el techo, $C$ tiene que ser pequeño. Para esquivar esa restricción, el artículo usa THINC/QQ en lugar de CICSAM.

Una tanh vuelve a dibujar la interfaz dentro de la celda#

THINC (Tangent of Hyperbola for INterface Capturing) se salta la elección de un valor de cara y dibuja directamente la distribución interna de la celda. Sobre una celda reescalada a $[0,1]$ con coordenada $\tilde{x}$ se plantea

\Phi(\tilde{x}) = \frac{1}{2}\left[1 + \gamma \tanh\!\big(\beta(\tilde{x} - \tilde{x}_c)\big)\right]

donde $\beta$ es la nitidez de la interfaz (habitualmente cerca de 2), $\gamma = \pm 1$ es la orientación leída en las celdas vecinas y $\tilde{x}_c$ ubica el salto de la tanh. $\tilde{x}_c$ queda fijado al exigir que se reproduzca exactamente el promedio de celda $\bar{\phi}$ , y se resuelve en forma cerrada.

\tilde{x}_c = \frac{1}{\beta}\,\mathrm{artanh}\!\left(\frac{\cosh\beta - e^{\beta\gamma(2\bar{\phi}-1)}}{\sinh\beta}\right)

La cantidad que cruza la cara se obtiene integrando esta curva sobre la región de partida.

F_{i+1/2} = \int_{1-C}^{1} \Phi(\tilde{x})\,\mathrm{d}\tilde{x}

THINC/QQ añade encima una reconstrucción de superficie cuadrática (quadratic surface), que captura mejor las interfaces con curvatura. El uso de la misma tanh del lado de la captura de choques aparece en la reconstrucción TENO-THINC.

Lo esencial es que esta integral lleva $C$ de forma explícita. Al crecer $C$ , la ventana de integración se acerca al ancho completo de la celda y el esquema acaba transportando nada más que el promedio de celda. Como la tanh nunca usa el valor de aguas abajo, el CBC se satisface solo, pero la caída de la capacidad compresiva con $C$ es idéntica.

El presupuesto de paso temporal tiene tres renglones#

Ahora conviene mirar el presupuesto completo. Resolver las ondas capilares de forma explícita cuesta

\Delta t_\sigma = \sqrt{\frac{(\rho_A + \rho_B)\,\Delta x^3}{2\pi\sigma}}

y junto a eso está la restricción CFL del flujo, $\Delta t_{\text{CFL}} = C_{\max}\Delta x / U$ . Lo que hizo el artículo fue tachar el primer renglón del presupuesto. Lo que queda es el $C_{\max}$ que permita el esquema de transporte de interfaz.

A continuación se hace correr dos solvers hasta el mismo tiempo físico.

D / dx120U (flow speed)2.0e-2

interface CFL cap0.05breach factor S8.0

Leave U small — capillary-driven flow — and lane A is bound by dt_sigma while lane B runs away: that is the paper’s selling point. Now drag the interface CFL cap down to 0.01, the CICSAM value: lane B collapses back onto lane A even though surface tension is still implicit. Push the cap up to 0.5 instead and the volume-error readout is what pays for it.

Con U pequeño — flujo dominado por capilaridad — el carril A queda atado a $\Delta t_\sigma$ mientras el carril B se adelanta: ese es el argumento de venta del artículo. Ahora conviene bajar interface CFL cap hasta 0.01: el carril B retrocede junto al A aunque la tensión superficial siga siendo implícita. Al subirlo hasta 0.5, en cambio, la lectura de error de volumen de abajo indica el precio.

Un vórtice y el error de volumen que mide#

¿Qué empeora realmente al crecer $C$ ? En el transporte con separación direccional (directional splitting), cada barrido no ve el campo de velocidad solenoidal, de modo que hace falta un término de corrección por dilatación.

\phi^{*}_{i} = \phi^{n}_{i} - \left(F_{i+1/2} - F_{i-1/2}\right) + \phi^{n}_{i}\left(C_{i+1/2} - C_{i-1/2}\right)

Ese término impide que una región uniforme con $\phi = 1$ se rompa en un solo barrido. Pero en las celdas de interfaz, $\phi^n_i$ difiere del valor realmente en juego a mitad del barrido, y esa diferencia sobrevive como error de volumen. La medición se hizo montando transporte THINC sobre el vórtice único de Rider–Kothe con inversión temporal, $T=2$ .

from math import atanh, cos, cosh, exp, log, log1p, pi, sin, sinh
 
N, BETA, EPS = 40, 2.0, 1e-6
H = 1.0 / N
 
 
def lncosh(z):
    a = abs(z)
    return a + log1p(exp(-2.0 * a)) - log(2.0)
 
 
def thinc_slab(pbar, g, a, b):
    """Integra la reconstruccion tanh de la celda donante sobre [a, b]."""
    s = g * (2.0 * pbar - 1.0)
    r = max(-0.999999, min(0.999999, (cosh(BETA) - exp(BETA * s)) / sinh(BETA)))
    xc = atanh(r) / BETA
    return 0.5 * ((b - a) + (g / BETA) * (lncosh(BETA * (b - xc)) - lncosh(BETA * (a - xc))))
 
 
def face_flux(pm, p0, pp, c):
    """Fraccion de color que cruza la cara. p0 es la celda donante, c su numero de Courant."""
    if abs(c) < 1e-14:
        return 0.0
    g = 1.0 if pp > pm else (-1.0 if pp < pm else 0.0)
    if g == 0.0 or p0 < EPS or p0 > 1.0 - EPS:
        return c * p0
    return thinc_slab(p0, g, 1.0 - c, 1.0) if c > 0 else -thinc_slab(p0, g, 0.0, -c)
 
 
def line(col, vel, k):
    """Un barrido 1-D periodico, con el termino de correccion por dilatacion."""
    n, out = len(col), [0.0] * len(col)
    for i in range(n):
        cw, ce = vel[i] * k, vel[i + 1] * k
        fw = face_flux(col[(i - 2) % n], col[(i - 1) % n], col[i], cw) if cw > 0 else \
            face_flux(col[(i - 1) % n], col[i], col[(i + 1) % n], cw)
        fe = face_flux(col[(i - 1) % n], col[i], col[(i + 1) % n], ce) if ce > 0 else \
            face_flux(col[i], col[(i + 1) % n], col[(i + 2) % n], ce)
        out[i] = col[i] - (fe - fw) + col[i] * (ce - cw)
    return out
 
 
def run(courant, tend=2.0):
    uf = [[-sin(pi * i * H) ** 2 * sin(2 * pi * (j + .5) * H) for i in range(N + 1)] for j in range(N)]
    vf = [[sin(pi * j * H) ** 2 * sin(2 * pi * (i + .5) * H) for i in range(N)] for j in range(N + 1)]
    nstep = max(1, int(tend * max(abs(x) for r in uf for x in r) / (courant * H)))
    dt = tend / nstep
    f = [[1.0 if ((i + .5) * H - .5) ** 2 + ((j + .5) * H - .75) ** 2 < .15 ** 2 else 0.0
          for i in range(N)] for j in range(N)]
    f0, m0 = [r[:] for r in f], sum(sum(r) for r in f)
    for n in range(nstep):
        w = cos(pi * (n + .5) * dt / tend)          # inversion temporal de Rider-Kothe
        for ax in ((0, 1) if n % 2 == 0 else (1, 0)):
            if ax == 0:
                f = [line(f[j], [x * w for x in uf[j]], dt / H) for j in range(N)]
            else:
                cols = [[f[j][i] for j in range(N)] for i in range(N)]
                vv = [[vf[j][i] * w for j in range(N + 1)] for i in range(N)]
                cols = [line(cols[i], vv[i], dt / H) for i in range(N)]
                f = [[cols[i][j] for i in range(N)] for j in range(N)]
    lo = min(min(r) for r in f)
    hi = max(max(r) for r in f)
    dm = (sum(sum(r) for r in f) - m0) / m0
    err = sum(abs(f[j][i] - f0[j][i]) for j in range(N) for i in range(N)) / m0
    return nstep, lo, hi, dm, err
 
 
print("   C   steps    min(f)     max(f)-1     dM/M    (dM/M)/C   shape err")
for c in (0.05, 0.1, 0.2, 0.4, 0.8):
    ns, lo, hi, dm, err = run(c)
    print(f"{c:5.2f} {ns:6d}  {lo:9.2e}  {hi - 1.0:9.2e}  {dm:8.2e}  {dm / c:8.4f}   {err:8.3e}")

   C   steps    min(f)     max(f)-1     dM/M    (dM/M)/C   shape err
 0.05   1595   2.23e-29  -6.03e-07  2.94e-03    0.0588   1.988e-01
 0.10    797  -6.51e-07  -6.33e-07  5.87e-03    0.0587   2.119e-01
 0.20    398  -1.62e-06  -5.27e-07  1.16e-02    0.0580   1.850e-01
 0.40    199  -3.87e-06   1.38e-07  2.33e-02    0.0583   1.771e-01
 0.80     99  -3.20e-06   1.43e-06  4.62e-02    0.0577   2.214e-01

Hay dos cosas para leer. Primero, la acotación está intacta: los sobrepasos por debajo quedan en el nivel $10^{-6}$ , así que THINC cumplió su promesa. Segundo, el error de volumen es exactamente proporcional a $C$ . La quinta columna, que es la cuarta dividida por $C$ , queda clavada cerca de 0.058. Se mantiene dentro de un 2% mientras $C$ crece por un factor de 16.

Un error de volumen de 0.3% con $C = 0.05$ pasa a 4.6% con $C = 0.8$ . Como se trata de un área bidimensional, eso equivale a 2.3% en el diámetro de la gota. En un cálculo con tensión superficial esto es fatal: la curvatura es el inverso del radio, así que el salto de presión de Laplace se desvía en ese mismo 2.3%.

La última columna, el error de forma, oscila entre 0.18 y 0.22 con independencia de $C$ . Esa la fija la resolución de malla.

¿Alcanza con refinar la malla?#

Para averiguar de dónde sale el 0.058, se repitió $C = 0.2$ sobre $60^2$ en vez de $40^2$ . El coeficiente baja de 0.0580 a 0.0384. La razón 0.66 es prácticamente la razón de espaciados de malla $40/60$ . Es decir,

\frac{\Delta M}{M} \approx 2.3\,C\,\Delta x \propto U\,\Delta t

El error de volumen es de primer orden en el tiempo. Refinar la malla manteniendo $C$ baja el error en proporción a $\Delta x$ . Refinarla manteniendo $\Delta t$ hace crecer $C$ en la misma proporción y deja el error donde estaba. Agrandar el paso temporal en el transporte de interfaz no es gratis, y la factura se emite exactamente en proporción a $\Delta t$ .

Esa factura se agrava al superponerse con el problema de las corrientes parásitas. Un volumen desviado en 0.5% da una curvatura desviada, y una curvatura errónea se vuelve una fuerza de tensión superficial sin equilibrar que contamina otra vez el campo de velocidad.

Qué debería cambiar para pasar de 0.05 a 0.5#

El propio artículo señala dos cosas en sus conclusiones. La primera es la robustez de la función de altura (height function) implícita: si la función de altura falla en una interfaz mal resuelta, la curvatura se derrumba por completo. La segunda es el esquema de transporte de interfaz; en palabras del artículo, las mejoras "deberían permitir simulaciones con números de CFL mayores, con potencial para mejorar enormemente el rendimiento del algoritmo propuesto".

Se ven tres caminos. Volver implícito el transporte mismo y escapar del techo $1/C$ del CBC; pasar al transporte no separado del VOF geométrico (PLIC) y eliminar por completo el término de dilatación; o separar la reconstrucción de interfaz del transporte, como hace el afilado por antidifusión. Los tres ceden parte del bajo costo de cálculo del VOF algebraico.

En resumen, lo que este artículo entrega no es un techo nuevo sino un cuello de botella nuevo. La restricción capilar dejó libre su asiento y el CFL del transporte de interfaz se sentó en él, y este último no se encoge como $\Delta x^{3/2}$ al modo de $\Delta t_\sigma$ , sino como $\Delta x$ . Lo cual significa que se vuelve relativamente más benévolo a medida que la malla se refina. Es información útil a la hora de elegir el próximo cuello de botella.