[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) aparece esta frase: "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 se añade aquí. La tensión superficial implícita soltó un grillete. El mayor CFL que el artículo llegó a correr fue 0.05. El esquema que transporta la función de color había tomado el paso temporal. El diagrama NVD muestra dónde se dobla ese techo, y un vórtice pone cifra al volumen que cobra.

La fórmula capilar aparece más adelante, intacta. Aquí solo queda el $C_{\max}$ que sobrevive después de tachar esa línea.

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.

Al subir Courant C de 0.05 a 0.9, el techo de la región admisible se dobla hacia la diagonal de aguas arriba.

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.

La integral mete $C$ en los límites, a la vista. 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.

La diferencia entre los dos carriles vuelve a verse en la carrera justo después de la tabla impresa. U e interface CFL cap son los mandos de ese sitio.

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

Al bajar U del lado dominado por capilaridad, el carril A queda atado a $\Delta t_\sigma$ y el B se adelanta. Al bajar interface CFL cap a 0.01, B vuelve al lado de A. Al subirlo a 0.5, el error de volumen de abajo crece.

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 el tope en 0.05 la tabla marca 1595 pasos y un error de volumen de $2.94\times 10^{-3}$ . Al empujar el tope hacia 0.5 bajan los pasos, y el libro, 0.058 por $C$ , crece en la misma proporción.

La acotación se mantiene. Los sobrepasos por debajo quedan en $10^{-6}$ , así que THINC se quedó entre 0 y 1. El error de volumen escala con $C$ . La quinta columna, la cuarta dividida por $C$ , queda clavada cerca de 0.058. De $C = 0.05$ a $0.80$ , un factor de 16, los coeficientes leen 0.0588, 0.0587, 0.0580, 0.0583, 0.0577. La dispersión cabe en un 2%.

La ventana de integración vuelve ordinaria esa proporción. La cantidad que cruza la cara es una integral cuyo límite inferior es $1 - C$ , de modo que el ancho de la ventana es $C$ . Con $C = 0.05$ solo viaja el último 5% de la celda. Con $C = 0.8$ viaja el 80%. El número de pasos baja con ella. Con $T = 2$ y $N = 40$ , $C = 0.05$ usa 1595 pasos, 0.1 usa 797, 0.2 usa 398, 0.4 usa 199 y 0.8 usa 99. Casi una mitad cada vez. Duplicar $C$ duplica la ventana, parte los pasos a la mitad y, aun así, duplica el error de volumen por unidad de tiempo físico.

$2.94\times 10^{-3}$ con $C = 0.05$ es 0.294%. $4.62\times 10^{-2}$ con $C = 0.8$ es 4.62%. El cociente coincide con el de los números de Courant, 16. En un área bidimensional eso es cerca de 2.3% en el diámetro de la gota. La curvatura es el inverso del radio, así que el salto de presión de Laplace se inclina ese mismo 2.3%. Con $\beta = 2$ el salto de la tanh queda cerrado solo por el promedio de celda, y una ventana más ancha mantiene el sobrepaso por debajo de $10^{-6}$ . La nitidez se queda. Solo el libro de volumen sigue a $C$ .

El techo del CBC lee el mismo $C$ . Con $C = 0.05$ el techo es $20\,\tilde{\phi}_D$ : basta que el donante pase apenas de 0.05 para subir el valor de cara hasta 1. Con $C = 0.8$ el techo es $1.25\,\tilde{\phi}_D$ . 20 dividido por 1.25 da 16, el mismo factor del error de volumen. El margen de compresión baja como $1/C$ mientras el error del libro sube como $C$ .

El error de forma de la última columna recorre 0.199, 0.212, 0.185, 0.177, 0.221, con independencia de $C$ . Ese suelo lo fija la malla $40^2$ .

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

Si esa factura se apila con un error de curvatura, el libro se equivoca dos veces. Un volumen desviado un 0.5% mueve el radio en esa vecindad, y la curvatura, inverso del radio, se inclina en la misma fracción. Una tensión superficial sin equilibrio ensucia de nuevo el campo de velocidad, y la velocidad sucia sacude el $C$ del barrido siguiente. El hueco que abre la tensión superficial implícita, donde agrandar $\Delta t$ igual cierra la ventana sigue cerrado si el transporte de interfaz deja $C$ en 0.05.

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

El artículo agarra primero la robustez de la función de altura (height function). En una interfaz mal resuelta, una función de altura caída mata la curvatura entera. Después va el esquema de transporte. En palabras del artículo, un transporte que admita un CFL mayor es donde todavía hay margen de rendimiento. Volver implícito el transporte y salir del techo $1/C$ del CBC, pasar al transporte no separado del VOF geométrico (PLIC) y borrar el término de dilatación, o elegir el valor de cara aparte del transporte como en la reconstrucción MUSCL-THINC-BVD. Los tres devuelven parte del bajo costo del VOF algebraico.

Lo que queda no es un techo nuevo sino un cuello de botella nuevo. La restricción capilar dejó el asiento libre y el CFL del transporte de interfaz se sentó. $\Delta t_\sigma$ baja como $\Delta x^{3/2}$ . Este cuello baja solo como $\Delta x$ . Una malla más fina le deja relativamente más aire. Aceptar diez veces el error de volumen de 0.294% en $C = 0.05$ , o conservar los 1595 pasos, es la decisión que sigue.