Cuando las características chocan, nace una onda de choque — Burgers y el plano (x,t)

En 1948 Burgers soltó una ecuación de una línea como modelo de turbulencia. Sin fuerza externa, una curva suave se vuelve acantilado en tiempo finito. El cumpleaños de ese acantilado es el instante en que se cruzan las características (curvas espacio-temporales por donde viaja la información). La solución débil — la forma integral que queda cuando la derivada muere — corta el solape con una sola línea.

Al subir la amplitud $A$ o el número de onda $k$, la línea roja punteada, el tiempo de ruptura, se desplaza hacia arriba.

Esa línea marca la primera vez que una recta de cresta pisa la recta de valle que va delante. Con $A$ mayor el encuentro se adelanta. Con $k = 2$ los cruces llegan de dos en dos.

Dos líneas se encuentran en un mismo punto#

Las ecuaciones hiperbólicas tratan de propagación de señales. Mientras la información fluya por características a velocidades fijas, todo es tranquilo. En la advección lineal $\partial_t q + u\,\partial_x q = 0$ con $u$ constante, las características son rectas paralelas y el perfil inicial se desplaza lateralmente. En esa imagen no cabe ninguna onda de choque.

Todo cambia cuando aparece la no linealidad. Si $u$ depende de $q$ mismo, los valores grandes corren rápido y los pequeños despacio. En el instante en que un valor rápido alcanza a uno lento, dos valores llegan al mismo punto $(x,t)$. La función se vuelve multivaluada, físicamente imposible — ese sinsentido es justo el nacimiento de un choque.

Una ecuación que se transporta a sí misma#

La forma conservativa de la ecuación de Burgers es:

$q$ es el escalar conservado y $f(q) = q^2/2$ es el flujo. La regla de la cadena da la forma no conservativa:

La velocidad de advección es $u(x,t) = q(x,t)$. La cantidad se transporta a sí misma. Esa única no linealidad es el origen de todo el drama. El jacobiano $\partial f/\partial q = q$ es positivo donde la señal va a la derecha y negativo donde va a la izquierda. Un perfil inicial que cambia de signo en el espacio — la onda sinusoidal es el caso típico — obliga a los puntos rápidos a perseguir a los lentos.

El plano (x,t) que dibujan las características#

Una característica es una curva que cumple $dx/dt = q$. Como sobre ella la derivada material $Dq/Dt$ es nula, $q$ queda congelado en su valor inicial a lo largo de la curva. La característica resulta ser una recta en el plano (x,t) con pendiente $1/q_0(x_0)$; el valor inicial $q_0(x_0)$ es la velocidad de esa recta.

Para la condición inicial $q_0(x_0) = 0.55 + A\sin(2\pi k x_0)$:

• Las rectas que parten de los picos ($q_0 > 0.55$) se inclinan fuertemente a la derecha (rápidas).
• Las que parten de los valles ($q_0 < 0.55$) se inclinan poco (lentas).

La línea del pico acabará pisando la punta de la línea del valle inmediatamente delante. El punto donde dos líneas se cruzan es el primer choque.

El tiempo de ruptura t* — el primer cruce#

Dos características vecinas, partiendo de $x_0$ y $x_0 + \Delta x$, se encuentran en el tiempo $t$ cuando:

Tomando el límite $\Delta x \to 0$:

Sólo las pendientes negativas producen un encuentro. El choque más temprano ocurre donde $q_0'(x_0)$ es más negativa:

Cuanto mayores la amplitud o el número de onda de $q_0$, antes se rompe. Para el seno anterior, $q_0'(x_0) = 2\pi k A\cos(2\pi k x_0)$, con mínimo $-2\pi k A$:

Con $A = 0.35$, $k = 1$, sale $t^* \approx 0.455$. Antes la curva se desfasa suavemente; al cruzar $t^*$ se levanta una pared vertical.

La pendiente más negativa es $-2\pi k A$. Como $2\pi \approx 6.28$, $2\pi \times 0.35 \approx 2.20$ y el recíproco es $0.455$. La cresta vale $0.55 + 0.35 = 0.90$ y el valle $0.20$, así que el salto inicial es $0.70$. Coincide con la altura de choque posterior, cerca de $0.7$. Hasta $t = 1.2$ el salto no crece. Solo se redondean las esquinas.

Si $A$ pasa a $0.45$, el denominador sube a unos $2.83$ y $t^*$ baja cerca de $0.35$. Con número de onda $2$ el denominador se duplica, $t^* \approx 0.227$. Misma amplitud, longitud de onda más corta, cruce más temprano.

Solución débil y Rankine–Hugoniot — el cuchillo que corta#

Más allá de $t^*$, la lectura ingenua escupe tres valores de $q$ en el mismo $x$ — imposible. La salida es la solución débil: dejar la forma diferencial y volver a la ley de conservación integral.

Esto tiene sentido aun cuando $q$ no sea derivable — es decir, aunque sea discontinuo. Si dos estados planos $q_L$ y $q_R$ están unidos por un choque que avanza a velocidad $s$:

Es la condición de Rankine–Hugoniot. La velocidad del choque es la media aritmética de los dos estados. Esta condición corta la zona multivaluada con una línea y exige que las áreas a ambos lados sean iguales — la regla de áreas iguales (equal-area rule).

Hace falta además una condición de entropía. El choque sólo es físico si $q_L > q_R$. Si $q_L < q_R$, la respuesta es una rarefacción en abanico (rarefaction). Si se ignora esta condición, queda como respuesta una expansión-choque matemáticamente válida pero no física.

