La movilidad subió 10× y el espesor de la interfaz no se movió — el término antidifusivo de Allen–Cahn conservativo

La difusión se multiplicó por diez y la interfaz no engordó#

En un código de campo de fase (phase field: la interfaz escrita como un escalar suave que pasa de 0 a 1) la movilidad $M$ pasó de 0.02 a 0.20. Eso es multiplicar por diez el coeficiente de difusión. La interfaz debería emborronarse diez veces más rápido. El espesor medido pasó de 4.44 a 4.24: bajó un poco, en lugar de subir.

No había nada roto. La ecuación está construida así. Allen–Cahn conservativo, la forma con la que Boltzmann en red (LBM) sostiene una interfaz, ata el espesor a $W$ y deja la antidifusión en una línea del equilibrio. Cómo se calcula $\nabla\phi$ sobre la red queda en los mismos números. Cada cifra de abajo salió de ejecutar el código impreso.

Al pulsar smear interface el perfil se espesa, y solo el lado con antidifusión vuelve al espesor de la línea verde discontinua.

Al pulsar smear interface el perfil se estropea a propósito; conviene entonces mirar la traza naranja de espesor de abajo. Con la antidifusión activada, cualquier $M$ devuelve el perfil a la altura de la línea verde discontinua; con ella desactivada, $M$ solo fija la velocidad a la que la interfaz sigue hinchándose.

Sujetar una interfaz con una sola ecuación#

VOF y level set resuelven primero la advección y después reconstruyen la interfaz o reinicializan la función distancia. Es una etapa aparte. Si a este mismo Allen–Cahn conservativo se le quita el término de curvatura, una gota en reposo se disuelve. Ese fallo está en el término de curvatura que borra Allen–Cahn conservativo.

Allen–Cahn conservativo pliega ese posproceso dentro de la misma ecuación.

Aquí $\phi$ es la variable de campo de fase (0 en una fase, 1 en la otra), $M$ es la movilidad, $W$ es el espesor objetivo de la interfaz en celdas de red y $\mathbf{n}=\nabla\phi/|\nabla\phi|$ es la normal a la interfaz. El primer término del corchete es la difusión ordinaria que extiende la interfaz; el segundo empuja en sentido contrario.

Que todo el miembro derecho quede dentro de una divergencia es lo importante. Con fronteras periódicas, $\int\phi\,dV$ se conserva hasta la precisión de máquina. Por eso una forma de Allen–Cahn —una ecuación célebre por no conservar la masa— puede llamarse conservativa.

El punto donde el término antidifusivo vale exactamente cero#

Sea $\xi$ la coordenada a lo largo de la normal a la interfaz. El perfil de equilibrio tiene esta forma.

Derivando se obtiene la pendiente.

Por otro lado, sobre ese mismo perfil $4\phi(1-\phi) = 1-\tanh^2(2\xi/W) = \mathrm{sech}^2(2\xi/W)$. Las dos expresiones son la misma función. Es decir, sobre el perfil de equilibrio

se cumple de forma idéntica. El corchete es cero. Y $M$ multiplica a ese cero.

Aquí se separan dos cosas. El espesor de equilibrio se cierra solo con $W$, y $M$ no entra en ese sitio. Lo que $M$ controla es el tiempo de recuperación cuando el perfil se aparta del equilibrio, y ese tiempo escala aproximadamente como $W^2/M$.

Sobre la red ese cero no se mantiene solo#

Dos términos que se cancelan de forma exacta en el continuo se separan al discretizar. Si se evalúan difusión y antidifusión por separado en los nodos con diferencias centradas, ambos esténciles ven $i-1$ e $i+1$ y saltan por encima de $i$. Los puntos pares y los impares dejan de hablarse, y la antidifusión amplifica ese modo de escala de malla. Escrito así, el perfil con $M=0.20$ se pasó del objetivo 4.0 y se quedó en 2.88.

La solución es construir el flujo en las caras y no en los nodos. Si el flujo en la cara $i+1/2$ se arma solo con $\phi_i$ y $\phi_{i+1}$, los puntos vecinos sí se ven entre sí. No queda hueco para el modo desacoplado.

Ocho medidas de espesor contra una misma regla#

