Una gota inmóvil se derritió a la mitad — el término de curvatura que borra el Allen

El término de afilado del Allen–Cahn conservativo no es un dispositivo para sostener el espesor de interfaz: es el término que cancela exactamente la desaparición impulsada por curvatura.

El campo de velocidad era cero y la gota igual se redujo a la mitad#

¿Debería moverse por sí sola una interfaz que tiene curvatura pero ningún flujo? Esa pregunta se salta con facilidad cuando un modelo phase-field (que trata la interfaz como una capa de transición con espesor finito) entra por primera vez en un código. Al avanzar 4000 pasos una gota 2D inicializada en reposo con la ecuación de Allen–Cahn, aparece la respuesta. El radio baja de 20.083 a 14.286 y el 49.4% de la fracción volumétrica total desaparece. El campo de velocidad se mantiene exactamente en cero todo el tiempo.

Esta entrada localiza esa pérdida término por término. Al poner el Allen–Cahn clásico al lado de la forma conservativa, la diferencia es un único término de curvatura. Después se sigue ese mismo término hasta un código de lattice Boltzmann (LB), hasta el momento en el que debe sentarse. Al terminar, la W del código y el término fuente quedan como dos promesas distintas, no una.

Comparación término por término de las dos formas#

Aspecto	Allen–Cahn clásico	Allen–Cahn conservativo
Forma	reacción–difusión	divergencia de flujo
Mantenimiento de interfaz	doble pozo $\phi(1-\phi)(1-2\phi)$	contraflujo $\frac{4}{W}\phi(1-\phi)\hat{\mathbf n}$
$\int\phi\,dV$	no se conserva	se conserva a nivel discreto
Equilibrio 1D	tanh	el mismo tanh
Respuesta a la curvatura	$v_n=-M\kappa$	cero

Las dos ecuaciones se escriben así.

\frac{\partial \phi}{\partial t} = M\left(\nabla^2\phi - \frac{16}{W^2}\phi(1-\phi)(1-2\phi)\right)

\frac{\partial \phi}{\partial t} + \nabla\cdot(\phi\mathbf u) = \nabla\cdot\left[M\left(\nabla\phi - \frac{4\phi(1-\phi)}{W}\hat{\mathbf n}\right)\right]

Aquí $\phi$ es la fracción volumétrica, $M$ la movilidad (el coeficiente que fija la rapidez de relajación de la interfaz), $W$ el espesor de interfaz, $\hat{\mathbf n}=\nabla\phi/|\nabla\phi|$ la normal a la interfaz y $\kappa=\nabla\cdot\hat{\mathbf n}$ la curvatura. Los controles M y W mueven por separado el radio del semicírculo izquierdo (clásico) y el del derecho (conservativo).

R0 = 20 luW = 4.0 luM = 0.050

t = 0  | classical R = 20.000, m/m0 = 1.0000 | conservative R = 20.000, m/m0 = 1.000000000

El semicírculo izquierdo es la forma clásica y el derecho la conservativa. Al subir M solo la mitad izquierda se derrite más rápido, siguiendo la curva punteada $\sqrt{R_0^2-2Mt}$ . Al cambiar W, la velocidad de desaparición de la izquierda casi no se mueve: el culpable es la curvatura, no el espesor.

Una identidad que la interfaz tanh regala#

Por qué $W$ entra con formas distintas en las dos ecuaciones sale del equilibrio 1D. Para una interfaz plana, la solución es este tanh.

\phi(x) = \frac{1}{2}\left(1 + \tanh\frac{2x}{W}\right)

Al derivarlo quedan dos identidades de regalo.

\frac{d\phi}{dx} = \frac{4}{W}\phi(1-\phi), \qquad \frac{d^2\phi}{dx^2} = \frac{16}{W^2}\phi(1-\phi)(1-2\phi)

