Pedí σ₁₂ = 1.00 y recibí 0.78 — dónde se filtra una fase ausente en el Allen–Cahn de tres componentes

Pedí σ₁₂ = 1.00 y recibí 0.78#

Entraron tres fases en un código phase-field y en el campo de tensión superficial 1-2 quedó escrito 1.00. Después vino la integral de la energía libre a través de la interfaz relajada. El resultado fue 0.78.

La fase que no debía contar era la 3. No estaba en ninguna parte del dominio: solo había quedado sembrada en $10^{-3}$. Y sin embargo, en equilibrio su fracción de volumen alcanzaba 0.21 justo en el centro de la interfaz 1-2.

Este artículo persigue ese 21% hasta el fondo. Invierte las $\sigma_{ij}$ escritas por pares en coeficientes de gradiente $\kappa_i$ por componente, muestra que el signo de $\kappa_i$ es la misma desigualdad que cierra el triángulo de Neumann, y localiza el único término que mantiene una fase ausente fuera de la interfaz. Al final, cuatro cierres corren sobre la misma malla y deciden los números.

N fracciones de volumen no son N grados de libertad#

Con $N$ fluidos inmiscibles, las fracciones de volumen $c_i$ arrastran una restricción.

Quedan $N-1$ grados de libertad. Pero el código reserva $N$ arreglos, y colgar una ecuación de Allen–Cahn de cada arreglo parece lo natural. Al copiar un solucionador de dos fases y subir el número de componentes sale exactamente esa forma.

Entonces la suma se rompe. Si cada potencial químico $\mu_i$ se escribe de forma independiente y se avanza con $\partial_t c_i = -M \mu_i$, resulta

y nada obliga al lado derecho a anularse. En la corrida de abajo, $\left|\sum c_i - 1\right|$ se abre hasta 0.859, la interfaz se disuelve y la tensión superficial cae a 0.026.

La receta de libro es un multiplicador de Lagrange. Al restar $\bar\mu = \frac{1}{N}\sum_j \mu_j$,

la suma se conserva exactamente. Conviene alternar entre los dos cierres en la simulación siguiente.

Con 3 equations, no closure, la traza de $\sum c_i - 1$ se desploma hacia abajo y las tres curvas se difuminan juntas. Al pasar a mean projection, la traza se clava en $10^{-13}$, pero un bulto rojo de $c_3$ se levanta en el centro de la interfaz y sigma_eff queda por debajo de 1.00. Subir el deslizador Lambda aplana ese bulto.

Invertir σᵢⱼ por pares en κᵢ por componente#

La energía libre se escribe normalmente por pares.

$\beta_{ij}$ es la profundidad del pozo doble del par $i$-$j$ y $\kappa_i$ el coeficiente de energía de gradiente del componente $i$. Los pozos van por pares; los coeficientes de gradiente, por componente. Esa asimetría genera todo lo que viene después.

Ni $\beta_{ij}$ ni $\kappa_i$ se eligen con libertad. El límite de dos fases los amarra a cantidades físicas. Sustituyendo $c_2 = 1 - c_1$,

que coincide con la forma habitual $\beta c^2(1-c)^2 + \frac{\lambda}{2}|\nabla c|^2$ si $\lambda = \kappa_1 + \kappa_2$. Integrar una vez el equilibrio 1D $\lambda c'' = dE_0/dc$ da $\frac{\lambda}{2}c'^2 = E_0$, es decir $c' = \sqrt{2\beta/\lambda}\, c(1-c)$. La solución es una tangente hiperbólica, y

donde $D$ es el espesor de interfaz en celdas y $\sigma$ la tensión superficial. Invirtiendo ambas relaciones,

$\beta_{ij}$ sale par por par. $\kappa_i$ exige resolver las tres ecuaciones a la vez. Para tres fases,

Un coeficiente por componente es una combinación de tres tensiones por pares. Esta es la identidad central del artículo. En el artículo sobre Allen–Cahn conservativo de dos fases, con la movilidad y el término antidifusivo había un solo $\lambda$. Con tres fases ese $\lambda$ se parte en tres, y se parte por resta.

