Una interfaz de dos celdas pierde el 39% de su tensión superficial — la anisotropía del esfuerzo de Korteweg

Un manómetro que cruza la interfaz#

Hay una interfaz plana líquido–vapor en reposo. Un manómetro se desplaza despacio desde el lado del vapor hacia el líquido. La aguja no se mueve. Una interfaz plana no tiene salto de presión de Laplace (el salto que produce la curvatura). Aun así, la tensión superficial existe. ¿Dónde se esconde?

Al girar el manómetro 90 grados aparece la respuesta. Medida en dirección paralela a la interfaz, la aguja baja, pero solo dentro de la interfaz. Este artículo deduce esa diferencia a partir del esfuerzo de Korteweg (el esfuerzo que generan los gradientes de densidad). Después mide, en una red LBM de espaciado unitario, cuánto de esa integral se pierde según el ancho de la interfaz. Con una interfaz de dos celdas desaparece el 39%.

El esfuerzo de un fluido en reposo no es una sola presión#

Al escribir la ecuación de momento de un LBM multifásico, el esfuerzo se divide en tres partes.

$p$ es la presión termodinámica, $\boldsymbol{\tau}_K$ el esfuerzo de Korteweg y $\boldsymbol{\tau}_{NS}$ el esfuerzo viscoso de Navier–Stokes. $\mathbf{S}^{d}$ es el tensor de tasa de deformación sin traza (trace-free rate-of-strain, deformación de corte pura sin cambio de volumen).

En reposo $\mathbf{S}^{d}=0$, así que el esfuerzo viscoso desaparece. Solo quedan $p$ y $\boldsymbol{\tau}_K$. En un fluido uniforme $\boldsymbol{\tau}_K$ también es cero. Solo aparece dentro de la interfaz, donde hay gradiente de densidad. Y distingue direcciones. Da un valor normal a la interfaz y otro paralelo a ella.

La simulación de abajo se puede manipular directamente. La línea punteada amarilla recorre la interfaz y muestra el esfuerzo sobre un elemento de fluido en esa posición.

Las flechas verdes (componente normal) miden lo mismo que el círculo punteado en toda la interfaz. Las flechas naranjas (componente tangencial) se acortan en el centro. Al subir κ la interfaz se ensancha y el área naranja, que es σ, crece como $\sqrt{\kappa}$.

De la energía libre al tensor de Korteweg#

Un modelo de interfaz difusa parte de un único funcional de energía libre.

$\psi(\rho)$ es la densidad de energía libre del fluido uniforme y $\kappa$ el coeficiente de energía de gradiente. La variación de este funcional da el tensor de presión.

$p_0=\rho\psi'-\psi$ es la presión de la ecuación de estado. El último término, $\kappa\nabla\rho\otimes\nabla\rho$, es el origen de la anisotropía. La divergencia de este tensor cumple $\nabla\cdot\mathbf{P}=\rho\nabla\mu$, donde $\mu=\psi'(\rho)-\kappa\nabla^2\rho$ es el potencial químico. La diferencia entre introducir una u otra forma como fuerza en LBM se trató en las dos formas del forcing en LBM no ideal.

Aquí se usa la ecuación de estado de van der Waals en unidades reducidas.

$\rho$, $T$ y $p$ están divididas por sus valores críticos. Para $T<1$ coexisten una densidad de vapor $\rho_g$ y una de líquido $\rho_l$. Las fijan las condiciones de Maxwell ($\mu$ y $p$ iguales en ambas fases).

La componente normal es plana; la tangencial se hunde#

Se toma la normal de la interfaz como eje $z$. La densidad depende solo de $z$. Las dos componentes del tensor son:

$p_N=P_{zz}$ es la componente normal y $p_T=P_{xx}$ la tangencial. El equilibrio estático exige $dP_{zz}/dz=0$. Por eso $p_N$ es igual a la presión del seno $p_b$ en toda la interfaz. Esa es la razón por la que el manómetro no se movió.

La diferencia es $p_N-p_T=\kappa\rho'^2\ge 0$. La presión tangencial es menor, y solo dentro de la interfaz. Una presión menor significa que la capa tira hacia los lados. Esa es la definición mecánica de la tensión superficial.

Al multiplicar la ecuación de equilibrio $\kappa\rho''=\psi'(\rho)-\mu_b$ por $\rho'$ e integrar se obtiene una primera integral: $\frac{\kappa}{2}\rho'^2 = W(\rho) \equiv \psi(\rho)-\mu_b\rho+p_b$. Al sustituirla, la integral espacial se convierte en una integral en densidad.

$W(\rho)$ es la colina de energía libre entre los dos estados del seno. La fórmula muestra de inmediato que $\sigma\propto\sqrt{\kappa}$. Es un valor de referencia que no sabe nada de la red.