La primera dice que la pendiente se puede escribir usando solo $\phi$ . No hace falta una coordenada de interfaz: con el valor local alcanza para construir el flujo de afilado. La segunda es el origen del coeficiente de doble pozo $16/W^2$ de la forma clásica. Al sustituir ambas en cada ecuación, el lado derecho queda en cero. En 1D los dos modelos comparten el mismo equilibrio.

¿El término de curvatura sobrevive o se cancela?#

La diferencia aparece en el momento en que la interfaz se curva. Con una coordenada normal $s$ , el laplaciano se parte en dos piezas.

\nabla^2\phi = \frac{\partial^2\phi}{\partial s^2} + \kappa\frac{\partial\phi}{\partial s}

Al sustituir las identidades anteriores, en el lado derecho clásico $\partial^2\phi/\partial s^2$ y el término de doble pozo se cancelan exactamente y solo sobrevive la pieza de curvatura.

\left.\frac{\partial\phi}{\partial t}\right|_{\text{classical}} = M\kappa\frac{\partial\phi}{\partial s} = M\kappa\frac{4\phi(1-\phi)}{W}

Ese único término es $v_n=-M\kappa$ , el movimiento impulsado por curvatura. Para una gota $\kappa=-1/R$ , entonces $dR/dt=-M/R$ , que integrado da $R^2=R_0^2-2Mt$ . Con radio 20 y $M=0.05$ , la vida útil es de 4000 unidades de tiempo.

La forma conservativa construye el mismo término una segunda vez dentro de la divergencia. Al expandir la divergencia del contraflujo queda

\nabla\cdot\left[\frac{4\phi(1-\phi)}{W}\hat{\mathbf n}\right] = \frac{16}{W^2}\phi(1-\phi)(1-2\phi) + \frac{4\phi(1-\phi)}{W}\kappa

así que entrega el término de doble pozo y un término de curvatura juntos. La forma clásica solo tenía la primera pieza. Al sumarse la segunda, $M\kappa\,\partial\phi/\partial s$ queda borrado. El lugar donde una fase ausente se filtraba a la interfaz en Allen–Cahn de tres fases tenía la misma estructura. Qué término se omite decide de antemano qué cantidad conservada se rompe.

El mismo tanh, dos reglas de actualización, 4000 pasos#

Ambas reglas corren sobre una malla 1D axisimétrica. Al resolver solo la dirección radial, la curvatura entra automáticamente como $1/r$ . El laplaciano y el flujo de afilado se escriben los dos como flujos de cara (face), así que la parte de divergencia es conservativa a nivel discreto en los dos modelos. Toda la pérdida de masa de la forma clásica sale entonces únicamente del término de doble pozo.

import math
 
NR, DR, R0, W, M, DT, NSTEP = 256, 0.5, 20.0, 4.0, 0.05, 0.5, 4000
rc = [(i + 0.5) * DR for i in range(NR)]          # centros de celda
rf = [(i + 1.0) * DR for i in range(NR - 1)]      # caras interiores
 
def tanh_disk(w):
    return [0.5 * (1.0 - math.tanh(2.0 * (r - R0) / w)) for r in rc]
 
def disk_mass(f):
    return sum(r * v for r, v in zip(rc, f)) * DR
 
def disk_radius(f):
    return math.sqrt(2.0 * disk_mass(f))
 
def peak_slope(f):
    return max(abs(f[i + 1] - f[i]) / DR for i in range(NR - 1))
 
def well_drive(v, w):
    return (16.0 / w ** 2) * v * (1.0 - v) * (1.0 - 2.0 * v)
 
def flux_div(F):
    out = [0.0] * NR
    for i in range(NR):
        lo = rf[i - 1] * F[i - 1] if i > 0 else 0.0
        hi = rf[i] * F[i] if i < NR - 1 else 0.0
        out[i] = (hi - lo) / (rc[i] * DR)
    return out
 
def grad_faces(f):
    return [(f[i + 1] - f[i]) / DR for i in range(NR - 1)]
 