Una lámina se inicializó con espesor 8 y se corrió 8.000 pasos hacia un objetivo de $W=4$. Antidifusión activada y desactivada, con cuatro valores de $M$ cada una.

import numpy as np
 
L, dx, dt = 384, 1.0, 0.2
W, u, NSTEP, R = 4.0, 0.05, 8000, 96.0
x = np.arange(L) * dx
 
def slab(width):
    d = (x - L / 2 + L / 2) % L - L / 2
    return 0.5 * (1.0 + np.tanh(2.0 * (R - np.abs(d)) / width))
 
def measured_width(phi):
    g = (np.roll(phi, -1) - np.roll(phi, 1)) / (2.0 * dx)
    return 1.0 / np.max(np.abs(g))
 
def face_flux(phi, M, sharp):
    # construye difusion y antidifusion juntas sobre la cara i+1/2
    g = (np.roll(phi, -1) - phi) / dx
    pf = 0.5 * (phi + np.roll(phi, -1))
    return M * (g - sharp * (4.0 * pf * (1.0 - pf) / W) * np.sign(g))
 
def step(phi, M, sharp):
    F = face_flux(phi, M, sharp)
    adv = u * (np.roll(phi, -1) - np.roll(phi, 1)) / (2.0 * dx)
    return phi + dt * (-adv + (F - np.roll(F, 1)) / dx)
 
def run(M, sharp, w_init=8.0):
    phi = slab(w_init)
    m0 = phi.sum() * dx
    for _ in range(NSTEP):
        phi = step(phi, M, sharp)
    return measured_width(phi), phi.sum() * dx - m0
 
print("ruler: an exact tanh of W=4.0 on this grid reads %.4f" % measured_width(slab(W)))
print("start W=8.0, %d steps at u=%.2f, target W=%.1f\n" % (NSTEP, u, W))
print("%-8s %-12s %-12s %-12s" % ("M", "anti-diff", "W_read", "mass drift"))
for M in (0.02, 0.05, 0.10, 0.20):
    for sharp in (1.0, 0.0):
        w, dm = run(M, sharp)
        print("%-8.2f %-12s %-12.4f %-+12.2e" % (M, "ON" if sharp else "OFF", w, dm))
w_lo, _ = run(0.05, 1.0, w_init=2.0)
print("\nfrom the thin side: W_init=2.0, M=0.05, ON -> %.4f" % w_lo)

ruler: an exact tanh of W=4.0 on this grid reads 4.3279
start W=8.0, 8000 steps at u=0.05, target W=4.0
 
M        anti-diff    W_read       mass drift  
0.02     ON           4.4403       -1.10e-11   
0.02     OFF          21.8974      -1.34e-11   
0.05     ON           4.2854       -2.95e-11   
0.05     OFF          32.9005      -6.76e-12   
0.10     ON           4.2442       -3.22e-11   
0.10     OFF          45.6971      -9.66e-13   
0.20     ON           4.2357       -3.60e-11   
0.20     OFF          64.0262      -2.27e-13   
 
from the thin side: W_init=2.0, M=0.05, ON -> 4.2854

La primera línea es la regla. El espesor se mide como $1/\max|\nabla\phi|$, y al colocar sobre la red una tanh exacta de ancho 4 y medirla del mismo modo se lee 4.3279 y no 4.0. Esa diferencia es el error de discretización de la propia regla.

Medidos contra esa regla, los cuatro valores con antidifusión activada, de 4.24 a 4.44, caen todos dentro de $\pm3\%$. La movilidad recorrió un factor diez y el espesor apenas se movió. Al desactivar la antidifusión, el mismo barrido se abre de 21.9 a 64.0. La deriva de masa se queda en torno a $10^{-11}$ en ambos lados.

Con $W=4$ y $M=0.02$, $W^2/M = 16/0.02 = 800$. Con $M=0.20$ vale 80. Ocho mil pasos por $dt=0.2$ dan un tiempo de integración de 1600. Eso cubre dos tiempos de recuperación con la movilidad pequeña y veinte con la grande. Por eso las cuatro filas ON se juntan en 4.4403, 4.2854, 4.2442 y 4.2357. La propia regla lee una tanh de ancho 4 como 4.3279. $1/4.3279$ es cerca de 0.231, por debajo de la pendiente máxima continua $1/W = 0.25$. Decir que 4.24 a 4.44 cae dentro de $\pm3\%$ se mide contra 4.3279, no contra 4.0.

