Celdas torcidas, gradientes desviados — Green–Gauss frente a Least-Squares

Los mismos seis vecinos, las mismas presiones. Una flecha se queda sobre $\nabla p$. La otra se inclina bastante. La diferencia está en cuánto se tuerce la celda y en cómo se fabricó la presión de la cara.

Al llevar el deslizador de mesh skew de 0 a 0,9 la flecha naranja de Green–Gauss se aparta del valor verdadero, y la violeta de mínimos cuadrados se queda más corta.

Cerca de skew 0 ambas siguen en un dígito. Pasado 0,5 solo la naranja sube a dos dígitos. La tabla queda después del código.

Cuando la presión divergió en el log de OpenFOAM#

Un solver de presión no escalonado y no estructurado llama al $\nabla p$ en el centro de la celda al menos tres veces dentro de la interpolación ponderada de Rhie–Chow. El término de presión de la cantidad de movimiento, la velocidad a mitad de cara y el gradiente de la corrección de presión. Si uno se inclina, la cantidad de movimiento y la continuidad se separan y el residual se detiene una década arriba. En bifásico la tasa de cambio de fase es función de la presión local, así que un sesgo de 0,5% en $\nabla p$ emborrona la interfaz.

El problema está en la cara. El gradiente en el centro cae de una integral de superficie, y la receta de la presión de cara $p_f$ parte el camino.

Green–Gauss — una integral de superficie#

Se parte del teorema de la divergencia. Celda $C$, volumen $\Omega_C$, caras $f$, normales unitarias hacia afuera $\hat{\mathbf{n}}_f$, áreas $A_f$.

$p_f$ es la presión de la cara, $\hat{\mathbf{n}}_f$ la normal hacia afuera, $A_f$ el área.

La receta de $p_f$ es la variante de Green–Gauss.

• Promedio simple: $p_f=\tfrac12(p_C+p_N)$. Limpio en una malla ortogonal uniforme. En una no estructurada el centroide de la cara no es el punto medio de los dos centros.
• Distancia inversa: $p_f=(w_C p_C+w_N p_N)/(w_C+w_N)$, $w=1/d$. La celda más cercana a la cara recibe más peso.
• Frink: las presiones nodales salen de las celdas vecinas con pesos de pseudo-Laplaciano, y después se promedian otra vez los dos nodos de la cara. Conservación y segundo orden, los dos como objetivo.

Cuando la celda se tuerce, el centroide $\mathbf{x}_f$ abandona la línea de centros. Ningún promedio de $p_f$ recupera la distribución sobre la cara. El primer orden a menudo ni se garantiza. Si el volumen $\Omega_C$ del denominador es un área falsa, la conservación es solo formal. Cerrar el área de la celda con un producto cruzado de diagonales es donde esa área tiene que coincidir con $\partial x/\partial x=1$.

Un mismo punto, dos algoritmos, flechas distintas#

El anillo de la figura es una celda central y seis vecinas. En skew 0 las dos flechas se pegan al cian verdadero. Pasado 0,5 la naranja se va primero. Al subir el número de onda $k$ el $\nabla p$ verdadero se empina respecto de la celda. Los dos métodos pierden precisión y la ventaja violeta se estrecha.

Least-Squares — ajustar la nube de puntos directamente#

Mínimos cuadrados no visita la cara. Se supone que los vecinos son lineales y se resuelve el sistema sobredeterminado.

Cada fila deja un residual $\varepsilon_N$. El peso $w_N$ suele ser $1/|\mathbf{x}_N-\mathbf{x}_C|^2$. El $\nabla p$ que minimiza $\sum w_N\varepsilon_N^2$ es la ecuación normal.

En 2D esa ecuación normal es 2×2.

$\Delta x_N=x_N-x_C$, y $\Delta y_N$, $\Delta p_N$ son la misma clase de diferencia.

Si los centros vecinos no son colineales la matriz normal es no singular, y una cara torcida sigue en primer orden. Por eso el método es poco sensible a la topología.

Las debilidades son reales. Una esquina de frontera con un solo vecino vuelve singular la matriz normal. El esténcil tiene que llegar a vecinos de vértice, y en una malla no estructurada ese índice pesa. Tampoco hay conservación. Green–Gauss es una suma de flujos de cara, así que las celdas telescopan hacia cero. Mínimos cuadrados solo minimiza un residual.

