Subir el número de Péclet de celda a 2 llevó la solución hasta -0.33 — donde desapareció la energía a minimizar

Cuatro métodos dieron la misma respuesta en una sola barra#

El capítulo 1 de unos apuntes de elementos finitos resuelve el mismo problema de barra cuatro veces: rigidez directa, mínima energía potencial total, residuos ponderados y Galerkin. Las cuatro producen la misma matriz de rigidez de 5×5. Los apuntes lo dicen sin rodeos: se use el método que se use, el resultado apenas cambia y el valor es preciso.

Esa frase lleva una condición adosada, y la condición se esconde dentro del problema de barra, donde nadie tropieza con ella. Las ecuaciones del flujo no la cumplen. Lo que se mide es qué garantía desaparece cuando la condición se rompe, y qué parte puede devolver la estabilización. Al subir advection a desde cero, el trabajo de la vuelta del punto rosa pasa de 0.000 a exactamente $2\pi a$.

La bola ámbar que se deslizaba hacia dentro sin cruzar las curvas de nivel grises empieza a remontarlas en espiral cuando a crece. Al empujar diffusion s hasta el tope, el total de la vuelta se queda quieto.

Para poder minimizar, la matriz tiene que ser simétrica#

La formulación de mínima energía potencial deriva la energía potencial total respecto de las incógnitas nodales y la iguala a cero. En forma discreta queda así.

$\boldsymbol{\phi}$ es el vector de incógnitas nodales, $\mathbf{K}$ la matriz de rigidez y $\mathbf{f}$ el vector de cargas. Ahora se deriva respecto de la componente $i$.

Conviene fijarse en lo que apareció: $K_{ij} + K_{ji}$. Para que esta expresión se reduzca a $\mathbf{K}\boldsymbol{\phi} = \mathbf{f}$ hace falta que $K_{ij} = K_{ji}$. Sin simetría, lo que la minimización resuelve en realidad es la matriz simetrizada $\tfrac{1}{2}(\mathbf{K}+\mathbf{K}^{T})$, que ya no es la ecuación original.

Hay una versión más profunda del mismo enunciado. Para que el campo de residuos $\mathbf{r}(\boldsymbol{\phi}) = \mathbf{f} - \mathbf{K}\boldsymbol{\phi}$ sea el gradiente de alguna función escalar, su jacobiano debe ser simétrico. Si no lo es, la función potencial sencillamente no existe. La forma más rápida de comprobarlo es recorrer un camino cerrado: un campo gradiente siempre devuelve trabajo cero tras una vuelta completa. Cuando a deja el cero y la bola remonta las curvas de nivel, esa vuelta ya no cierra en cero.

La advección rompe exactamente una cantidad $u$#

Al escribir la forma débil de la ecuación de advección-difusión unidimensional $u\,\phi' = \epsilon\,\phi''$, los dos términos muestran caracteres distintos.

$w$ es la función de prueba, $\epsilon$ la difusividad y $u$ la velocidad de advección. El primer término no cambia al intercambiar $w$ y $\phi$: es simétrico. El término de advección invierte el signo al integrar por partes. Para funciones de prueba que se anulan en los bordes, $\int w\,u\,\phi'\,dx = -\int \phi\,u\,w'\,dx$, de modo que la advección es una contribución puramente antisimétrica.

Al ensamblar con elementos lineales, esa estructura sobrevive en los coeficientes. Con tamaño de elemento $h$, una fila interior queda así.

La difusión deposita el mismo valor en ambos vecinos; la advección añade $+u/2$ de un lado y $-u/2$ del otro. Por eso la desviación de la simetría es exactamente un número.

Por mucho que se refine la malla, este valor sigue siendo $u$, porque $h$ nunca entró en él. Mientras haya flujo, el principio de mínima energía potencial no va a volver. Plegar la rigidez en una matriz de Voigt de 6×6 solo funcionaba mientras la matriz seguía simétrica. Cuando la advección hace oscilar un esquema de segundo orden junto a un escalón, la vía está en la advección TVD con un limitador de flujo.

