[Reseña de artículo] Los dos fragmentos no se repelieron: volaron pegados — el ángulo de equilibrio de 144° de un par de esferas a Mach 20

Supongamos que un meteorito se parte en dos dentro de la atmósfera. ¿Los dos trozos se alejan uno del otro? Si la respuesta fuera "siempre", reconstruir la entrada a partir de los campos de cráteres en tierra sería mucho más fácil. El cálculo a Mach 20 que Whalen, Deiterding y Laurence publicaron en JFM 2026 (vol. 1029, A35) sacude esa respuesta. Para ciertos ángulos de alineación las dos esferas vuelan pegadas y el par entero llega incluso a generar sustentación. Este texto reproduce de dónde sale ese ángulo con una aproximación newtoniana de 30 líneas. El valor obtenido fue 143.6°, dentro del intervalo estable de 132°–145.7° del artículo.

Q1. Por qué sigue abierto el problema de solo dos fragmentos#

El cálculo de fragmentos meteoríticos tiene dos ramas antiguas. Con unos pocos trozos se usa el enfoque discrete-fragment (fragmento discreto), que sigue cada trozo por separado. Con trozos prácticamente infinitos se usa el enfoque debris-cloud (nube de escombros), que aplasta y dispersa el conjunto como si fuera un líquido. El modelo de dos cuerpos que Passey y Melosh propusieron en 1980 es el prototipo del primero. Si se supone que las dos esferas solo se empujan lateralmente, la velocidad transversal final se resume en esta proporcionalidad.

$V_T \propto V\sqrt{\rho_a/\rho_m}$

$V$ es la velocidad de entrada, $\rho_a$ la densidad atmosférica y $\rho_m$ la del meteorito. Solo queda una constante de proporcionalidad, y los valores retrodeducidos de campos de cráteres se dispersan entre 0.03 y 2.28. Dos órdenes de magnitud.

El problema está en medio. La zona de poblaciones intermedias —ni 2 trozos ni varios miles— estaba casi vacía. El artículo apunta de frente a ese tramo. Con esferas iguales ordenadas en arreglos de 2, 4 y 13, cambiando solo la actitud inicial, se calcularon 83 casos. El flujo lo resuelve AMROC, que integra las ecuaciones de Euler con embedded boundary, y las esferas las resuelve DYNA3D incluyendo el contacto. El gas es perfecto con $\gamma=1.4$, la corriente entrante es Mach 20 y la relación de densidades entre esfera y gas es $10^4$. La malla se refina sola siguiendo el gradiente de densidad. La estructura de costo de este esquema es la misma que se trató en etiquetado AMR y refluxing.

Conviene manipular directamente la simulación de abajo.

Basta con dejar el ángulo de alineación $\theta_0$ en cualquier punto y pulsar release. Si la flecha curva sale verde, el par se abre; si sale roja, se cierra. Todos los ángulos iniciales terminan absorbidos por el único punto donde la curva de la derecha cruza el cero.

Q2. Por qué un par pegado se va justo a un solo ángulo#

El ángulo de alineación $\theta$ se define entre la dirección de la corriente libre y el segmento "esfera secundaria → esfera primaria". Con $\theta = 180^\circ$ la esfera trasera se asienta exactamente en la estela de la delantera. Con $\theta = 90^\circ$ quedan lado a lado.

Mientras las dos esferas mantienen el contacto se las puede tratar como un solo cuerpo rígido. Ahí la fuerza de contacto es interna y desaparece del momento del par completo. Solo queda el momento aerodinámico.

$S_i$ es la superficie de cada esfera, $\hat{n}$ la normal exterior y $\mathbf{x}_{\rm com}$ el centro de masa del par. La condición de equilibrio solo exige que esto valga 0; la estabilidad la decide el signo en su entorno.

Empecemos por $\theta = 180^\circ$. La esfera de aguas arriba recibe la corriente de frente y la de aguas abajo queda completamente tapada, con resistencia casi nula. La diferencia de resistencia entre ambas actúa en una posición desplazada del centro de masa, así que cualquier giro pequeño genera un momento que lo amplifica. Es un equilibrio inestable. El artículo también anota que a 180° no se observó cambio de actitud, y añade que debería ser inestable porque la resistencia de la esfera de aguas arriba es mayor.

En el otro extremo, $\theta = 90^\circ$ es simétrico y el momento vale exactamente 0. Entonces en algún punto intermedio tiene que invertirse el signo. Ese punto es el ángulo de equilibrio estable.