Con $D=4$ el factor $\frac{3D}{4}$ vale $3$. En la segunda fila de la tabla, $\sigma=(1.00,\,0.80,\,0.60)$ da $\kappa_1=3\times1.20=3.60$, $\kappa_2=3\times0.80=2.40$ y $\kappa_3=3\times0.40=1.20$. Los tres coeficientes suman $7.20$. El par $1$-$2$ solo suma $6.00$. La relación de dos fases $\kappa_i+\kappa_j=\frac{3}{2}\sigma_{ij}D$ con $1.00$ y $4$ también devuelve $6.00$. Al meter $\lambda=6$ en $\sigma=\frac{2\lambda}{3D}$ el cociente es $\frac{12}{12}=1.00$. La fórmula continua coincide con la entrada. El caso A en la malla devuelve $0.9665$, y ese $3.3\%$ es lo que cuestan cuatro celdas a través de la interfaz. El mismo factor $3$ sobre $\sigma=(1.00,\,0.40,\,0.40)$ da $\kappa_3=-0.60$. Esa es la fila en la que el triángulo no cierra. En $(1.00,\,0.55,\,0.55)$, $\kappa_3=0.30$ sigue positivo y las cuñas son $155.4^\circ$, $155.4^\circ$, $49.2^\circ$. El margen antes de que el signo muera es $0.30$, y el ángulo de la fase 3 ya bajó a $49^\circ$. $\kappa_3=1.20$ no entra en la cuenta de $\sigma_{12}$. Aparece como energía de gradiente solo cuando la fase 3 ocupa la interfaz. Tensiones iguales $\sigma=(1.00,\,1.00,\,1.00)$ ponen cada $\kappa$ en $3.00$ y cada cuña en $120^\circ$. La suma es $360^\circ$. La cuña de $49.2^\circ$ en la fila $0.55$ es $70.8^\circ$ más estrecha que esos $120^\circ$. El ángulo se aprieta antes de que el signo cambie.

Dónde κ₃ cruza a negativo#

Restar significa que el signo puede perderse. Si $\sigma_{13} + \sigma_{23} < \sigma_{12}$ entonces $\kappa_3 < 0$ y la energía de gradiente se abre hacia abajo. Afilar la interfaz ahora baja la energía, así que la solución colapsa hasta el espaciado de malla. El modelo no se degrada: el funcional no está acotado por debajo desde el principio.

La misma condición se conoce con otro nombre. En una unión triple las tres tensiones se equilibran y fijan el ángulo de cuña $\theta_i$ de cada fase.

$\theta_1$ es el ángulo que ocupa la fase 1; para los otros dos basta rotar los índices. La condición para que este coseno salga de $[-1, 1]$ es exactamente la condición $\kappa_i < 0$. Es decir, que el triángulo de Neumann cierre y que el modelo esté bien planteado son la misma desigualdad. Con $\sigma_{12}=1.00$ fijo, cuando $\sigma_{13}+\sigma_{23}$ cae por debajo de $\sigma_{12}$ el triángulo se rompe y la cuña de la fase 3 se abre a 180°. Eso es la fase 3 mojando por completo la interfaz 1-2. En $\sigma=(1.00, 0.80, 0.60)$ las cuñas salen 143.1° / 126.9° / 90.0°. Un triángulo rectángulo 3-4-5, así que el ángulo de la fase 3 es exactamente recto.

Cuatro cierres sobre la misma malla#

Malla 1D de 200 celdas, $D=4$, $\sigma=(1.00, 0.80, 0.60)$. Fase 1 a la izquierda, fase 2 a la derecha, fase 3 sembrada en $10^{-3}$ solo en la interfaz. Tras 8000 pasos de relajación se miden $\left|\sum c_i - 1\right|$, $\max c_3$ y la integral de energía libre a través de la interfaz, que es la tensión superficial realizada.

import math
import numpy as np
 
D, MOB, DT, NX = 4.0, 1.0 / 3.0, 0.2, 200
 
 
def pair_to_component(s12, s13, s23):
    """Tensiones superficiales por pares -> coeficientes de gradiente kappa_i y pozos beta_ij."""
    f = 0.75 * D
    kap = np.array([f * (s12 + s13 - s23),
                    f * (s12 + s23 - s13),
                    f * (s13 + s23 - s12)])
    bet = 12.0 / D * np.array([[0, s12, s13], [s12, 0, s23], [s13, s23, 0]])
    return kap, bet
 
 