La última línea muestra que esto no es un simple "no cambió nada". Partiendo de un espesor de 2, el perfil sube hasta 4.2854, coincidiendo en cuatro decimales con el valor alcanzado desde arriba. El perfil de equilibrio es un atractor que tira desde los dos lados.

Llevado a LBM, el término no es un término fuente#

Se le da a $\phi$ su propia función de distribución $h_i$. Aquí empieza un rodeo frecuente: programar la antidifusión como una fuente que se suma después de la colisión. No hace falta.

El equilibrio se elige así.

$w_i$ son los pesos de la red, $c_s^2=1/3$ es la velocidad del sonido de red al cuadrado y $\tau_\phi$ es el tiempo de relajación adimensional. El corchete es el primer momento del equilibrio, es decir, la velocidad que transporta $\phi$. La antidifusión se sienta ahí, junto a la velocidad física $\mathbf{u}$, disfrazada de velocidad.

Un desarrollo de Chapman–Enskog saca el coeficiente de difusión del tiempo de relajación.

Juntando ambas piezas, la ecuación de Allen–Cahn conservativa se recupera intacta. La relajación produce $+M\nabla^2\phi$; el primer momento extra del equilibrio produce $-\nabla\cdot(M\,4\phi(1-\phi)\mathbf{n}/W)$. Al dividir el $M=0.10$ impreso por $c_s^2=1/3$ sale $(\tau_\phi - 1/2)\delta t = 0.30$. Con $M=0.02$ ese producto es 0.06. Con $M=0.20$ es 0.60. Subir la movilidad por diez estira este producto. No toca $W$.

La ganancia en implementación es clara. Streaming y colisión quedan intactos. Todo el trabajo de captura de interfaz cabe en una línea del equilibrio. La tensión superficial del lado del momento es otra perilla. No se mezcla con la perilla del espesor.

El precio es calcular $\mathbf{n}$ en cada nodo#

Nada es gratis. El equilibrio contiene $\mathbf{n}=\nabla\phi/|\nabla\phi|$, así que hay que evaluar $\nabla\phi$ en todas partes, y la elección del esténcil decide el resultado.

La opción evidente es una diferencia centrada con los cuatro vecinos axiales. Boltzmann en red ofrece otra: ponderar todo el conjunto de velocidades.

En D2Q9 los pesos son $1/9$ sobre los ejes y $1/36$ en las diagonales. Al entrar los vecinos diagonales se cancela el error de anisotropía de cuarto orden. El error de ambos esténciles sobre la misma gota queda impreso en la salida de abajo.

Por qué conviene elegir el error medio más grande#

Se colocó una gota tanh de radio 24 y espesor $W$ sobre una red de 96×96, y el error relativo de $|\nabla\phi|$ en el anillo de la interfaz se agrupó en 24 sectores angulares.

import numpy as np
 
N, R = 96, 24.0
cx = np.array([0, 1, 0, -1, 0, 1, -1, -1, 1])
cy = np.array([0, 0, 1, 0, -1, 1, 1, -1, -1])
w = np.array([4/9] + [1/9]*4 + [1/36]*4)
 
def fields(W):
    i = np.arange(N)
    X, Y = np.meshgrid(i - N/2 + 0.5, i - N/2 + 0.5, indexing='ij')
    r = np.sqrt(X**2 + Y**2)
    phi = 0.5 * (1.0 + np.tanh(2.0 * (R - r) / W))
    exact = (1.0 / W) / np.cosh(2.0 * (R - r) / W)**2
    return X, Y, r, phi, exact
 
def grad_central(phi):
    gx = (np.roll(phi, -1, 0) - np.roll(phi, 1, 0)) / 2.0
    gy = (np.roll(phi, -1, 1) - np.roll(phi, 1, 1)) / 2.0
    return np.sqrt(gx**2 + gy**2)
 
def grad_lattice(phi):
    # incorpora tambien los vecinos diagonales con los pesos de la red
    gx = np.zeros_like(phi); gy = np.zeros_like(phi)
    for k in range(9):
        s = np.roll(np.roll(phi, -cx[k], 0), -cy[k], 1)
        gx += w[k] * cx[k] * s
        gy += w[k] * cy[k] * s
    return np.sqrt((3.0 * gx)**2 + (3.0 * gy)**2)
 