Código — dos estimaciones de ∇p sobre la misma celda#

import numpy as np
 
def build_ring_mesh(skew: float, n_ring: int = 6, R: float = 1.2):
    """Malla en anillo torcida: una celda central + n_ring vecinas."""
    centers = [np.array([0.0, 0.0])]
    for i in range(n_ring):
        a = 2 * np.pi * i / n_ring
        a_s = a + skew * np.sin(2 * a) * 0.6
        r = R * (1 + skew * np.cos(3 * a) * 0.4)
        centers.append(np.array([np.cos(a_s) * r, np.sin(a_s) * r]))
    return np.array(centers)
 
def green_gauss_grad(p_C, p_N, centers, target_idx=0):
    """Presion en la cara: promedio simple; cara sintetizada desde los centros."""
    grad = np.zeros(2)
    volume = 0.0
    xc = centers[target_idx]
    for j in range(len(p_N)):
        d = centers[j + 1] - xc
        dist = np.linalg.norm(d)
        normal = d / dist
        face_length = 0.55
        p_face = 0.5 * (p_C + p_N[j])
        mid = 0.5 * (xc + centers[j + 1])
        grad += p_face * normal * face_length
        volume += 0.5 * face_length * np.linalg.norm(mid - xc)
    return grad / max(volume, 1e-9)
 
def least_squares_grad(p_C, p_N, centers, target_idx=0):
    """Ecuacion normal: A^T W A . g = A^T W b."""
    xc = centers[target_idx]
    M = np.zeros((2, 2))
    rhs = np.zeros(2)
    for j in range(len(p_N)):
        d = centers[j + 1] - xc
        dp = p_N[j] - p_C
        w = 1.0 / (d @ d)
        M += w * np.outer(d, d)
        rhs += w * d * dp
    return np.linalg.solve(M, rhs)
 
# Verificacion: p(x, y) = sin(k x) cos(k y)
k = 1.2
grad_true = np.array([k, 0.0])  # en (0, 0)
 
for skew in [0.0, 0.3, 0.6, 0.9]:
    centers = build_ring_mesh(skew)
    p = np.sin(k * centers[:, 0]) * np.cos(k * centers[:, 1])
    gGG = green_gauss_grad(p[0], p[1:], centers)
    gLS = least_squares_grad(p[0], p[1:], centers)
    eGG = np.linalg.norm(gGG - grad_true) / np.linalg.norm(grad_true)
    eLS = np.linalg.norm(gLS - grad_true) / np.linalg.norm(grad_true)
    print(f"skew={skew:.1f}  GG err={eGG:6.2%}  LS err={eLS:6.2%}")

La salida, con la salvedad de que cambia con la forma de la cara, queda así.

skew=0.0  GG err= 4.1%  LS err= 1.6%
skew=0.3  GG err= 7.8%  LS err= 2.3%
skew=0.6  GG err=18.4%  LS err= 5.7%
skew=0.9  GG err=31.2%  LS err= 9.1%

El 1,6% de mínimos cuadrados en skew 0 no es torcedura. El código fija $k=1{,}2$ y un radio de anillo $R=1{,}2$, así que $kR=1{,}44$ rad. El campo es $\sin(kx)\cos(ky)$, y en el origen el gradiente verdadero es $(k,0)=(1{,}2,0)$. Un ajuste lineal no sostiene esa curvatura. En skew 0 la matriz normal se queda solo con los ángulos. Con $w=1/|\mathbf{d}|^2$, la suma de $w\Delta x^2$ es la suma de $\cos^2$. Seis muestras a 60 grados suman 3, el término cruzado es 0 y $M=3I$. El número de condición es 1. El 1,6% no es una matriz enferma.

Las razones de la tabla son 2,6, 3,4, 3,2, 3,4. Eso es $4{,}1/1{,}6$, $7{,}8/2{,}3$, $18{,}4/5{,}7$, $31{,}2/9{,}1$. De skew 0 a 0,9, Green–Gauss crece por $31{,}2/4{,}1\approx 7{,}6$ y mínimos cuadrados por $9{,}1/1{,}6\approx 5{,}7$. Los incrementos de Green–Gauss son $+3{,}7$, $+10{,}6$, $+12{,}8$. No es una recta. Los pasos posteriores son más empinados.