def neumann_wedges(s12, s13, s23):
    """Angulos de cuna de las tres fases en grados si el triangulo de Neumann cierra, si no None."""
    args = [(s23 ** 2 - s12 ** 2 - s13 ** 2) / (2 * s12 * s13),
            (s13 ** 2 - s12 ** 2 - s23 ** 2) / (2 * s12 * s23),
            (s12 ** 2 - s13 ** 2 - s23 ** 2) / (2 * s13 * s23)]
    if any(abs(a) > 1 for a in args):
        return None
    return [math.degrees(math.acos(a)) for a in args]
 
 
def _ext(a):
    return np.concatenate(([a[0]], a, [a[-1]]))
 
 
def lap1d(a):
    e = _ext(a)
    return e[:-2] - 2 * a + e[2:]
 
 
def grad1d(a):
    e = _ext(a)
    return (e[2:] - e[:-2]) / 2.0
 
 
def chem_potential(c, kap, bet, lam3):
    c1, c2, c3 = c
    dbulk = np.array([
        2 * bet[0, 1] * c1 * c2 ** 2 + 2 * bet[0, 2] * c1 * c3 ** 2 + 2 * lam3 * c1 * (c2 * c3) ** 2,
        2 * bet[0, 1] * c2 * c1 ** 2 + 2 * bet[1, 2] * c2 * c3 ** 2 + 2 * lam3 * c2 * (c1 * c3) ** 2,
        2 * bet[0, 2] * c3 * c1 ** 2 + 2 * bet[1, 2] * c3 * c2 ** 2 + 2 * lam3 * c3 * (c1 * c2) ** 2,
    ])
    return np.array([dbulk[i] - kap[i] * lap1d(c[i]) for i in range(3)])
 
 
def line_tension(c, kap, bet, lam3):
    """Integral de energia libre a traves de la interfaz = tension superficial realizada."""
    c1, c2, c3 = c
    e = (bet[0, 1] * (c1 * c2) ** 2 + bet[0, 2] * (c1 * c3) ** 2
         + bet[1, 2] * (c2 * c3) ** 2 + lam3 * (c1 * c2 * c3) ** 2)
    for i in range(3):
        e = e + 0.5 * kap[i] * grad1d(c[i]) ** 2
    return float(e.sum())
 
 
def initial_profile(seed):
    x = np.arange(NX) - NX / 2
    c3 = seed * np.exp(-(x / D) ** 2)
    c1 = 0.5 * (1 - np.tanh(2 * x / D)) * (1 - c3)
    c2 = 0.5 * (1 + np.tanh(2 * x / D)) * (1 - c3)
    return np.array([c1, c2, c3])
 
 
def relax(kap, bet, closure, lam3=0.0, steps=8000, seed=1e-3):
    c = initial_profile(0.0 if closure == 'two-phase' else seed)
    for _ in range(steps):
        mu = chem_potential(c, kap, bet, lam3)
        if closure == 'two-phase':
            m = 0.5 * (mu[0] + mu[1])
            rhs = np.array([mu[0] - m, mu[1] - m, np.zeros(NX)])
        elif closure == 'independent':
            rhs = mu
        else:
            rhs = mu - mu.mean(axis=0)
        c = c - MOB * DT * rhs
        if not np.isfinite(c).all() or np.abs(c).max() > 1e3:
            return None, None, None
    return c, float(np.abs(c.sum(axis=0) - 1.0).max()), line_tension(c, kap, bet, lam3)
 
 
print("== pairwise sigma_ij -> per-component kappa_i (D = %.1f) ==" % D)
print("%-20s %-26s %-8s %s" % ("sigma_12,13,23", "kappa_1,2,3", "closes", "wedge angles [deg]"))
for s in [(1.0, 1.0, 1.0), (1.0, 0.8, 0.6), (1.0, 0.55, 0.55), (1.0, 0.4, 0.4)]:
    kap, _ = pair_to_component(*s)
    th = neumann_wedges(*s)
    print("%-20s %-26s %-8s %s" % (
        "%.2f %.2f %.2f" % s,
        "%+.2f %+.2f %+.2f" % tuple(kap),
        "yes" if kap.min() > 0 else "NO",
        "-" if th is None else "%.1f %.1f %.1f (sum %.1f)" % (th[0], th[1], th[2], sum(th))))
 