print("%-6s %-12s %-14s %-14s" % ("W", "stencil", "mean err %", "anisotropy %"))
for W in (5.0, 4.0, 3.0):
    X, Y, r, phi, exact = fields(W)
    ring = np.abs(r - R) < 0.6
    theta = np.arctan2(Y, X)[ring]
    for name, g in (("central", grad_central(phi)), ("D2Q9", grad_lattice(phi))):
        err = (g[ring] - exact[ring]) / exact[ring] * 100.0
        b = ((theta + np.pi) / (np.pi / 12)).astype(int) % 24
        binned = np.array([err[b == j].mean() for j in range(24) if np.any(b == j)])
        print("%-6.1f %-12s %-14.3f %-14.3f" % (W, name, np.abs(err).mean(),
                                                binned.max() - binned.min()))

W      stencil      mean err %     anisotropy %  
5.0    central      3.641          1.812         
5.0    D2Q9         4.675          0.075         
4.0    central      5.338          2.575         
4.0    D2Q9         6.803          0.155         
3.0    central      8.369          3.749         
3.0    D2Q9         10.539         0.402         

Al cambiar el botón del esténcil de diferencia centrada a D2Q9, la rosa de la derecha se pliega de cuatro pétalos a un círculo.

Con $W$ en 4, la anisotropía 2.575% frente a 0.155% es esa diferencia. Al bajar $W$ hacia 2 empeoran las dos, pero solo la diferencia centrada se deforma.

Con $W=4$ el error medio de D2Q9 es 6.803% frente al 5.338% de la diferencia centrada. Solo en exactitud, pierde. Sin embargo la dispersión con el ángulo cae de 2.575% a 0.155%, cerca de un factor diecisiete. $2.575/0.155 = 16.6$. La diferencia de error medio es $6.803 - 5.338 = 1.465$, cerca de 1.5 puntos. En $W=5$ la anisotropía baja de 1.812 a 0.075. En $W=3$ baja de 3.749 a 0.402. El error medio se mueve al revés: diferencia centrada 3.641, 5.338, 8.369, y D2Q9 4.675, 6.803, 10.539. Una interfaz más fina conserva la ganancia de isotropía y más que duplica el error de magnitud respecto de $W=5$.

Estos dos números difieren en naturaleza, no solo en tamaño. Un error repartido de forma uniforme por toda la interfaz equivoca la magnitud de $\nabla\phi$ en una proporción fija. Eso equivale a ajustar mal el coeficiente de tensión superficial, y una gota en reposo sigue en reposo. Un error que varía con el ángulo es otra cosa. Puntos distintos del mismo círculo sienten fuerzas de distinto tamaño. Ese desequilibrio agita la superficie de una gota que debería estar quieta: son las corrientes espurias de las que todo el mundo se queja.

Así que aquí se entregan 1,5 puntos porcentuales de error medio para comprar un factor diecisiete de isotropía. El trato conviene.

Qué se está regateando realmente al elegir $W$#

La tercera y la cuarta columna de esa tabla se mueven en sentidos opuestos. Bajar $W$ de 5 a 3 afila físicamente la interfaz, pero el error medio en $\nabla\phi$ sube de 3.6% a 8.4% y la anisotropía de la diferencia centrada de 1.8% a 3.7%. Subir $W$ reduce los errores, pero embota la interfaz hasta que las gotas pequeñas y las películas delgadas dejan de resolverse.

Esa tabla explica por qué $W=4$ a $5$ se volvió el valor habitual en la práctica. Recurrir a $M$ cuando se quiere una interfaz más fina no sirve de nada. El 4.44 y el 4.24 del inicio eran esa historia. Girar $M$ cambia solo la velocidad de recuperación, y bajarlo demasiado deja a la antidifusión sin capacidad de alcanzar a la difusión numérica que la advección no deja de generar. La perilla del espesor es $W$ y nada más. Cuando tres fases reciben tensiones distintas y la interfaz se escapa, esa fuga queda aparte en la tensión por pares del Allen–Cahn trifásico.