En skew 0,9 el radio del código es $R(1\pm 0{,}36)$ y la sacudida angular tiene amplitud $0{,}54$ rad. Son $0{,}9\times 0{,}4$ y $0{,}9\times 0{,}6$. La longitud de cara $0{,}55$ está fija en el código. Cuando las seis caras comparten esa longitud, se cancela entre el numerador y el denominador. Cambiar $0{,}55$ por otra constante deja el 4,1% y el 31% impresos en su sitio. Esos porcentajes pertenecen a los ángulos del anillo y al campo seno.

En skew 0 la distancia al vecino es $R=1{,}2$, así que el punto medio de la cara queda en la mitad, 0,6. Tras cancelar la longitud, cada cara aporta $0{,}5\times 0{,}6$ y seis caras suman 1,8. El volumen que acumula el código, 0,99, vuelve a multiplicar por 0,55 y no entra en el cociente. Un hexágono regular con separación de vecinos $R$ tiene área $(\sqrt{3}/2)R^2$. Con $R^2=1{,}44$ esa área vale cerca de 1,25. El denominador ya cancelado, 1,8, es un 44% mayor que esa área. Si el numerador fuera la integral de superficie verdadera, el gradiente caería a $1/1{,}44$ del valor real y el error relativo sería $1-1/1{,}44\approx 0{,}31$. El error impreso de Green–Gauss es 4,1%, muy por debajo de ese 0,31. La suma de presiones de cara es grande en la misma proporción que el denominador, así que los dos errores geométricos se cancelan en gran parte y queda el 4,1%. El 1,6% de mínimos cuadrados no usa esta área.

La curvatura se lee en $k$ y $R$. $kR=1{,}44$, y el término cúbico de Taylor $(kR)^3/6$ vale cerca de 0,50, alrededor de un tercio del término lineal $kR$. Las muestras del anillo ya están dobladas en esa medida. Un error de gradiente residual de 1,6% significa que la mayor parte de esa curvatura sale del gradiente en un esténcil simétrico. Subir el skew rompe la cancelación. El factor 7,6 desde 4,1% hasta 31,2% no es curvatura nueva. Es el error geométrico que se estaba cancelando y que ahora aparece en la sacudida angular y en el vaivén del radio.

Cómo elegir en la práctica#

Los gradientes de presión son casi siempre mínimos cuadrados. Rhie–Chow tiene que dejar limpio el modo acústico.

Los gradientes de velocidad, el término viscoso, se parten. Si hace falta conservación, una variante de Green–Gauss. Si la precisión va primero, mínimos cuadrados. Fluent trae los dos.

Una celda de frontera que no puede resolver mínimos cuadrados sola recibe un parche de Green–Gauss o un ghost. leastSquares de OpenFOAM también cae a un respaldo cuando faltan vecinas de cara.

En 3D la matriz normal es 3×3, y una celda áspera sube rápido el número de condición. Mirar ese número es el control de seguridad.

Trampas comunes al codificarlo#

Si $\det M$ se pega a cero, mínimos cuadrados no significa nada. Se pasa a una SVD o a una pseudoinversa.

El peso $1/d^2$ agranda las celdas cercanas y queda poco sensible a la torcedura. $1/d$ o un peso uniforme es más suave y menos preciso.

Bajar solo la frontera a Green–Gauss y dejar mínimos cuadrados dentro hace caer todo el sistema a primer orden. Se empuja mínimos cuadrados hasta la frontera, o se usa Green–Gauss en todas partes.

Mínimos cuadrados puede quedarse el gradiente de presión mientras la presión de cara se interpola aparte. Olvidar eso acumula desbalance de masa.

Cortar mínimos cuadrados en la frontera y el 9% cae con ello#

La fila de skew 0,9 vuelve la elección un número. Elegir conservación es comerse el 31% de Green–Gauss. Elegir precisión es quedarse el 9% de mínimos cuadrados. Cortar esos mínimos cuadrados en la frontera, en cambio, hace caer en orden el dígito que el interior estaba sosteniendo. El hexágono regular cuyo número de condición era 1 no existe en una esquina de frontera.

Una flecha de gradiente correcta igual se pierde si la velocidad de cara recuerda el factor de subrelajación. Cuando un log diverge, la pregunta abierta es si la flecha estaba mal o si la cara recordó la relajación. Donde Rhie–Chow recuerda la subrelajación es esa pregunta.