Tres esquemas sobre la misma malla en Python#

Diez elementos, $u = 1$, condiciones de contorno $\phi(0)=0$ y $\phi(1)=1$. Un único coeficiente de difusión artificial $\beta$ genera los tres esquemas mediante $\epsilon_{\text{eff}} = \epsilon + \beta\,u\,h/2$: con $\beta = 0$ se tiene Galerkin puro, con $\beta = 1$ upwind completo y con $\beta = \coth(Pe_h) - 1/Pe_h$, SUPG.

import numpy as np
 
def assemble_ad(n, vel, eps, beta):
    """Adveccion-difusion 1D en elementos lineales. beta = difusion artificial."""
    h = 1.0 / n
    eps_eff = eps + beta * vel * h / 2.0
    kd = (eps_eff / h) * np.array([[1.0, -1.0], [-1.0, 1.0]])   # difusion: simetrica
    ka = (vel / 2.0) * np.array([[-1.0, 1.0], [-1.0, 1.0]])     # adveccion: antisimetrica
    K = np.zeros((n + 1, n + 1))
    for e in range(n):
        K[e:e + 2, e:e + 2] += kd + ka
    return K
 
def solve_bvp(K):
    n = K.shape[0] - 1
    A, b = K.copy(), np.zeros(n + 1)
    A[0, :], A[0, 0], b[0] = 0.0, 1.0, 0.0
    A[n, :], A[n, n], b[n] = 0.0, 1.0, 1.0
    return np.linalg.solve(A, b)
 
def exact_ad(x, pe):
    return (np.exp(pe * (x - 1.0)) - np.exp(-pe)) / (1.0 - np.exp(-pe))
 
def skew_ratio(K):
    return np.linalg.norm(K - K.T) / np.linalg.norm(K + K.T)
 
def loop_work(K, m=20000):
    """Trabajo del campo residual -K.phi al recorrer el circulo unidad en el espacio de GDL."""
    t = np.linspace(0.0, 2.0 * np.pi, m, endpoint=False)
    path = np.stack([np.cos(t), np.sin(t)])          # posicion
    tang = np.stack([-np.sin(t), np.cos(t)])         # dl / dt
    return float(np.sum(np.sum(-(K @ path) * tang, axis=0)) * (2.0 * np.pi / m))
 
n, vel = 10, 1.0
x = np.linspace(0.0, 1.0, n + 1)
print(f"{'Pe_h':>5} {'scheme':>9} {'max|err|%':>10} {'min phi':>9} {'skew/sym':>9} {'loop W':>8}")
for pe_h in [0.5, 1.0, 2.0, 5.0]:
    eps = vel / (2.0 * n * pe_h)
    ex = exact_ad(x, vel / eps)
    for name, beta in [("Galerkin", 0.0), ("upwind", 1.0),
                       ("SUPG", 1.0 / np.tanh(pe_h) - 1.0 / pe_h)]:
        K = assemble_ad(n, vel, eps, beta)
        phi = solve_bvp(K)
        print(f"{pe_h:5.1f} {name:>9} {100 * np.max(np.abs(phi - ex)):10.2f} "
              f"{phi.min():9.4f} {skew_ratio(K):9.4f} {loop_work(K[4:6, 4:6]):8.4f}")
 