Código: una interfaz de van der Waals en una red dz=1#

LBM fija el espaciado de red en 1. Entonces el ancho de la interfaz depende solo de κ y $T$. El código siguiente construye la interfaz de equilibrio por flujo de gradiente (gradient flow, una iteración que baja un poco en la dirección del residuo). Luego calcula $p_N$ y $p_T$ en la red y compara la σ de la integral de esfuerzos con la σ de la primera integral.

import math
 
# van der Waals reducido: rho_c = T_c = p_c = 1
def p0(r, T):  return 8*T*r/(3 - r) - 3*r*r
def mu0(r, T): return 8*T/3*(math.log(r/(3 - r)) + r/(3 - r)) - 6*r
def psi(r, T): return r*mu0(r, T) - p0(r, T)
 
def maxwell(T):
    # igual potencial químico e igual presión (Newton, jacobiano numérico)
    rg, rl = 1 - 2*math.sqrt(1 - T), 1 + 2*math.sqrt(1 - T)
    h = 1e-7
    for _ in range(50):
        f1 = mu0(rl, T) - mu0(rg, T)
        f2 = p0(rl, T) - p0(rg, T)
        a = (mu0(rl+h, T) - mu0(rl, T))/h; b = -(mu0(rg+h, T) - mu0(rg, T))/h
        c = (p0(rl+h, T) - p0(rl, T))/h;   d = -(p0(rg+h, T) - p0(rg, T))/h
        det = a*d - b*c
        rl -= ( d*f1 - b*f2)/det
        rg -= (-c*f1 + a*f2)/det
    return rg, rl
 
def korteweg_profile(T, kappa, n=100, dz=1.0, tol=1e-7):
    rg, rl = maxwell(T)
    mub, pb = mu0(rl, T), p0(rl, T)
    L = n*dz                                         # unidades de red: dz = 1 como en una red LB
    z = [(i + 0.5)*dz - L/2 for i in range(n)]
    r = [0.5*(rl + rg) + 0.5*(rl - rg)*math.tanh(zi*math.sqrt(1 - T)/math.sqrt(kappa)) for zi in z]
    dt = 0.2*min(dz*dz/kappa, 1.0)
    for it in range(1, 30001):
        e = [r[0]] + r + [r[-1]]                     # extremos con gradiente nulo
        res = [mu0(r[i], T) - kappa*(e[i+2] - 2*r[i] + e[i])/dz**2 - mub for i in range(n)]
        r = [r[i] - dt*res[i] for i in range(n)]
        if max(abs(x) for x in res) < tol:
            break
    e = [r[0]] + r + [r[-1]]
    pN, pT = [], []
    for i in range(n):
        d1 = (e[i+2] - e[i])/(2*dz)
        lap = (e[i+2] - 2*r[i] + e[i])/dz**2
        iso = p0(r[i], T) - kappa*r[i]*lap
        pN.append(iso + 0.5*kappa*d1*d1)             # componente normal del tensor de presión de Korteweg
        pT.append(iso - 0.5*kappa*d1*d1)             # componente tangencial
    sigma_stress = sum(a - b for a, b in zip(pN, pT))*dz
    # primera integral kappa/2 rho'^2 = W(rho)  ->  sigma = int sqrt(2 kappa W) d rho
    m = 4000
    sigma_exact = 0.0
    for k in range(m):
        rr = rg + (k + 0.5)*(rl - rg)/m
        W = max(psi(rr, T) - mub*rr + pb, 0.0)
        sigma_exact += math.sqrt(2*kappa*W)*(rl - rg)/m
    w = (rl - rg)/max(abs(e[i+2] - e[i])/(2*dz) for i in range(n))
    return dict(w=w, it=it, rg=rg, rl=rl, spread=max(abs(p - pb) for p in pN),
                dip=min(pT) - pb, s_st=sigma_stress, s_ex=sigma_exact)
 
print("  T   kappa  width  iters  max|pN-pb|  min(pT-pb)  sigma_stress  sigma_exact  ratio")
for T in (0.85, 0.95):
    for kappa in (0.1, 0.5, 2.0, 8.0):
        d = korteweg_profile(T, kappa)
        print(f"{T:.2f}  {kappa:4.1f}  {d['w']:5.2f}  {d['it']:5d}   {d['spread']:.2e}   {d['dip']:+.4f}"
              f"     {d['s_st']:.5f}       {d['s_ex']:.5f}    {d['s_st']/d['s_ex']:.4f}")

La salida es:

  T   kappa  width  iters  max|pN-pb|  min(pT-pb)  sigma_stress  sigma_exact  ratio