Q3. ¿Sale ese ángulo solo con la aproximación newtoniana?#

En régimen hipersónico la distribución de presión se puede aproximar con la teoría de impacto newtoniana. Si se supone que la partícula fluida golpea la superficie y pierde todo su momento normal, el coeficiente de presión queda así.

$C_p = C_{p,\max}\,(\hat{u}\cdot\hat{n})^2, \qquad C_{p,\max} = 1.8394$

$\hat{u}$ es el vector unitario de la corriente libre y $\hat{n}$ la normal exterior de la superficie; solo vale en las caras a barlovento, donde $\hat{u}\cdot\hat{n}<0$. $C_{p,\max}$ es el valor de punto de remanso para $\gamma=1.4$ y $M\to\infty$. A esto se le suma una única línea de prueba de ocultamiento: si una cara queda escondida detrás de la otra esfera, su presión se pone a 0.

La superficie esférica se dividió en 24 mil paneles con una retícula de Fibonacci y se barrió de 90° a 180° en pasos de 1°.

import math
 
CP_MAX = 1.8394          # teoría newtoniana modificada, gamma=1.4, M -> infinito
R = 1.0                  # radio de la esfera
N_PANEL = 24000
 
 
def fib_sphere(n):
    """Puntos casi uniformes sobre la esfera unitaria + área de un panel."""
    pts, ga = [], math.pi * (3.0 - math.sqrt(5.0))
    for i in range(n):
        z = 1.0 - 2.0 * (i + 0.5) / n
        rho = math.sqrt(max(0.0, 1.0 - z * z))
        a = ga * i
        pts.append((rho * math.cos(a), rho * math.sin(a), z))
    return pts, 4.0 * math.pi / n
 
 
PANELS, DA = fib_sphere(N_PANEL)
 
 
def pair_geometry(theta_deg):
    """La esfera secundaria (aguas abajo) en el origen; la primaria a 2R según n_hat.
    theta = ángulo entre el x_hat de la corriente libre y el segmento 'secundaria -> primaria'."""
    t = math.radians(theta_deg)
    n_hat = (math.cos(t), math.sin(t), 0.0)
    return (0.0, 0.0, 0.0), tuple(2.0 * R * c for c in n_hat), n_hat
 
 
def shadowed(px, py, pz, cx, cy, cz):
    """¿Es un punto tapado por la esfera en c, que no ve la corriente libre?
    El flujo va en +x, así que se remonta en -x con una prueba cilíndrica."""
    if cx >= px:
        return False
    return (py - cy) ** 2 + (pz - cz) ** 2 < R * R
 
 
def newtonian_cp(nx):
    """nx = componente x de la normal exterior. Si nx < 0, cara a barlovento."""
    return CP_MAX * nx * nx if nx < 0.0 else 0.0
 
 
def pair_loads(theta_deg):
    """Coeficientes de fuerza de cada esfera y momento respecto al centro de masa del par.
    Referencias: la fuerza con q_inf * pi * R^2; el momento, eso mismo por 2R."""
    cs, cp, _ = pair_geometry(theta_deg)
    com = tuple(0.5 * (a + b) for a, b in zip(cs, cp))
    out = []
    for me, other in ((cs, cp), (cp, cs)):
        fx = fy = mz = 0.0
        for ux, uy, uz in PANELS:
            cpv = newtonian_cp(ux)
            if cpv == 0.0:
                continue
            x, y, z = me[0] + R * ux, me[1] + R * uy, me[2] + R * uz
            if shadowed(x, y, z, other[0], other[1], other[2]):
                continue
            dfx, dfy = -cpv * ux * DA * R * R, -cpv * uy * DA * R * R
            fx += dfx
            fy += dfy
            mz += (x - com[0]) * dfy - (y - com[1]) * dfx
        s = math.pi * R * R
        out.append((fx / s, fy / s, mz / (s * 2.0 * R)))
    return out
 
 
def sweep_alignment(lo, hi, step):
    rows, t = [], lo
    while t <= hi + 1e-9:
        (dxs, dys, ms), (dxp, dyp, mp) = pair_loads(t)
        cd, cl, cm = dxs + dxp, dys + dyp, ms + mp
        rows.append((t, cd, cl, cm, dxp, dxs))
        t += step
    return rows
 
 
