El flujo supersónico no sabe lo que ocurre aguas arriba — curvas características y ondas sonoras en las ecuaciones de Euler

La cuña del dominio de influencia cruza la vertical en Mach 1. Por encima, una salida supersónica ya no tiene una característica entrante que fijar.

Un caza supersónico llega antes que el sonido de su propio motor. Primero pasa por encima y el estruendo viene después. La razón está en las velocidades características de las ecuaciones de Euler. Una perturbación pequeña se separa en tres señales, y sus trazas en el plano $(x,t)$ son las curvas características (el rastro que deja cada señal).

Con el control M en 0.4, la onda azul hacia la izquierda alcanza aguas arriba de la fuente. Al pasarlo de 1, esa traza desaparece de la figura.

background flow u₀ (Mach, c₀=1) = 0.40

The blue tint marks the upstream region the left-going sound has reached. Push u₀ past 1 and the blue tint vanishes: the C− signal now has a positive speed, so every signal runs downstream and the source can no longer talk to what lies upstream.

Ahí cambia de signo la pendiente de la línea azul. Conviene mirar si las tres señales corren solo hacia la derecha y si el lado izquierdo queda vacío.

Las ecuaciones compresibles son propagación de señales#

El flujo compresible no viscoso en una dimensión se escribe con tres leyes de conservación.

\partial_t \rho + \partial_x(\rho u) = 0

\partial_t(\rho u) + \partial_x(\rho u^2 + P) = 0

\partial_t(\rho e_{\text{tot}}) + \partial_x\big[(\rho e_{\text{tot}} + P)\,u\big] = 0

$\rho$ es la densidad, $u$ la velocidad, $P$ la presión y $e_{\text{tot}}$ la energía específica total. A primera vista parecen ecuaciones de un fluido que "fluye". Pero la esencia de estas ecuaciones no es el transporte, sino la propagación de señales. Toda la dinámica del flujo se descompone en unas pocas señales que viajan a velocidades determinadas. Para ver cuáles son esas velocidades basta con linealizar la ecuación por un momento.

Una pequeña perturbación genera una ecuación de onda#

Sobre un estado de fondo $(\rho_0, u_0, P_0)$ se superpone una perturbación muy pequeña. Tomando $\rho = \rho_0 + \rho_1$ , $u = u_0 + u_1$ y usando $P = K\rho^\gamma$ , se conservan solo los términos de primer orden. Definiendo la derivada que se mueve con el fondo $D_t \equiv \partial_t + u_0\,\partial_x$ (derivada co-móvil), al combinar las dos ecuaciones de perturbación resulta una única ecuación de onda.

D_t^2\,\rho_1 - c_0^2\,\partial_x^2\,\rho_1 = 0

Aquí $c_0$ es la velocidad del sonido.

c_0 = \sqrt{\gamma\,\frac{P_0}{\rho_0}}

$\gamma$ es la razón de calores específicos. Al insertar $\rho_1 \sim e^{i(kx-\omega t)}$ en esta ecuación se obtiene la relación de dispersión.

\frac{\omega}{k} = u_0 \pm c_0

$\omega$ es la frecuencia angular y $k$ el número de onda. La perturbación se deja arrastrar por el flujo de fondo $u_0$ y se propaga hacia ambos lados a la velocidad del sonido $c_0$ . Es decir, la perturbación en un punto se separa en dos velocidades: $u_0 + c_0$ y $u_0 - c_0$ .

Las curvas características que se separan en tres ramas#

Las dos ondas sonoras no son todo. La tercera señal es el propio fluido. Supongamos que se pinta el fluido con color. Si en $t=0$ se tiñe la izquierda de azul y la derecha de rojo, la frontera de color $\phi$ simplemente se deja arrastrar por el flujo.

\partial_t \phi + u_0\,\partial_x \phi = 0

Esta señal se desplaza a la velocidad $u_0$ . Esta es la característica de entropía (o superficie de contacto). En resumen, el flujo unidimensional tiene tres velocidades características.

\lambda_- = u_0 - c_0, \qquad \lambda_0 = u_0, \qquad \lambda_+ = u_0 + c_0