def sharp_faces(f, w, eps=1e-12):
    F = []
    for i in range(NR - 1):
        g = (f[i + 1] - f[i]) / DR
        pf = 0.5 * (f[i] + f[i + 1])
        F.append(g - (4.0 / w) * pf * (1.0 - pf) * g / (abs(g) + eps))
    return F
 
def step_classical(f, w, dt):
    lap = flux_div(grad_faces(f))
    return [v + dt * M * (lap[i] - well_drive(v, w)) for i, v in enumerate(f)]
 
def step_conservative(f, w, dt):
    div = flux_div(sharp_faces(f, w))
    return [v + dt * M * div[i] for i, v in enumerate(f)]
 
f0 = tanh_disk(W)
m0 = disk_mass(f0)
fk, fc = list(f0), list(f0)
print("    t   R_classic  R_conserv   mass_classic  mass_conserv")
for s in range(NSTEP + 1):
    if s % 1000 == 0:
        print("%5.0f  %9.3f  %9.3f  %13.9f  %13.9f"
              % (s * DT, disk_radius(fk), disk_radius(fc),
                 disk_mass(fk) / m0, disk_mass(fc) / m0))
    fk = step_classical(fk, W, DT)
    fc = step_conservative(fc, W, DT)
print("curvature law  sqrt(R0^2-2*M*t) = %.3f" % math.sqrt(R0 ** 2 - 2 * M * NSTEP * DT))
print("interface width  init %.3f  classic %.3f  conserv %.3f"
      % (1 / peak_slope(f0), 1 / peak_slope(fk), 1 / peak_slope(fc)))
print("mass drift  classic %+.3e  conserv %+.3e"
      % (disk_mass(fk) / m0 - 1, disk_mass(fc) / m0 - 1))

El radio se queda en 20.083 y el error de masa en 1e-14#

    t   R_classic  R_conserv   mass_classic  mass_conserv
    0     20.083     20.083    1.000000000    1.000000000
  500     18.803     20.083    0.876607772    1.000000000
 1000     17.429     20.083    0.753192944    1.000000000
 1500     15.936     20.083    0.629689171    1.000000000
 2000     14.286     20.083    0.506056483    1.000000000
curvature law  sqrt(R0^2-2*M*t) = 14.142
interface width  init 4.021  classic 4.032  conserv 4.000
mass drift  classic -4.941e-01  conserv -2.631e-14

El 14.286 clásico está 1.0% por debajo del 14.142 analítico. Esa brecha viene del espaciamiento de malla 0.5 y de la definición de radio equivalente en área. La tendencia sigue $R^2=R_0^2-2Mt$ con precisión. La forma conservativa no se mueve de 20.083 hasta la tercera decimal.

Los números de masa son más contundentes. La clásica da $-4.941\times10^{-1}$ y la conservativa $-2.631\times10^{-14}$ . Lo segundo es redondeo acumulado en doble precisión y no se reduce al refinar la malla. No es un error que el refinamiento disminuya: es un error que nunca estuvo. El espesor de interfaz se queda cerca de 4 en los dos modelos, así que el espesor no es lo que los separa.

En t=0 ambos radios son 20.083 y ambas razones de masa son 1.000000000. En t=500 solo se mueve la clásica, a 18.803 y 0.876607772. En t=1000 quedan 17.429 y 0.753192944. En t=1500, 15.936 y 0.629689171. En t=2000, 14.286 y 0.506056483. La columna conservativa se queda en 20.083 y 1.000000000. El radio analítico en ese mismo instante, $\sqrt{R_0^2-2Mt}$ , es 14.142. La diferencia 14.286 menos 14.142 es 0.144, el 1.0% de arriba. La deriva de masa $-4.941\times10^{-1}$ es 1 menos 0.506056483. Junto a $-2.631\times10^{-14}$ , las escalas impresas son 1e-1 y 1e-14. Los anchos de interfaz son 4.021 al inicio, 4.032 clásico y 4.000 conservativo. Radio 20 y $M=0.05$ dan una vida de 4000 unidades de tiempo porque $R^2=R_0^2-2Mt$ se lee $400=0.1\,t$ . El t=2000 impreso es la mitad de esa vida, y 14.142 es el radio analítico ahí.