kap, bet = pair_to_component(1.0, 0.8, 0.6)
print()
print("== 1-2 interface at equilibrium, phase 3 seeded at 1e-3 (input sigma_12 = 1.00) ==")
print("%-30s %-13s %-9s %s" % ("closure", "max|sum c-1|", "max c_3", "sigma_eff"))
cases = [("A two-phase (c_3 frozen)", 'two-phase', 0.0),
         ("B three equations, no closure", 'independent', 0.0),
         ("C mean projection", 'projected', 0.0),
         ("D mean projection + Lambda=400", 'projected', 400.0)]
for tag, closure, lam3 in cases:
    c, drift, sig = relax(kap, bet, closure, lam3)
    if c is None:
        print("%-30s diverged" % tag)
        continue
    print("%-30s %-13.2e %-9.4f %.4f" % (tag, drift, c[2].max(), sig))
 
print()
print("== triple-penalty Lambda scan (mean projection) ==")
print("%-10s %-10s %-10s %s" % ("Lambda", "max c_3", "sigma_eff", "vs input"))
for lam3 in (0.0, 10.0, 30.0, 60.0, 150.0, 400.0):
    c, drift, sig = relax(kap, bet, 'projected', lam3)
    print("%-10.0f %-10.4f %-10.4f %+.1f%%" % (lam3, c[2].max(), sig, (sig - 1.0) * 100))

== pairwise sigma_ij -> per-component kappa_i (D = 4.0) ==
sigma_12,13,23       kappa_1,2,3                closes   wedge angles [deg]
1.00 1.00 1.00       +3.00 +3.00 +3.00          yes      120.0 120.0 120.0 (sum 360.0)
1.00 0.80 0.60       +3.60 +2.40 +1.20          yes      143.1 126.9 90.0 (sum 360.0)
1.00 0.55 0.55       +3.00 +3.00 +0.30          yes      155.4 155.4 49.2 (sum 360.0)
1.00 0.40 0.40       +3.00 +3.00 -0.60          NO       -
 
== 1-2 interface at equilibrium, phase 3 seeded at 1e-3 (input sigma_12 = 1.00) ==
closure                        max|sum c-1|  max c_3   sigma_eff
A two-phase (c_3 frozen)       1.78e-15      0.0000    0.9665
B three equations, no closure  8.59e-01      0.0000    0.0257
C mean projection              2.59e-13      0.2114    0.7809
D mean projection + Lambda=400 8.77e-15      0.0512    0.9715
 
== triple-penalty Lambda scan (mean projection) ==
Lambda     max c_3    sigma_eff  vs input
0          0.2114     0.7809     -21.9%
10         0.1887     0.8031     -19.7%
30         0.1448     0.8377     -16.2%
60         0.0963     0.8704     -13.0%
150        0.0515     0.9111     -8.9%
400        0.0512     0.9715     -2.9%

Mover los tres deslizadores de $\sigma$ cambia a la vez el signo de las barras de $\kappa$ y los ángulos de cuña. Al bajarlos a $(1.00,\,0.40,\,0.40)$, la columna closes de la tabla se da la vuelta.

Con $\sigma_{12}$ en $1.00$, al bajar $\sigma_{13}$ y $\sigma_{23}$ el triángulo se rompe en el instante en que la $\kappa_3$ roja cruza $0$. Esa es la fila $\kappa_3=-0.60$ de $(1.00,\,0.40,\,0.40)$. En $(1.00,\,0.80,\,0.60)$ los ángulos se quedan en $143.1^\circ$ / $126.9^\circ$ / $90.0^\circ$.

Σc aguanta a 1e-13 y la interfaz pierde 21.9%#

A congela la fase 3 por completo, así que es una corrida de dos fases. $\sigma_{\rm eff} = 0.9665$, un 3.3% por debajo del 1.00 de entrada. No es un defecto del modelo sino error de discretización con cuatro celdas a través de la interfaz. Confirma que la relación $\kappa_1 + \kappa_2 = \frac{3}{2}\sigma_{12}D$ se cumple.

B resuelve las tres ecuaciones de forma independiente. La suma deriva 0.859 y la tensión superficial se desvanece a 0.0257. La interfaz se derritió.

C es la línea interesante. La proyección promedio clava la suma en $2.6\times10^{-13}$. Pero $\max c_3 = 0.2114$ y $\sigma_{\rm eff} = 0.7809$, un 21.9% por debajo de la entrada. La restricción se cumple y la física está mal.

Una violación de restricción se anuncia sola. Un error del 21.9% no. Por eso es el caro.