K0 = assemble_ad(n, 0.0, 0.1, 0.0)
print(f"\nvel = 0 : skew/sym = {skew_ratio(K0):.2e},  loop W = {loop_work(K0[4:6, 4:6]):.2e}")
print(f"pi * u  = {np.pi * vel:.4f}")

 Pe_h    scheme  max|err|%   min phi  skew/sym   loop W
  0.5  Galerkin       3.45    0.0000    0.2924   3.1416
  0.5    upwind      13.17    0.0000    0.1954   3.1416
  0.5      SUPG       0.00   -0.0000    0.2704   3.1416
  1.0  Galerkin      13.53    0.0000    0.5774   3.1416
  1.0    upwind      19.80    0.0000    0.2924   3.1416
  1.0      SUPG       0.00   -0.0000    0.4428   3.1416
  2.0  Galerkin      35.17   -0.3334    1.1010   3.1416
  2.0    upwind      18.17   -0.0000    0.3885   3.1416
  2.0      SUPG       0.00    0.0000    0.5572   3.1416
  5.0  Galerkin      69.61   -0.6961    2.1517   3.1416
  5.0    upwind       9.09    0.0000    0.4836   3.1416
  5.0      SUPG       0.00    0.0000    0.5773   3.1416
 
vel = 0 : skew/sym = 0.00e+00,  loop W = -4.29e-16
pi * u  = 3.1416

Al empujar el deslizador Pe_h más allá de 1, la barra a_E cruza en rojo y los nodos se hunden por debajo de cero.

Con beta en cero y Pe_h en 2, el mínimo nodal cae cerca del -0.3334 de la tabla. Al pulsar SUPG, el error se queda en 0 % y la línea rosa antisimétrica no se mueve de 3.1416.

Los valores de contorno están fijados en 0 y 1, y aun así la solución de Galerkin con $Pe_h = 2$ baja hasta $-0.3334$. Con $Pe_h = 5$ llega a $-0.6961$. Ninguno de los dos valores puede existir físicamente.

En $Pe_h = 1$ un coeficiente cambia de signo#

Definiendo el número de Péclet de celda como $Pe_h = u h / (2\epsilon)$, el coeficiente derecho de esa fila se reescribe así.

En cuanto $Pe_h > 1$, el término fuera de la diagonal se vuelve positivo. En ese instante la matriz deja de ser una M-matriz, y con ella se va el principio del máximo discreto: la garantía de que los valores interiores permanecen dentro del rango fijado por los contornos. Por eso min phi en la tabla se mantiene en cero hasta $Pe_h = 1.0$ y se vuelve negativo en 2.0. El instante en que a_E cruzó en rojo es este cambio de signo.

La estabilización no devuelve la simetría#

Aquí aparece una lectura equivocada muy común. Cuando se dice que el upwind o SUPG "recuperan la estabilidad", eso no significa que vuelvan la simetría y el principio de minimización. Basta mirar la columna loop W: tres esquemas y cuatro valores de $Pe_h$, y las doce filas marcan 3.1416. Ese número es $\pi u$, sin rastro de $\epsilon$ ni de $\beta$.

La razón es sencilla. Lo único que $\beta$ agranda es $\epsilon_{\text{eff}}$, que vive en la parte simétrica de la matriz. La parte antisimétrica, $u$, queda intacta. Que la razón skew/sym baje de 2.1517 a 0.4836 con $Pe_h = 5$ no es que la parte antisimétrica encoja: es el denominador simétrico el que crece.

Con $Pe_h = 0.5$, Galerkin imprime max|err|% 3.45, min phi 0.0000 y skew/sym 0.2924. Con $Pe_h = 1.0$ el error sube a 13.53 y el mínimo sigue en 0.0000. De $K_{i,i+1} = (\epsilon/h)(Pe_h - 1)$, el coeficiente derecho es cero en $Pe_h = 1$ y se vuelve positivo en cuanto el Péclet de celda pasa de 1. Con $Pe_h = 2$ el mínimo de Galerkin es -0.3334, el error 35.17 % y skew/sym 1.1010. En la misma fila el upwind tiene mínimo -0.0000, error 18.17 % y skew/sym 0.3885. SUPG tiene error 0.00 y skew/sym 0.5572. Las tres entradas de loop W marcan 3.1416. Con $u = 1$, $\pi u$ es 3.1416, así que el trabajo de la vuelta es el mismo si $\beta$ vale 0, 1 o $\coth(Pe_h) - 1/Pe_h$. Con $Pe_h = 5$ Galerkin llega a mínimo -0.6961, error 69.61 % y skew/sym 2.1517. El upwind baja el error a 9.09 % y skew/sym a 0.4836. Lo que encogió es el denominador simétrico. El control con velocidad nula imprime skew/sym 0.00e+00 y loop W -4.29e-16. Diez elementos y $u = 1$ dan $h = 0.1$, y $Pe_h = uh/(2\epsilon)$ pone $\epsilon = 0.025$ en $Pe_h = 2$. Inflar ese $\epsilon$ con $\beta$ deja $K_{i,i+1} - K_{i+1,i} = u$ igual a 1.