0.85   0.1   2.08     47   2.09e-02   -0.0309     0.10067       0.16545    0.6085
0.85   0.5   2.63    275   1.08e-02   -0.1550     0.32864       0.36995    0.8883
0.85   2.0   4.35  30000   9.00e-03   -0.2429     0.71074       0.73990    0.9606
0.85   8.0   8.21    414   2.52e-03   -0.2650     1.46491       1.47980    0.9899
0.95   0.1   2.39    254   2.67e-03   -0.0112     0.02759       0.03239    0.8518
0.95   0.5   3.99  30000   9.96e-04   -0.0255     0.06887       0.07243    0.9508
0.95   2.0   7.38    290   3.61e-04   -0.0290     0.14301       0.14485    0.9873
0.95   8.0  14.48   1168   9.36e-05   -0.0298     0.28878       0.28971    0.9968

El mismo código puede ejecutarse directamente en el navegador. El valor por defecto es la segunda fila de la salida ($T=0.85$, $\kappa=0.5$).

Con κ en 0.1 el punto de abajo a la derecha cae cerca de un cociente de 0.61. Con κ en 8 se queda en 0.99. Cada corrida terminada deja un punto, así que conviene cambiar $T$ y completar combinaciones que la tabla de salida no tiene.

Una interfaz de dos celdas pierde el 39% de σ#

La columna ratio es la clave. Con un ancho de interfaz (diferencia de densidad entre pendiente máxima) de 2.08 celdas vale 0.6085. La red perdió el 39% de la tensión superficial. Con 2.63 celdas vale 0.888, con 4.35 celdas 0.961 y con 8.21 celdas 0.990. El cociente lo fija casi solo el ancho, sin importar la temperatura. A $T=0.95$, un ancho de 3.99 celdas da 0.951, cerca del caso de 4.35 celdas a $T=0.85$.

La columna min(pT-pb) cuenta lo mismo. A $T=0.85$, el hundimiento tangencial se profundiza de -0.031 a -0.265 al crecer κ. Se acerca al valor del continuo. En una interfaz delgada, el fondo del hundimiento queda escondido entre nodos.

La columna max|pN-pb| indica cuán lejos está la presión normal de ser plana. En el continuo sería cero. Con 2.08 celdas de ancho vale $2.1\times10^{-2}$, dos tercios del propio hundimiento tangencial. Ese gradiente residual actúa como una fuerza falsa sobre una interfaz estática. Es una de las fuentes de las corrientes espurias (spurious currents, flujo artificial alrededor de una interfaz que debería estar quieta).

Las dos filas con iters en 30000 no son fallas de convergencia. Una interfaz en red no es simétrica ante traslaciones. El perfil se desliza muy despacio hacia una posición anclada a los nodos. Mientras el residuo se estanca cerca de $10^{-7}$, σ no cambia en sus cinco primeras cifras.

√κ y (1−T)^1.5#

En la columna sigma_exact, cada vez que κ se cuadruplica σ se duplica: 0.165 → 0.370 → 0.740 → 1.480. Es el $\sqrt{\kappa}$ de la fórmula de la primera integral. El ancho también escala con $\sqrt{\kappa}$. Un κ mayor da una interfaz más gruesa, pero también sube la tensión superficial.

La dependencia con la temperatura sigue el exponente crítico de campo medio. Con κ=0.5, el cociente de σ entre $T=0.85$ y $T=0.95$ es 0.36995/0.07243 = 5.11. Frente a $(0.15/0.05)^{3/2}=5.20$, la coincidencia es buena. Cerca del punto crítico la diferencia de densidades disminuye y la interfaz se ensancha. Por eso las filas de $T=0.95$ dan mejor cociente con el mismo κ.

En la práctica ambos efectos chocan. Elegir una σ y una razón de densidades objetivo fija κ y $T$, y con ello el ancho. Por debajo de cuatro celdas, σ cae más de un 5%. Los modelos de campo de fase conservativos tienen el mismo límite de ancho. El espesor de interfaz W de el término de curvatura de Allen–Cahn conservativo se mantiene en cuatro o cinco celdas o más por la misma razón.

Medir el ancho antes de fijar κ#

Una vez elegido κ para una corrida LBM multifásica, conviene relajar primero una sola interfaz plana en 1D. Hay que revisar dos números. El primero es el cociente entre el $\kappa\sum\rho'^2$ de la red y $\int\sqrt{2\kappa W}\,d\rho$. El segundo es $\max|p_N-p_b|$. Si el cociente queda por debajo de 0.95, la σ de la prueba de Laplace saldrá baja en la misma proporción. En ese caso es mejor subir κ o acercar $T$ al punto crítico. Si la presión normal no es plana, en la prueba de gota estática lo primero que hay que sospechar son las corrientes espurias.