En LB este término se sienta en la fuente, no en el equilibrio#

La manera estándar de resolver la forma conservativa en LB es llevar una segunda función de distribución $h_i$ para $\phi$ . En una malla D2Q9 la actualización se escribe así.

h_i(\mathbf x+\mathbf c_i\Delta t,\ t+\Delta t) = h_i(\mathbf x,t) - \frac{h_i-h_i^{eq}}{\tau_\phi} + F_i

h_i^{eq} = \phi\,w_i\left(1+\frac{\mathbf c_i\cdot\mathbf u}{c_s^2}\right), \qquad M = c_s^2\left(\tau_\phi-\frac{1}{2}\right)\Delta t

donde $h_i$ es la distribución que transporta $\phi$ , $\tau_\phi$ el tiempo de relajación adimensional, $w_i$ los pesos de malla y $c_s^2=1/3$ . El equilibrio contiene solo advección. El término de afilado lo suministra la fuente $F_i$ , no el equilibrio.

F_i = \Delta t\,w_i\left(1-\frac{1}{2\tau_\phi}\right)\frac{4\phi(1-\phi)}{W}\frac{\mathbf c_i\cdot\hat{\mathbf n}}{c_s^2}

La expansión de Chapman–Enskog (ordenar por orden las desviaciones pequeñas del equilibrio en la razón de espaciamiento de malla) pide dos condiciones de momento. El momento de orden cero debe cumplir $\sum_i F_i = 0$ . Como $\phi=\sum_i h_i$ , un momento de orden cero distinto de cero significa que la fuente crea o destruye fracción volumétrica. El primer momento debe reproducir el contraflujo con exactitud. Como $\sum_i w_i(\mathbf c_i\cdot\hat{\mathbf n})\mathbf c_i = c_s^2\hat{\mathbf n}$ , el $c_s^2$ de $F_i$ se cancela y sale el valor buscado. Es la misma especie de restricción que por qué el equilibrio D2Q9 es un polinomio, solo que aquí la lista de líneas por respetar termina en dos.

Si se toca el momento de orden cero, la masa se fuga#

Un accidente frecuente es escribir $F_i$ como si fuera un equilibrio: meter un 1 dentro del paréntesis, como $w_i\,\mathrm{amp}\,(1+\mathbf c_i\cdot\hat{\mathbf n}/c_s^2)$ . El primer momento sigue saliendo bien. Por eso el espesor de interfaz se mantiene sano y el código corre en silencio.

CX = [0, 1, -1, 0, 0, 1, -1, -1, 1]
CY = [0, 0, 0, 1, -1, 1, 1, -1, -1]
WT = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]
CS2 = 1.0 / 3.0
W, TAU, PHI, NX, NY = 4.0, 0.8, 0.5, 0.6, 0.8
 
def sharp_source(phi, nx, ny, w, tau, spill=False):
    amp = (1.0 - 0.5 / tau) * 4.0 * phi * (1.0 - phi) / w
    out = []
    for k in range(9):
        proj = (CX[k] * nx + CY[k] * ny) / CS2        # proyeccion para el primer momento
        out.append(WT[k] * amp * (proj + 1.0 if spill else proj))
    return out
 
amp = (1.0 - 0.5 / TAU) * 4.0 * PHI * (1.0 - PHI) / W
print("target            sum c_i F_i = [%+.6f %+.6f]" % (amp * NX, amp * NY))
for tag, spill in (("moment-safe", False), ("0th-moment spill", True)):
    F = sharp_source(PHI, NX, NY, W, TAU, spill)
    sx = sum(CX[k] * F[k] for k in range(9))
    sy = sum(CY[k] * F[k] for k in range(9))
    print("%-17s sum F_i = %+.3e   sum c_i F_i = [%+.6f %+.6f]" % (tag, sum(F), sx, sy))