En el plano $(x,t)$ , las líneas que satisfacen $\dfrac{dx}{dt} = \lambda$ se llaman curvas características. $\lambda_\pm$ son las ondas sonoras y $\lambda_0$ es el transporte de materia. A lo largo de cada característica acústica se conservan constantes los invariantes de Riemann $u \pm \dfrac{2c}{\gamma-1}$ , y a lo largo de la característica material se conserva la entropía.

Con $M = u_0/c_0$ en 0.4, $\lambda_-$ sigue negativa y la onda sonora izquierda todavía llega aguas arriba. Cuando $M$ supera 1, $\lambda_- = u_0 - c_0$ se vuelve positiva y las tres señales viajan solo aguas abajo.

La misma descomposición, vista en la estructura propia de Euler para una normal arbitraria, añade una raíz repetida y componentes tangenciales. En una dimensión esa repetición se pliega en la única onda de entropía.

En régimen supersónico el cono de señal se vuelca#

La misma historia se ve con más claridad en un diagrama espacio-temporal. La cuña entre las tres curvas características $C_-,\,C_0,\,C_+$ que salen de la fuente es el dominio de influencia (domain of influence). Es el conjunto de eventos espacio-temporales a los que la fuente puede enviar una señal.

Si $M<1$ , esta cuña se abre hacia ambos lados de la línea vertical. La señal va tanto aguas arriba como aguas abajo. Si $M>1$ , la cuña entera se vuelca hacia la derecha. Incluso la característica $C_-$ se inclina hacia la derecha, y ningún evento aguas arriba queda dentro del dominio de influencia. Esta es la razón de que en el flujo supersónico "aguas arriba no sepa lo que ocurre aguas abajo". También es el fundamento numérico de por qué no se deben imponer condiciones de contorno en una salida supersónica.

Se puede arrastrar la sonda redonda en el diagrama de abajo.

background flow u₀ (Mach, c₀=1) = 0.40

Drag the probe across the plane. At M = 0.4 the wedge straddles the vertical, so events both upstream and downstream are reachable. Slide u₀ above 1 and the whole wedge leans right — drop the probe anywhere to the left of the source and it stays red.

Si la sonda se coloca dentro de la cuña se vuelve verde (alcanzable) y, si se coloca fuera, roja (inalcanzable). Al subir $u_0$ por encima de 1 y colocar la sonda en cualquier punto a la izquierda de la fuente, permanece roja. Aguas arriba queda completamente desconectado de la señal.

Verificación de las velocidades características con Python#

Las velocidades características se comparan con $u \pm c$ y $u$ . Al escribir el Euler unidimensional con las variables primitivas $W=(\rho,u,P)$ se obtiene la forma cuasilineal $\partial_t W + A(W)\,\partial_x W = 0$ , y los autovalores de la matriz de coeficientes $A$ son precisamente las velocidades características.

import numpy as np
 
def characteristic_speeds(rho, u, p, gamma=1.4):
    """Jacobiano de variables primitivas del Euler 1D A(W), W=(rho,u,p).
    Autovalores de A = las tres velocidades características."""
    c = np.sqrt(gamma * p / rho)              # velocidad del sonido
    A = np.array([[u,      rho,        0.0  ],
                  [0.0,    u,          1/rho],
                  [0.0,    rho*c**2,   u    ]])
    lam = np.sort(np.linalg.eigvals(A).real)
    return c, lam
 
rho, u, p = 1.225, 220.0, 1.0e5               # aire, u = 220 m/s
c, lam = characteristic_speeds(rho, u, p)
print(f"c0         = {c:8.2f} m/s")
print(f"eig(A)     = {np.round(lam, 2)}")
print(f"u-c,u,u+c  = {np.round([u-c, u, u+c], 2)}")
print(f"Mach       = {u/c:5.3f}  -> {'supersonic' if u > c else 'subsonic'}")

La salida es la siguiente.

c0         =   338.06 m/s
eig(A)     = [-118.06  220.    558.06]
u-c,u,u+c  = [-118.06  220.    558.06]
Mach       =  0.651  -> subsonic

Los autovalores coinciden exactamente con $u-c,\,u,\,u+c$ . La perturbación se avanza en el tiempo y se mide la velocidad del frente de la señal. Los invariantes de Riemann acústicos $w_\pm = u_1 \pm (c_0/\rho_0)\,\rho_1$ son escalares que se dejan arrastrar respectivamente a $u_0 \pm c_0$ , así que basta con transportarlos por separado con esquema upwind.