SLOPE_FLOOR = 5.0e-5     # ruido de integración. Un cruce real es mucho más pronunciado
 
 
def stable_window(rows):
    """Conserva solo los ceros de C_M cuya pendiente supera el ruido."""
    hits = []
    for (t0, _, _, m0, _, _), (t1, _, _, m1, _, _) in zip(rows, rows[1:]):
        if m0 * m1 >= 0.0 or abs(m1 - m0) < SLOPE_FLOOR * (t1 - t0):
            continue
        te = t0 + (t1 - t0) * (-m0) / (m1 - m0)
        hits.append((te, "stable" if m1 < m0 else "unstable", m1 - m0))
    return hits
 
 
rows = sweep_alignment(90.0, 180.0, 1.0)
 
print("theta   C_D    C_L     C_M      L/D    C_D,up  C_D,down")
for t, cd, cl, cm, dxp, dxs in rows:
    if abs(t % 7.5) < 1e-9:
        print(f"{t:5.1f} {cd:6.3f} {cl:7.3f} {cm:8.4f} {cl / cd:7.3f} "
              f"{dxp:7.3f} {dxs:8.3f}")
 
print()
for te, kind, slope in stable_window(rows):
    print(f"C_M = 0 at theta = {te:6.2f} deg  slope {slope:+.2e}/deg  -> {kind}")
 
eq = [h[0] for h in stable_window(rows) if h[1] == "stable"][0]
near = min(rows, key=lambda r: abs(r[0] - eq))
print(f"\nat the stable angle: L/D = {near[2] / near[1]:.3f}, C_D = {near[1]:.3f}")
lo = min(r[0] for r in rows if r[3] > SLOPE_FLOOR)
print(f"restoring sign holds from theta = {lo:.1f} deg up to the equilibrium")
print(f"tandem (180 deg): C_D,up = {rows[-1][4]:.3f}, C_D,down = {rows[-1][5]:.3f}")

theta   C_D    C_L     C_M      L/D    C_D,up  C_D,down
 90.0  1.839   0.000   0.0000   0.000   0.920    0.920
105.0  1.839   0.003   0.0000   0.002   0.920    0.919
120.0  1.819   0.037   0.0003   0.020   0.920    0.899
135.0  1.714   0.128   0.0007   0.075   0.920    0.794
150.0  1.459   0.215  -0.0020   0.147   0.920    0.539
165.0  1.112   0.170  -0.0122   0.152   0.920    0.193
180.0  0.920   0.000  -0.0000   0.000   0.920    0.000
 
C_M = 0 at theta = 143.61 deg  slope -2.09e-04/deg  -> stable
 
at the stable angle: L/D = 0.119, C_D = 1.580
restoring sign holds from theta = 111.0 deg up to the equilibrium
tandem (180 deg): C_D,up = 0.920, C_D,down = 0.000

Salió 143.61°. El intervalo estable que el artículo encontró rastreando historias de coeficientes y momentos es 132°–145.7°, y el máximo de eficiencia aerodinámica del par en contacto está en 141.5°. Treinta líneas sin malla y sin ondas de choque cayeron dentro.

La eficiencia aerodinámica es 0.119, por debajo del promedio 0.197 y del máximo 0.22 del intervalo estable del artículo. La razón de la diferencia es clara: la aproximación newtoniana no sabe nada de la interferencia entre ondas de choque. La mancha de alta presión que aparece en las caras interiores entre ambas esferas, visible en las figuras del artículo, nace del encuentro de las dos ondas de proa, y en este modelo ese lugar sencillamente no existe. Como se vio en el tren de choques oblicuos, la presión donde se superponen ondas de choque supera la suma simple. Así que lo que este código acierta no es la magnitud, sino el signo y la posición.

También vale la pena leer dónde cambia el signo en la tabla. $C_L$ alcanza su máximo de 0.215 a 150°, pero $C_M$ pasa por cero antes, a 143.6°. La actitud de mayor sustentación y la actitud que se sostiene sola no coinciden.

Q4. De dónde sale el empuje lateral una vez separadas#

Cuando el contacto se rompe la historia cambia. Si las dos esferas estaban lado a lado desde el inicio (cerca de 90°), una onda de proa común presiona las caras interiores y las separa. Esa repulsión mutua desaparece rápido al abrirse el hueco, y la velocidad transversal adimensional del artículo se detiene alrededor de 0.2.

