Cuando la malla se sesga, el laplaciano se fuga — corrección de flujo difusivo no ortogonal

Un solver de difusión que va bien en cubos empieza a mentir sobre triángulos sesgados. Misma ecuación. Mismo código. El campo de temperatura se emborrona y la convergencia se arrastra. No son los datos de contorno. No es el material. Es la forma de medir el gradiente en la cara.

Al empujar el deslizador de no ortogonalidad de 20° a 65°, el vector cian de la sobrerrelajación se alarga más que el área, y el vector rosa de la corrección mínima se alarga todavía más.

Cuando el cian queda más corto que el rosa, la lectura de ρ bajo el lienzo cruza 1. A 45° solo la corrección mínima debería ponerse roja primero.

El gradiente en la cara es el problema#

En el método de volúmenes finitos, la difusión es una suma de flujos de cara. El flujo por una cara $f$ se escribe así.

$\Gamma_f$ es la difusividad en la cara. $\nabla\phi_f$ es el gradiente de cara. $\mathbf{S}_f$ es el vector de área saliente, área por normal.

Lo difícil es bajar $\nabla\phi_f \cdot \mathbf{S}_f$ a valores de celda. En una malla ortogonal, la línea $\mathbf{d}$ entre centros va con la normal. $(\phi_N-\phi_P)/|\mathbf{d}|$ es el gradiente normal.

En una malla no estructurada esa hipótesis se rompe. $\mathbf{d}$ y $\mathbf{S}_f$ se desalinean. El ángulo es la no ortogonalidad $\theta$. $(\phi_N-\phi_P)/|\mathbf{d}|$ pasa a ser el gradiente a lo largo de $\mathbf{d}$, no a lo largo de la normal. Ignorar la brecha hace que la difusión se fugue. El campo se emborrona.

Partir S_f en dos piezas#

El vector de área se parte en dos. Una parte a lo largo de $\mathbf{d}$, otra con el resto.

$\mathbf{E}_f$ es la parte a lo largo de $\mathbf{d}$. $\mathbf{T}_f$ es el resto tangencial. El flujo queda partido.

El primer término solo usa $\phi_N-\phi_P$. Entra en la matriz de forma implícita. Dos centros de celda, nada más. El segundo término necesita el gradiente promedio de cara $\overline{\nabla\phi}_f$, un promedio de gradientes vecinos. Es un mal coeficiente de matriz, así que va al lado derecho como corrección explícita. Rehacer ese gradiente por mínimos cuadrados es otro trabajo, escrito en pesos que mantienen vivo un gradiente en una malla no estructurada. Cuando un lazo diferido oscila, el culpable a veces es ese gradiente, no ρ.

Si $\theta=0$, entonces $\mathbf{T}_f=0$. Vuelve la vieja diferencia central.

Mínima, ortogonal, sobrerrelajada: tres particiones#

La longitud de $\mathbf{E}_f$ a lo largo de $\mathbf{d}$ no es única. Con $\hat{\mathbf{d}}=\mathbf{d}/|\mathbf{d}|$ y $\cos\theta=\hat{\mathbf{d}}\cdot\mathbf{S}_f/|\mathbf{S}_f|$, hay tres opciones.

• Mínima: $|\mathbf{E}_f|=|\mathbf{S}_f|\cos\theta$. El $\mathbf{T}_f$ más corto.
• Ortogonal: $|\mathbf{E}_f|=|\mathbf{S}_f|$. La magnitud implícita se queda en el área de la cara.
• Sobrerrelajada: $|\mathbf{E}_f|=|\mathbf{S}_f|/\cos\theta$. La parte implícita crece cuando la malla se sesga.

En la figura de arriba, θ a 60° hace que el $\mathbf{E}_f$ sobrerrelajado sea más largo que el vector de área. El coeficiente diagonal creció. La corrección mínima hace lo contrario: $\mathbf{E}_f$ se encoge y $\mathbf{T}_f$ se estira.

Código: un solver de Laplace con corrección diferida#

El término explícito usa el gradiente del pase anterior. Una sola resolución no basta. La corrección se refresca a lo largo de varios barridos. Eso es la corrección diferida.

import numpy as np
 
def decompose_face(S_f, d, scheme="over"):
    """Vector de area S_f = E_f (implicito, segun d) + T_f (correccion explicita)."""
    d_hat = d / np.linalg.norm(d)
    Smag = np.linalg.norm(S_f)
    dS = d_hat @ S_f                      # = |S_f| cos(no ortogonalidad)
    if scheme == "min":                   # correccion minima
        E = dS * d_hat
    elif scheme == "orth":                # correccion ortogonal
        E = Smag * d_hat
    else:                                 # sobrerrelajada (over-relaxed)
        E = (Smag**2 / dS) * d_hat
    return E, S_f - E                     # E_f, T_f
 