target            sum c_i F_i = [+0.056250 +0.075000]
moment-safe       sum F_i = +8.674e-19   sum c_i F_i = [+0.056250 +0.075000]
0th-moment spill  sum F_i = +9.375e-02   sum c_i F_i = [+0.056250 +0.075000]

W_init = 3.0 luW_model = 6.0 lu

t = 0  |  W_eff = 3.000 lu  |  max |net flux| = 0.00e+0  | mass / mass0 = 1.000000

Al activar source, el momento de orden cero se rompe. Con valores distintos en W_init y W_model, la curva roja de net flux sube y luego cae hacia el tanh punteado. W_eff sigue pegado a W_model mientras la curva de masa se aleja de 1. Un número de espesor solo no muestra ese hueco.

Los primeros momentos de las dos variantes coinciden hasta la sexta decimal. Solo se mueve el momento de orden cero, de cero a $9.375\times10^{-2}$ . Ese valor es exactamente $(1-1/2\tau_\phi)\cdot4\phi(1-\phi)/W$ , lo que significa que se está inyectando fracción volumétrica en cada nodo de interfaz y en cada paso.

El primer momento impreso es +0.056250 y +0.075000 en los dos modos. Con $\phi=0.5$ , $\tau_\phi=0.8$ y $W=4$ , $(1-1/(2\tau_\phi))\cdot 4\phi(1-\phi)/W = 0.09375$ . Por $n_x=0.6$ sale 0.056250. Por $n_y=0.8$ sale 0.075000. El momento de orden cero se parte en +8.674e-19 y +9.375e-02. 9.375e-02 es 0.09375 escrito con exponente. Un ancho de 4.000 puede cuadrar y, si esa suma no es cero, la masa se fuga.

Una tabla para elegir entre los dos modelos#

Situación	Clásico	Conservativo
Relajación de interfaz plana	apto	apto
Sostener una gota estática en corridas largas	no apto	apto
Modelar crecimiento de grano o movimiento por curvatura	apto	no apto
Reportar la fracción volumétrica total	requiere corrección	tal cual
Costo de implementación en LB	sin término fuente	requiere $\hat{\mathbf n}$

La tercera fila es la que importa. La forma conservativa no siempre es la correcta. Si el movimiento impulsado por curvatura es justo lo que se quiere modelar, la forma clásica es la física correcta. Ese término solo se vuelve un error cuando el objetivo es transportar una interfaz de fluido. En cualquiera de los dos casos hay que calcular $\hat{\mathbf n}$ , así que la forma conservativa deja pendiente proteger una división donde $|\nabla\phi|$ se acerca a cero en el volumen.

Lo que se firma al escribir W en un capturador de interfaz#

Escribir W = 4 en un código fija más que un espesor. El coeficiente de doble pozo $16/W^2$ y el coeficiente de contraflujo $4/W$ descienden los dos de ese mismo número. Si se cambia solo uno, el equilibrio deja de ser un tanh. Que ambos coeficientes se deriven de la misma $W$ es lo primero que hay que revisar.

La revisión siguiente es la cantidad conservada. Al imprimir la fracción volumétrica total en cada paso, 1e-14 se separa de 1e-1 a primera vista. El espesor se ve sano en los dos modelos, así que vigilar el espesor no filtra nada. Una forma conservativa cuyo espesor no se movió al subir la movilidad por 10 ya separaba espesor y masa. En una implementación LB la revisión cabe en una línea: imprimir $\sum_i F_i$ y confirmar que es cero.

Una gota inmóvil se derritió a la mitad — el término de curvatura que borra el Allen–Cahn conservativo