Pero cuando la esfera trasera se monta sobre la onda de choque de la delantera ocurre otra cosa. Justo detrás del choque la presión es alta y por fuera está la corriente libre. Si la esfera avanza mordiendo esa frontera por la mitad, recibe de forma continua un empuje hacia afuera. Es el fenómeno que Laurence y Deiterding bautizaron en 2011 como shock surfing (surfeo de onda de choque). En el artículo, este efecto sube la velocidad transversal final de 0.2 a 0.25 y alarga el tiempo hasta la separación de $1.7\tau_s$ a $4\tau_s$.

Las referencias de tiempo y velocidad son estas.

$\tau_s = \sqrt{\frac{\rho_{\rm sph}}{\rho_\infty}}\,\frac{r_c}{u_\infty}$

$r_c$ es el radio circunscrito del cúmulo. Como la relación de densidades es $10^4$, $\tau_s$ vale 100 veces el tiempo de paso. Adimensionalizando con esta escala, la ecuación de movimiento relativo se queda con un solo coeficiente.

$\frac{\mathrm{d}V'}{\mathrm{d}t'} = \frac{3}{8}\bigl(C_{F,\rm trail} - C_{F,\rm lead}\bigr)$

El $3/8$ sale de la relación entre el volumen y el área frontal de la esfera, y la relación de densidades quedó absorbida por completo dentro de la escala. La animación de abajo integra esta ecuación tal cual.

Vale la pena arrastrar el release angle hacia arriba desde 90°. Cerca de 90° la línea continua y la punteada se superponen: la onda común hace todo el trabajo y el surfeo no aporta nada. Pasados los 135° ambas curvas se separan, y el borde de la esfera se ilumina solo mientras muerde la línea rosa del choque.

El comportamiento durante el contacto también está en el modelo. Al soltar a 165°, las dos esferas ruedan pegadas durante 3–4$\tau_s$ y se separan recién cuando el ángulo de alineación baja a unos 120°. Coincide a grandes rasgos con los 4–6$\tau_s$ de contacto y el ángulo de desprendimiento cercano a 130° que reporta el artículo. La forma en que la resistencia se desploma en la estela es otro asunto distinto de la curva de resistencia de una esfera aislada que se vio en la crisis de resistencia. Aquí la causa no es el punto de desprendimiento, sino el ocultamiento.

Q5. Entonces, ¿por qué a 180° no pasa nada?#

Un tándem perfectamente alineado es simétrico. La fuerza transversal es 0 y el momento también. El cálculo termina con las dos esferas pegadas durante 8$\tau_s$ sin cambiar la actitud en absoluto. La observación del artículo es la misma.

Pero este es el equilibrio de un lápiz puesto de punta. Como la resistencia de la esfera de aguas arriba supera la de aguas abajo (0.920 contra 0.000 en la tabla), la menor inclinación trae consigo un momento que la amplifica. En el cálculo real, incluso soltando a 172.5° casi no hay movimiento al principio, luego aparece un giro "rodante" cada vez más rápido y el contacto acaba rompiéndose.

Aquí surge un valor de importancia práctica. En un par que mantuvo el contacto y generó sustentación, el propio centro de masa se desplaza lateralmente. El artículo anota que la velocidad transversal del centro de masa de los pares soltados a 150°, 157.5° y 172.5° fue de 0.42, 0.32 y 0.28 respectivamente, mayor que la velocidad relativa entre fragmentos. Contar solo la repulsión de cada fragmento por separado deja fuera esta componente entera.

Lo que queda al llegar a 13#

El artículo continúa el mismo estudio con 38 arreglos tetraédricos de 4 esferas y 34 arreglos cúbicos centrados en las caras de 13 esferas. La velocidad transversal final de cada esfera se ordena bastante bien con un solo ángulo polar inicial, y esa tendencia aparece de forma parecida con 4 y con 13. Eso sí, la influencia del achatamiento global del cúmulo sobre el comportamiento en volumen se diluye a medida que crece la población. Cuantos más trozos, más homogénea es la separación, lo cual es una pista sobre a partir de dónde el enfoque debris-cloud empieza a servir.

En resumen, lo que gobierna la separación de un par de fragmentos no es la magnitud de la repulsión. Es cuándo se rompe el contacto, y eso lo decide el signo del momento. Y ese signo se puede conocer sin malla alguna. Con la presión superficial puesta en aproximación newtoniana y el ocultamiento bien contado, la respuesta de que las dos esferas encuentran solas un ángulo de 143.6° sale en 30 líneas.