[Revisión de artículo] El tetraedro que adelantó una esfera no se dispersó — un índice de achatamiento para cúmulos de fragmentos a Mach 20

La velocidad lateral de un fragmento la fija casi por completo su ángulo polar inicial en el cúmulo, y la del cúmulo entero, el achatamiento B que pliega esos ángulos en un solo número.

El tetraedro que adelantó una esfera no se dispersó#

Cuatro esferas de hierro iguales forman un tetraedro dentro de un flujo a Mach 20. Un vértice apunta directo aguas arriba y se las suelta. La intuición dice que el cúmulo debería estallar. En las simulaciones de Whalen, Deiterding y Laurence (JFM 2026), apenas se separó. La presión que empujaba hacia fuera a las tres esferas traseras quedó compensada exactamente por las fuerzas de contacto.

Este artículo recorre la sección del tetraedro de cuatro esferas (§4). ¿Qué fija la velocidad lateral de un fragmento? ¿Qué número resume un cúmulo entero? Se recalculan en Python el índice de achatamiento $B$ y el índice de asimetría $A$ del artículo. Luego se comprueban su rango y su correlación sobre 40 000 orientaciones aleatorias.

El caso de dos esferas que vuelan pegadas cerca de 144° generando sustentación se trató en la revisión anterior. Hoy el par se convierte en un grupo de cuatro.

Los deslizadores de pitch y yaw mueven el círculo blanco en el plano $A$ – $B$ , y el roll deja ese círculo en su sitio.

pitch 36°yaw 18°roll 0°

B = 0.859 · A = 0.347 · blunt (face into the flow)
θ0 = [39.7, 88.0, 149.2, 86.9]° · V'_i (sketch) = [0.44, 0.70, 0.45, 0.70] · t' = 0.0
color = sketch V'_T (blue slow → red fast). Roll leaves B, A and θ0 unchanged.

vertex-first parte los ángulos polares en 0° y 109.5°. face-first junta tres de ellos en la meseta cercana a 70°. Conviene mirar si el círculo recorre el tramo entre $B=1/3$ y $B=1$ .

1980, Passey–Melosh: dos fragmentos, una constante#

La modelación de fragmentos de meteoroides parte del modelo binario de Passey y Melosh (1980). Supusieron dos esferas empujadas solo de forma lateral y obtuvieron una fórmula para la velocidad lateral final. Escala como $V_T \propto V\sqrt{\rho_a/\rho_m}$ . La constante adimensional restante $C$ variaba mucho entre observaciones.

La otra rama es el modelo pancake (una masa sin resistencia que se aplana y se extiende). Supone que los fragmentos se expanden de forma puramente radial. Casi ningún cálculo había puesto a prueba esa hipótesis. El artículo plantea ambas preguntas. ¿Qué fija la velocidad de un fragmento? ¿Su dirección es de verdad radial?

Escalas de tiempo y velocidad#

El artículo adimensionaliza con dos escalas.

\tau_s = \sqrt{\frac{\rho_{sph}}{\rho_\infty}}\,\frac{r_c}{u_\infty}, \qquad v_s = \sqrt{\frac{\rho_\infty}{\rho_{sph}}}\,u_\infty

$\rho_{sph}$ es la densidad de la esfera, $\rho_\infty$ y $u_\infty$ la densidad y la velocidad de la corriente libre, y $r_c$ el radio circunscrito del cúmulo. La razón de densidades es $10^4$ , así que $v_s = u_\infty/100$ . Para un meteoroide a 20 km/s, $V_T' = 0.35$ equivale a 70 m/s laterales.

La velocidad lateral $V_{T,i}$ es la velocidad relativa al centro de masa, proyectada sobre la dirección radial en el plano $y$ – $z$ . La velocidad colectiva de separación $V_T$ es el promedio de las cuatro esferas. AMROC (un solver de Euler con malla adaptativa) resuelve el flujo. DYNA3D (un código explícito de elementos finitos) resuelve el contacto y las colisiones elásticas.

B y A en lugar de pitch y yaw#

Dos ángulos, pitch y yaw, fijan la orientación del tetraedro. El artículo sostiene que no tienen significado físico y define dos índices nuevos.

B = \max_i\left(\hat{n}_i\cdot\langle -1,0,0\rangle\right)