def sweep_residual(S_f, d, scheme, sweeps=25):
    """E_f implicito, T_f diferido. El residual se contrae a razon rho = |T|/|E|."""
    E, T = decompose_face(S_f, d, scheme)
    rho = np.linalg.norm(T) / np.linalg.norm(E)
    r, hist = 1.0, [1.0]
    for _ in range(sweeps):
        r *= rho                          # una barrida de correccion diferida = residual x rho
        hist.append(r)
    return rho, hist
 
S_f = np.array([1.0, 0.0])                # normal de cara (vector de area, |S|=1)
for ang in (20, 45, 65):
    d = np.array([np.cos(np.radians(ang)), np.sin(np.radians(ang))])
    for s in ("min", "orth", "over"):
        rho, hist = sweep_residual(S_f, d, s)
        tag = "converge" if rho < 1 else "diverge"
        print(f"theta={ang:2d}  {s:5s}  rho={rho:.3f}  {tag}  (tras 25: {hist[-1]:.1e})")

La corrida dice que los tres esquemas se separan al crecer el ángulo. A 20° los tres son suaves. A 65° la corrección mínima diverge con $\rho=\tan 65°\approx 2.14$, mientras la sobrerrelajada todavía contrae con $\rho=\sin 65°\approx 0.91$. Misma malla, misma ecuación. La partición decide la convergencia.

Los tres cocientes se cierran solo con el ángulo. Con área 1, la corrección mínima tiene longitud de E igual a cosθ y longitud de T igual a sinθ, así que ρ es tanθ. La tangente de 45° es 1, y el lazo diferido multiplica el residuo por 1. Veinticinco barridos no cambian nada. La tangente de 20° anda cerca de 0.364. La de 65°, cerca de 2.14.

La corrección ortogonal tiene longitud de E igual a 1 y longitud de T igual a la raíz de 2 menos 2cosθ. ρ llega a 1 cuando cosθ es un medio, a 60°. A 65° ya anda cerca de 1.07. Llamarla insegura solo a 70° llega tarde.

El ρ sobrerrelajado es sinθ, cerca de 0.906 a 65°. El fragmento multiplica el residuo por ese cociente 25 veces. Elevado a 25, 0.906 da cerca de 0.085, así que queda cerca de un 8 por ciento del error explícito. A 20°, con seno cerca de 0.342, la misma potencia cae alrededor de $10^{-12}$. Dos o tres correctores a 65° son 0.906 al cubo, cerca de 0.74. El muro de 45° es exacto, no empírico.

¿Por qué ρ es la longitud de T sobre la longitud de E? El coeficiente implícito escala con E. Cada barrido devuelve un error que escala con T. El residuo se encoge o crece por ese cociente. El lazo se cierra solo cuando el cociente queda por debajo de 1.

Qué pasa cuando crece la no ortogonalidad#

Los tres residuos corren en un solo lienzo. El deslizador del ángulo fija ρ, y cada barrido baja la curva por ese cociente.

Por debajo de 30° las tres curvas caen juntas. A 45° solo la corrección mínima se pone horizontal. Pasados los 60°, la corrección ortogonal también cruza ρ de 1, y solo queda la sobrerrelajada, sentada en sinθ.

Hay un costo. La sobrerrelajación agranda a la vez la parte implícita y la corrección explícita. Llegar a un residuo estrecho pide más barridos. nNonOrthogonalCorrectors es ese conteo. Una malla limpia quiere 0 o 1. Una malla muy sesgada se empuja a 2 o 3. A 65° esos dos o tres barridos todavía dejan cerca de tres cuartos del error. Sumar más encarece el paso de tiempo. Recortar el ángulo casi siempre sale más barato.

Por qué la corrección mínima ya no sirve a 45°#

A 45° la corrección mínima deja de reducir el residuo, porque ρ es exactamente 1. La corrección ortogonal choca con el mismo muro a 60° y ya lo pasó a 65°. La sobrerrelajada se mantiene bajo 1 hasta cerca de 90°, porque ρ es sinθ. El error que queda dentro de un paso es el precio.

¿Sumar barridos, o rehacer la malla antes de que el ángulo pase de 60°? El orden en que los triángulos cubren los puntos fija ese ángulo. Dos formas de cubrir una nube de puntos con triángulos pregunta dónde ese orden deja no ortogonalidad atrás. ¿Sale más barato alargar el lazo de correctores, o empujar el frente otra vez?