N, L = 400, 1.0
x = np.linspace(0, L, N, endpoint=False)
dx = L / N
rho0, c0 = 1.0, 1.0
u0 = 0.6 * c0                                  # Mach de fondo 0.6
 
bump = np.exp(-((x - 0.5) ** 2) / (2 * 0.01 ** 2))
wp = bump.copy()                              # invariante derecho (u0 + c0)
wm = bump.copy()                              # invariante izquierdo (u0 - c0)
 
def upwind_step(w, a, dt):
    if a > 0:
        dwdx = (w - np.roll(w, 1)) / dx
    else:
        dwdx = (np.roll(w, -1) - w) / dx
    return w - a * dt * dwdx
 
CFL = 0.4
dt = CFL * dx / (abs(u0) + c0)
steps = int(0.25 / dt)
for _ in range(steps):                        # bucle de avance temporal
    wp = upwind_step(wp, u0 + c0, dt)
    wm = upwind_step(wm, u0 - c0, dt)
 
T = steps * dt
def front_speed(w):
    return (x[np.argmax(w)] - 0.5) / T
print(f"u0+c0: target {u0 + c0:+.3f}, measured {front_speed(wp):+.3f}")
print(f"u0-c0: target {u0 - c0:+.3f}, measured {front_speed(wm):+.3f}")

u0+c0: target +1.600, measured +1.600
u0-c0: target -0.400, measured -0.400

Cada uno de los dos frentes se mueve a $u_0+c_0$ y $u_0-c_0$ . El invariante derecho va aguas abajo y el izquierdo aguas arriba. Al subir $u_0$ por encima de $c_0$ , el valor medido de la segunda línea también se vuelve positivo.

Con los números del aire el hueco se ve mejor. Con $\rho=1.225$ , $u=220$ m/s, $P=10^5$ y $\gamma=1.4$ , la velocidad del sonido es 338.06 m/s. Los autovalores impresos son $-118.06$ , $220$ y $558.06$ , idénticos a $u-c$ , $u$ y $u+c$ . El tramo $558.06-(-118.06)$ es 676.12 m/s, el doble de 338.06. El valor central, 220, solo transporta materia y no sabe de la velocidad del sonido. Mach 0.651 está por debajo de 1, así que $\lambda_-$ sigue negativa y todavía queda una señal que viaja aguas arriba.

El segundo cálculo fija la velocidad del sonido en 1. Con $u_0=0.6$ los objetivos son $+1.600$ y $-0.400$ , y los frentes medidos dan esos números. $-0.400$ es $0.6-1$ . El signo negativo es el alcance aguas arriba. Cuando $u_0$ supera a $c_0$ , esa línea se vuelve positiva, el mismo estado en que la cuña espacio-temporal se vuelca a la derecha.

En el invariante de Riemann $u \pm \dfrac{2c}{\gamma-1}$ , con $\gamma=1.4$ se tiene $\gamma-1=0.4$ y $2/0.4=5$ . Lo que permanece constante a lo largo de una característica acústica es $u\pm 5c$ . Los escalares que transporta el código, $w_\pm = u_1 \pm (c_0/\rho_0)\,\rho_1$ , son ese invariante reducido a una perturbación pequeña. Con $\rho_0=c_0=1$ el coeficiente queda en 1.

Fijar la presión en una salida supersónica devuelve la onda#

El instante en que la cuña del dominio de influencia cruza la vertical es $M=1$ . Por encima no hay ninguna característica que vaya aguas arriba. Si la salida es supersónica, los tres signos apuntan hacia afuera y no queda ningún valor por imponer. Fijar la presión inventa una característica entrante que no existía, y una onda reflejada regresa al dominio.

Una salida subsónica es distinta. Solo $\lambda_-$ apunta hacia adentro, así que basta una condición. Esa característica entrante es la que NSCBC multiplica por un factor de reflexión σ en la salida. Imponer la presión antes de contar los signos mezcla una señal de aguas arriba en el tratamiento de salida, incluso en el aire a Mach 0.651.