$\hat{n}_i$ es la normal unitaria exterior de una cara que pasa por tres centros de esfera. El flujo va en $+x$ . Con un vértice al frente, $B=1/3$ . Con una cara de frente al flujo, $B=1$ . Un $B$ mayor significa que más esferas comparten el choque de proa (bow shock).

A = \frac{\left\|\sum_i w_i'\,(\mathbf{x}_{i,l}-\mathbf{x}_{com,l})\right\|}{r_v},\qquad w_i = \begin{cases}(x_i-x_{com})^2 & x_i<x_{com}\\ 0 & x_i\ge x_{com}\end{cases}

$l$ indica la componente lateral, $w_i'$ los pesos normalizados para sumar uno y $r_v$ el radio de la esfera virtual circunscrita. $A$ crece cuanto más se alejan del eje las esferas de aguas arriba.

Un ángulo polar decide un fragmento#

El artículo agrupa cada esfera según su ángulo polar inicial $\theta_0$ (el ángulo respecto al eje aguas arriba, visto desde el centro de masa) en intervalos de 15°. El resultado es claro. La velocidad lateral final sube casi linealmente de 0° a unos 70°. Entre 70° y 90° forma una meseta breve cerca de 0.7. Esa franja es donde el shock surfing (una esfera que se desliza a lo largo del choque de otra) es más probable. No dejar que ese choque se emborrone en varias celdas es el problema de reconstrucción de TENO de Fu (2019).

Pasados los 90°, más esferas quedan arrastradas en la estela de una esfera delantera y la velocidad cae. Por encima de 120° vuelve a subir. Son esferas ocultas detrás de un cuerpo aguas arriba hasta que empieza una interacción entre dos esferas. Cerca de 180°, la simetría mantiene las velocidades bajas.

También se revisó la dirección. El 75 % de las esferas se mantuvo a menos de 20° de su azimut inicial. Para fragmentos de igual tamaño, la hipótesis de expansión radial del modelo pancake se sostiene en buena medida.

Código: 40 000 orientaciones del tetraedro#

$B$ , $A$ y los ángulos polares dependen solo de la orientación. Se pueden contar sin resolver ningún flujo. El código siguiente recorre 40 000 rotaciones aleatorias uniformes y una malla pitch–yaw de 36 puntos como la del artículo. La velocidad de cada esfera sale de la curva aproximada CURVE. Está transcrita de la descripción del texto, no de los datos de la Fig. 16(a).

import math, random
 
R_SPH = 1.0
 
def sub(a, b): return [a[i] - b[i] for i in range(3)]
def dot(a, b): return sum(a[i] * b[i] for i in range(3))
def cross(a, b): return [a[1]*b[2]-a[2]*b[1], a[2]*b[0]-a[0]*b[2], a[0]*b[1]-a[1]*b[0]]
def norm(a): return math.sqrt(dot(a, a))
def rotate(R, p): return [dot(R[i], p) for i in range(3)]
 
def tetra_centres():
    rc = 2 * R_SPH * math.sqrt(6) / 4      # cuatro esferas en contacto: circunradio del tetraedro de centros
    s = rc * math.sqrt(8) / 3              # radio del círculo de la base
    X = [[-rc, 0.0, 0.0]]                  # un vértice apunta aguas arriba (-x)
    for t in (0, 2*math.pi/3, 4*math.pi/3):
        X.append([rc/3, s*math.cos(t), s*math.sin(t)])
    return X, rc
 
def attitude(pitch, yaw):
    p, y = math.radians(pitch), math.radians(yaw)
    Ry = [[math.cos(p), 0, math.sin(p)], [0, 1, 0], [-math.sin(p), 0, math.cos(p)]]
    Rz = [[math.cos(y), -math.sin(y), 0], [math.sin(y), math.cos(y), 0], [0, 0, 1]]
    return [[sum(Rz[i][k]*Ry[k][j] for k in range(3)) for j in range(3)] for i in range(3)]
 
def random_rotation(rng):
    w, x, y, z = (rng.gauss(0, 1) for _ in range(4))
    n = math.sqrt(w*w + x*x + y*y + z*z); w, x, y, z = w/n, x/n, y/n, z/n
    return [[1-2*(y*y+z*z), 2*(x*y-z*w), 2*(x*z+y*w)],
            [2*(x*y+z*w), 1-2*(x*x+z*z), 2*(y*z-x*w)],
            [2*(x*z-y*w), 2*(y*z+x*w), 1-2*(x*x+y*y)]]
 
FACES = [(0, 1, 2), (0, 1, 3), (0, 2, 3), (1, 2, 3)]
 
def centre(X): return [sum(p[i] for p in X) / len(X) for i in range(3)]
 
def bluntness(X):
    com = centre(X); best = -1.0
    for a, b, c in FACES:
        n = cross(sub(X[b], X[a]), sub(X[c], X[a]))
        if dot(n, sub(X[a], com)) < 0: n = [-v for v in n]   # normal exterior
        best = max(best, -n[0] / norm(n))                    # n · <-1,0,0>
    return best
 
def asymmetry(X, rv):
    com = centre(X); d = [sub(p, com) for p in X]
    w = [v[0]**2 if v[0] < 0 else 0.0 for v in d]          # peso cuadrático solo para las esferas aguas arriba
    sw = sum(w)
    ly = sum(wi * v[1] for wi, v in zip(w, d)) / sw
    lz = sum(wi * v[2] for wi, v in zip(w, d)) / sw
    return math.hypot(ly, lz) / rv
 
def polar_angles(X):
    com = centre(X)
    return [math.degrees(math.acos(-(p[0]-com[0]) / norm(sub(p, com)))) for p in X]
 
# curva aproximada transcrita de la descripción de la Fig. 16(a): ángulo polar inicial θ0 -> velocidad lateral final V'_T
CURVE = [(0, 0.10), (70, 0.70), (90, 0.70), (120, 0.30), (160, 0.50), (180, 0.20)]
def v_of_theta(th):
    for (t0, v0), (t1, v1) in zip(CURVE[:-1], CURVE[1:]):
        if th <= t1: return v0 + (v1 - v0) * (th - t0) / (t1 - t0)
    return CURVE[-1][1]
 
X0, rc = tetra_centres()
rv = rc + R_SPH                            # radio de la esfera virtual que envuelve el cúmulo
print(f"rc(centres) = {rc:.4f} r,  rv = {rv:.4f} r")
 
rng = random.Random(7)
N = 40000
B, A, V, band = [], [], [], 0
for _ in range(N):
    R = random_rotation(rng)
    X = [rotate(R, p) for p in X0]
    th = polar_angles(X)
    B.append(bluntness(X)); A.append(asymmetry(X, rv))
    V.append(sum(v_of_theta(t) for t in th) / 4)
    band += sum(60 < t < 90 for t in th)
print(f"random {N}:  B in [{min(B):.4f}, {max(B):.4f}],  A in [{min(A):.4f}, {max(A):.4f}]")
print(f"  P(B > 0.9) = {sum(b > 0.9 for b in B)/N:.3f},  P(B < 0.4) = {sum(b < 0.4 for b in B)/N:.3f}")
print(f"  spheres with 60-90 deg polar angle: {band/(4*N):.3f}")
 
print("\npitch yaw |   B      A    | theta0 of 4 spheres         | mean V'")
gB, gA = [], []
for p in range(-36, 55, 18):
    for y in range(0, 91, 18):
        X = [rotate(attitude(p, y), q) for q in X0]
        b, a, th = bluntness(X), asymmetry(X, rv), sorted(polar_angles(X))
        gB.append(b); gA.append(a)
        if y in (0, 36, 90):
            print(f"{p:5d} {y:3d} | {b:.3f}  {a:.3f} | " + " ".join(f"{t:6.1f}" for t in th)
                  + f" | {sum(v_of_theta(t) for t in th)/4:.3f}")
uniq = {(round(b, 3), round(a, 3)) for b, a in zip(gB, gA)}
print(f"\ngrid 36 attitudes: {len(uniq)} distinct (B, A),  B in [{min(gB):.3f}, {max(gB):.3f}],  A max {max(gA):.3f}")
 
def corr(x, y):
    mx, my = sum(x)/len(x), sum(y)/len(y)
    sxy = sum((a-mx)*(b-my) for a, b in zip(x, y))
    return sxy / math.sqrt(sum((a-mx)**2 for a in x) * sum((b-my)**2 for b in y))
 
print("\nB bin        n     mean A   mean V'")
edges = [1/3, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0001]
for lo, hi in zip(edges[:-1], edges[1:]):
    idx = [k for k in range(N) if lo <= B[k] < hi]
    print(f"{lo:.2f}-{min(hi, 1):.2f} {len(idx):6d}   {sum(A[k] for k in idx)/len(idx):.3f}    "
          f"{sum(V[k] for k in idx)/len(idx):.3f}")
print(f"corr(B, V') = {corr(B, V):.3f},  corr(A, V') = {corr(A, V):.3f}")
print(f"orientation-averaged V' = {sum(V)/N:.3f}")

rc(centres) = 1.2247 r,  rv = 2.2247 r
random 40000:  B in [0.3363, 1.0000],  A in [0.0020, 0.3637]
  P(B > 0.9) = 0.204,  P(B < 0.4) = 0.009
  spheres with 60-90 deg polar angle: 0.251
 
pitch yaw |   B      A    | theta0 of 4 spheres         | mean V'
  -36   0 | 0.750  0.274 |   36.0   77.9  105.6  138.6 | 0.498
  -36  36 | 0.884  0.255 |   49.1   70.4   96.1  152.1 | 0.575
  -36  90 | 0.471  0.184 |   19.5   90.0  118.1  118.1 | 0.404
  -18   0 | 0.569  0.170 |   18.0   93.7  108.5  124.7 | 0.420
  -18  36 | 0.738  0.244 |   39.7   72.7  109.2  137.5 | 0.493
  -18  90 | 0.471  0.184 |   19.5   90.0  118.1  118.1 | 0.404
    0   0 | 0.333  0.000 |    0.0  109.5  109.5  109.5 | 0.355
    0  36 | 0.547  0.230 |   36.0   73.5  123.1  123.1 | 0.435
    0  90 | 0.471  0.184 |   19.5   90.0  118.1  118.1 | 0.404
   18   0 | 0.569  0.170 |   18.0   93.7  108.5  124.7 | 0.420
   18  36 | 0.738  0.244 |   39.7   72.7  109.2  137.5 | 0.493
   18  90 | 0.471  0.184 |   19.5   90.0  118.1  118.1 | 0.404
   36   0 | 0.750  0.274 |   36.0   77.9  105.6  138.6 | 0.498
   36  36 | 0.884  0.255 |   49.1   70.4   96.1  152.1 | 0.575
   36  90 | 0.471  0.184 |   19.5   90.0  118.1  118.1 | 0.404
   54   0 | 0.856  0.153 |   54.0   62.3  101.3  148.9 | 0.548
   54  36 | 0.970  0.201 |   61.6   66.7   84.3  165.9 | 0.603
   54  90 | 0.471  0.184 |   19.5   90.0  118.1  118.1 | 0.404
 
grid 36 attitudes: 21 distinct (B, A),  B in [0.333, 0.970],  A max 0.347
 
B bin        n     mean A   mean V'
0.33-0.50   2466   0.098    0.379
0.50-0.60   5438   0.184    0.422
0.60-0.70   7932   0.193    0.456
0.70-0.80   8052   0.224    0.500
0.80-0.90   7965   0.255    0.554
0.90-1.00   8147   0.243    0.600
corr(B, V') = 0.987,  corr(A, V') = 0.511
orientation-averaged V' = 0.504

Lo que dice 0.987 y lo que se le escapa#

Sobre orientaciones aleatorias, $B$ cubre $[0.336, 1.000]$ . Coincide con el límite tetraédrico $(1/3, 1)$ que indica el artículo. Solo el 0.9 % de las orientaciones tiene $B<0.4$ . Un vértice que apunte exactamente aguas arriba es raro.

El máximo de $A$ resultó 0.364. El artículo da cerca de 0.32. La causa probable es la elección de $r_v$ : radio circunscrito de los centros más un radio de esfera. Una pequeña diferencia de definición produce fácilmente una brecha de ese tamaño.

La velocidad media del cúmulo obtenida con la curva polar se correlaciona con $B$ en 0.987. Con $A$ , solo 0.511. En el fondo, $B$ es la distribución de cuatro ángulos polares plegada en un número. Con una cara de frente, tres esferas quedan a 70.5°, justo en la franja de surfing. Por eso la Fig. 18(a) del artículo muestra la velocidad colectiva subiendo de 0 a 0.5 con $B$ .

Leídas hacia abajo, las medias por intervalo suben con $B$ en una sola dirección: 0.379, 0.422, 0.456, 0.500, 0.554 y 0.600. La diferencia entre extremos es 0.221. Esa monotonía es la correlación de 0.987. El promedio sobre orientaciones, 0.504, queda justo encima del 0.500 del intervalo 0.70-0.80. El peso de la muestra aleatoria no está en el vértice afilado, sino un poco más romo que el medio.

Los conteos no son parejos. $B \lt 0.5$ tiene 2466 orientaciones y $B \gt 0.9$ tiene 8147. Sobre 40 000 eso es 0.062 y 0.204. Coincide con el $P(B \gt 0.9)=0.204$ impreso, y $P(B \lt 0.4)=0.009$ es la cola de un vértice casi de frente. La media de $A$ por intervalo recorre 0.098, 0.184, 0.193, 0.224, 0.255, 0.243. Hace pico en 0.80-0.90 y baja en el último, así que la correlación de $A$ con $V'$ se queda en 0.511. Cuanto más de frente queda la cara, más se reduce otra vez el desplazamiento lateral.

La fila $(0°,0°)$ muestra la simetría que el promedio tapa. Ángulos polares 0.0, 109.5, 109.5, 109.5, con $B=0.333$ , $A=0.000$ y $V'$ media de 0.355. En CURVE, 0° vale 0.10, 90° vale 0.70 y 120° vale 0.30. A 0.65 del tramo entre 90° y 120°, 109.5° baja de 0.70 hacia 0.30 y cae en 0.44. Entonces $(0.10+3\times0.44)/4=0.355$ , el valor impreso. El cálculo real del artículo deja esta orientación cerca de cero. La curva cobra 0.44 a cada una de las tres esferas y no ve el contacto que anula al trío. El mayor $V'$ de la malla es 0.603 en pitch 54°, yaw 36°, con $B=0.970$ . $B=1$ sigue fuera de alcance.

La curva también deja algo fuera. En la orientación $(0°, 0°)$ predice 0.355. En la simulación real del artículo, este cúmulo casi no se separó. Las tres esferas traseras quedaron simétricas a 109.5° y las fuerzas de contacto anularon la presión. El ángulo polar cuenta el destino de un fragmento. No ve el bloqueo colectivo que crea la simetría.

Los valores absolutos tampoco son fiables. El promedio de la curva sobre orientaciones es 0.504, mientras que la media ponderada por área del artículo ronda 0.35. Aun así, supera con claridad el promedio de los pares, cerca de 0.22. Con una curva aproximada, conviene leer solo orden y correlación.

Solo 21 de los 36 puntos de la malla son distintos#

Con el orden de rotación de este artículo (primero pitch, luego yaw), las seis orientaciones de la fila yaw 90° colapsan en el mismo cúmulo. Los $\pm$ pitch son imágenes especulares con el mismo $(B, A)$ . Así, solo 21 de 36 puntos difieren. La malla tampoco llega a $B=1$ ; su máximo es 0.970.

La convención de ángulos del artículo puede ser otra. La lección es la misma. Una malla de ángulos de orientación no cubre el espacio físico de forma pareja. Por eso el artículo trazó el borde del mapa con muestras aleatorias y volvió a dibujar los resultados en el plano $A$ – $B$ . Antes conviene contar el intervalo que las muestras cubren de verdad.

Con 13 esferas, la orientación pesa menos#

La §5 pasa a un cúmulo cúbico centrado en las caras de 13 esferas. Ahora una esfera queda totalmente encerrada. La forma general de la tendencia con el ángulo polar se conserva. La media de $V_T'$ es 0.474 tras la fase primaria y 0.551 al final. El índice de asimetría cae al orden de $10^{-7}$ y deja de servir.

La conclusión: los efectos de orientación se debilitan al crecer la población. Trece sigue siendo poco. Los autores esperan que la dependencia de la orientación casi desaparezca en cúmulos mayores. En los tetraedros, la velocidad lateral máxima de 0.82 vino de un trío con sustentación, unas 2.5 veces la media.

Al llevar estas cifras a un modelo de dispersión en tierra#

Con pocos fragmentos, la orientación decide el resultado. Una distribución condicionada a $B$ es mejor que una sola media. El propio centro de masa se desplaza lateralmente unos 0.13, lo que cambia dónde se deposita la energía.

La velocidad de un fragmento se puede aproximar con la curva polar. Las disposiciones simétricas requieren una restricción de contacto aparte. La hipótesis de dirección radial se sostiene en buena medida. Dejar la cara de contacto dentro de una celda es el problema de reconstrucción de MUSCL-THINC-BVD, que se queda con el salto menor. Más allá de 13 fragmentos, el tamaño de la población y una sola media pueden bastar.