Así que lo que la estabilización compra de verdad es una garantía más débil. Se renuncia al principio de minimización — la mejor aproximación en norma de energía — y a cambio se obtiene la propiedad de M-matriz y el principio del máximo discreto. El pago se hace en precisión. Con $Pe_h = 2$ el upwind eliminó la oscilación, pero arrastra un error del 18.17 %. SUPG, en cambio, añade solo el $\beta$ justo y resulta exacto en los nodos.

Este valor tiende a 0 cuando $Pe_h \to 0$ y a 1 cuando $Pe_h \to \infty$: si domina la difusión, se apaga la estabilización; si domina la advección, se va al upwind completo. Ahora bien, la exactitud nodal es un privilegio del problema unidimensional con coeficientes constantes. En dos dimensiones hace falta la forma original de SUPG, que añade difusión artificial solo en la dirección de la línea de corriente, y ni siquiera entonces queda nada de "exacto".

El nombre que el método de volúmenes finitos usa en el mismo sitio#

El cálculo anterior habló en lenguaje de elementos finitos, pero la conclusión no depende de la discretización. En volúmenes finitos, un término convectivo con diferencias centradas produce exactamente el mismo esténcil y cambia de signo exactamente en el mismo $Pe_h = 1$. Ahí es donde entra el upwind de primer orden, y la difusión numérica que introduce vale lo mismo que el $u h / 2$ correspondiente a $\beta = 1$.

Los dos mundos se separan en la conservación. El upwind de volúmenes finitos modifica el flujo en la cara, así que el balance global se mantiene intacto. Es el mismo tipo de bifurcación que la forma conservativa y la primitiva al discrepar en la velocidad del choque. La difusión artificial de elementos finitos añade un término a la matriz de rigidez, de modo que lo que sigue conservando hay que mirarlo aparte.

Elegir otro miembro de la familia de residuos ponderados cuesta algo distinto. El método de mínimos cuadrados minimiza la norma del residuo, así que siempre produce una matriz simétrica definida positiva: el principio de minimización regresa. A cambio, el número de condición se eleva al cuadrado, y en elementos lineales la segunda derivada se anula dentro de cada elemento, con lo que el término difusivo desaparece por completo. La base y la regla de integración cambian en silencio el carácter de la formulación.

Con $Pe_h = 2$, la cota se compra antes que la malla fina#

La frase de los apuntes sobre que todos los métodos coinciden vale solo sobre operadores autoadjuntos. La difusión, la elasticidad y el flujo potencial viven ahí dentro. La advección trae $K_{i,i+1} - K_{i+1,i} = u$, la simetría se rompe y loop W se queda en 3.1416, o sea $\pi u$, haga lo que haga $\epsilon$. El óptimo en norma de energía que prometía CG ya no está. Una matriz no simétrica pasa al lado de Krylov, de Jacobi a BiCGSTAB.

Un min phi negativo significa que el Péclet de celda pasó de 1. Con $Pe_h = 2$ Galerkin llega a -0.3334 con 35.17 % de error. El upwind borra la oscilación y deja 18.17 % de error. SUPG es exacto en los nodos, pero skew/sym solo baja de 1.1010 a 0.5572. La parte antisimétrica no se movió. En tres dimensiones, bajar $Pe_h$ de 5 a 1 multiplica las celdas por 125. Cuando ese factor no cabe, lo que se compra es la cota, no la precisión.