La aritmética de una fase ausente que se filtra#

Quitando el término de gradiente, la fuga se ve de inmediato. Suponiendo simetría en el centro de la interfaz 1-2 y poniendo $c_3 = p$, $c_1 = c_2 = (1-p)/2$,

Derivando en $p$ y evaluando en $p=0$: el segundo término arranca en $p^2$, así que no aporta a la pendiente. Solo sobrevive el primero.

Siempre negativo, sin importar cuánto valgan $\sigma_{13}$ y $\sigma_{23}$. En un pozo cuártico por pares, añadir una tercera fase a una interfaz binaria es siempre cuesta abajo. Con $\beta_{12} = 3.0$ la pendiente es $-0.75$.

Los números concuerdan. En $p=0$, $E_0 = \beta_{12}/16 = 0.1875$. Con el $p = 0.2114$ medido, $E_0 = 0.1018$: una caída del 45.7%. La energía de volumen recibe con gusto a la fase 3. Lo único que la frena es $\kappa_3 |\nabla c_3|^2$, y con $\kappa_3 = 1.20$ el equilibrio solo llega en 0.21.

Hay un parecido con las corrientes parásitas, donde una gota que debería quedarse quieta empieza a fluir. Si la discretización de la tensión se desvía, la suma puede cuadrar y la interfaz equivocarse en silencio. Aquí $\sum c_i$ queda en $2.59\times10^{-13}$ y $\sigma_{\rm eff}$ en $0.7809$.

Dónde vive el penalizador triple Λ#

El remedio es un término que penalice directamente el solapamiento de tres fases en un punto.

En el barrido, subir $\Lambda$ empuja $\max c_3$ de 0.2114 a 0.0512 y devuelve $\sigma_{\rm eff}$ de 0.7809 a 0.9715. En $\Lambda = 400$ el déficit es −2.9%, prácticamente el mismo sitio que el error de discretización de A, −3.3%. Es indistinguible de congelar la fase 3.

Pero el término no tapa el agujero. Derivar $\Lambda c_1^2c_2^2c_3^2$ en $p$ también da cero en $p=0$. La pendiente cuesta abajo $-\beta_{12}/4$ sigue intacta. Lo que hace $\Lambda$ es levantar un muro después de que el bulto ya creció, y por eso $c_3$ se detiene en 0.05 en lugar de ir a cero por más que se suba $\Lambda$.

Para cerrar el agujero en sí hay que reducir a $N-1$ grados de libertad con $c_N = 1 - \sum_{i<N} c_i$, o simplemente no resolver la ecuación de una fase localmente ausente. A es ese caso, y por eso su $\max c_3$ es exactamente cero.

Si creer el 1.00 escrito o el 0.7809 integrado#

Meter la tabla por pares fija a la vez el signo de $\kappa_i = \frac{3D}{4}(\sigma_{ij}+\sigma_{ik}-\sigma_{jk})$ y los ángulos de cuña. Con $\sigma=(1.00,\,0.80,\,0.60)$ los coeficientes son $+3.60$, $+2.40$, $+1.20$ y los ángulos $143.1^\circ$ / $126.9^\circ$ / $90.0^\circ$, con suma $360^\circ$. Con $\sigma=(1.00,\,0.40,\,0.40)$, $\kappa_3=-0.60$ y el triángulo no cierra. Cinco fases dan diez ternas así. Una entrada negativa y la energía se abre hacia abajo en la celda donde esa terna se encuentra. Preguntar si los tres $\cos\theta_i$ suman $360^\circ$ es la misma pregunta que si cada $\kappa_i$ es positivo. Un ángulo se puede superponer al dibujo.

El valor de verificación no es la $\sigma_{12}=1.00$ escrita. Es la integral de C, $0.7809$. La suma aguanta en $2.59\times10^{-13}$ y la tensión queda un $21.9\%$ por debajo. Igual que en un paso de tiempo capilar que sigue presente cuando la tensión superficial se trata de forma implícita, ese $21.9\%$ entra directo en la presión de Laplace y en la velocidad de las ondas capilares. Una comprobación que solo mira $\sum c_i=1$ pasa. Incluso con $\Lambda=400$ el déficit es $-2.9\%$, o $0.9715$, a $0.0050$ del $0.9665$ de A con la fase congelada. Si el $1.00$ del campo de entrada queda como referencia, esa diferencia no llega al registro.