Al meter $0.90$ y $0.20$ en $s = (q_L + q_R)/2$ sale $s = 0.55$, el fondo de la senoide. El choque se desliza a esa velocidad media y deja solo los dos estados. Invertir el orden no cambia la fórmula: sigue dando $0.55$. Cambia la respuesta, de choque a rarefacción. Una media aritmética no distingue las dos ondas.

El paso de malla es $1/256 \approx 0.00391$. Con CFL $0.4$ y $\max|q| = 0.90$ inicial, $\Delta t \approx 0.00174$. Llegar a $t = 1.2$ pide unos $690$ pasos. Tras el choque, un pico local de $\max|q|$ hace que la misma fórmula parta el paso más fino. Ese cuchillo en cada interfaz es el mismo que en un tubo de choque de una celda.

Python: del trazado de características al choque de Godunov#

Trazamos las características y atrapamos el choque débil con un esquema de Godunov de primer orden.

import numpy as np
 
def burgers_breaking_time(q0_prime_min):
    # q0_prime_min: minimo de q0'(x) en x (debe ser negativo para tener choque)
    if q0_prime_min >= 0:
        return float('inf')
    return -1.0 / q0_prime_min
 
def trace_characteristics(q0_func, x0_arr, t_arr):
    # x[i, j] = x0_arr[i] + q0(x0_arr[i]) * t_arr[j]
    return x0_arr[:, None] + q0_func(x0_arr)[:, None] * t_arr[None, :]
 
def rankine_hugoniot_speed(qL, qR):
    # velocidad del choque s = (qL + qR) / 2
    return 0.5 * (qL + qR)
 
def godunov_burgers_step(q, dx, dt):
    # forma conservativa: q^{n+1} = q^n - (dt/dx) * (F_{i+1/2} - F_{i-1/2})
    qL = q                          # estado izquierdo de la interfaz derecha
    qR = np.roll(q, -1)             # estado derecho de la interfaz derecha
    s = 0.5 * (qL + qR)
    f_shock = np.where(s >= 0, 0.5 * qL**2, 0.5 * qR**2)
    f_rare  = np.where(qL >= 0, 0.5 * qL**2,
              np.where(qR <= 0, 0.5 * qR**2, 0.0))
    F_iph = np.where(qL > qR, f_shock, f_rare)   # flujo en la interfaz derecha
    F_imh = np.roll(F_iph, 1)                    # interfaz izquierda = desplazar
    return q - (dt / dx) * (F_iph - F_imh)
 
# Parametros de la simulacion --------------------------------
N, L = 256, 1.0
A, k = 0.35, 1
x = (np.arange(N) + 0.5) * (L / N)
q0 = lambda xx: 0.55 + A * np.sin(2 * np.pi * k * xx)
q  = q0(x)
 
t_star = burgers_breaking_time(-2 * np.pi * k * A)
print(f"breaking time t* = {t_star:.3f}")
 
t_end, t = 1.2, 0.0
while t < t_end:
    dt = 0.4 * (L / N) / np.max(np.abs(q))    # CFL = 0.4
    q  = godunov_burgers_step(q, L / N, dt)
    t += dt
print(f"shock height (max - min) at t = {t:.2f}: {q.max() - q.min():.3f}")

Aparece breaking time t* = 0.455 y una altura de choque cercana a 0.7. El esquema de Godunov corta automáticamente la región multivaluada en una sola línea — el cuchillo es la línea np.where(qL > qR, f_shock, f_rare). Esa línea resuelve el problema de Riemann exacto en cada interfaz y garantiza la solución débil.

Mira chocar las líneas#

En la figura de arriba, la curva cian y el trazo inicial tenue encierran áreas iguales a cada lado del choque. Es la regla de áreas iguales. Con $k = 2$, los dos cruces se funden pronto y el salto crece.

Áreas iguales quiere decir que, tras el corte, $\int q\,dx$ sigue coincidiendo con la diferencia de flujo en los extremos. Si el choque camina a la media de $0.90$ y $0.20$, esto es $0.55$, los dos triángulos recortados se cancelan. La altura se queda en $0.70$. Solo se mueve la posición, a $s = 0.55$.

Dos motivos por los que falla la captura numérica del choque#

(1) Violación de CFL. Si se rompe $\Delta t \le \Delta x / \max|q|$, el dominio numérico adelanta al físico. Si la amplitud parece crecer con el tiempo, el primer sospechoso es $dt$. El código de arriba recalcula $dt$ en cada paso con np.max(np.abs(q)), lo cual importa aún más justo después de formarse el choque, porque $q$ puede picar localmente.

(2) Entropía mal tratada. Sólo con Rankine–Hugoniot caben choques de expansión no físicos. Solvers linealizados sencillos como Roe producen un "entropy glitch" en rarefacciones transónicas (el abanico cruza el cero). Godunov, HLLE y los Riemann solvers exactos imponen la condición de entropía automáticamente; con Roe hace falta una corrección de entropía explícita estilo Harten.

Cuando un código LES compresible muestra posiciones de choque oscilantes o presiones negativas, casi siempre es uno de estos dos.

Si un atasco cuenta como un choque#

Los coches se aprietan en un carril y bajan la velocidad. El de delante va lento, el de detrás rápido. El contacto de parachoques es un cruce en el plano (x,t). Esa cabeza es un choque de Burgers y, mientras no se deshace, deriva aguas arriba a la velocidad de Rankine–Hugoniot.

Dejar el atasco como curva continua conserva el pliegue multivaluado. Cortarlo con una línea o abrirlo como rarefacción lo decide qué estado es mayor. Si un escalón empieza a oscilar, el primer sitio donde mirar es un limitador de flujo.