[Paper Review] The Tetrahedron That Led With One Sphere Never Broke Apart — A Bluntness Index for Mach 20 Fragment Clusters

The tetrahedron that led with one sphere never broke apart#

Glue four equal iron spheres into a tetrahedron and put them in a Mach 20 stream. Point one vertex straight upstream and let go. Intuition says the cluster should burst apart. In the simulations of Whalen, Deiterding and Laurence (JFM 2026), it barely moved apart. The pressure pushing the three rear spheres outward was balanced exactly by contact forces.

This post reads the four-sphere tetrahedron section (§4) of that paper. What sets the lateral speed of one fragment? What single number summarizes a whole cluster? We recompute the paper's bluntness index $B$ and asymmetry index $A$ in Python. Then we check their range and correlation over 40,000 random attitudes.

The story of two spheres flying stuck together near 144° while generating lift was covered in the previous review. Today the pair becomes a group of four.

The pitch and yaw sliders move the white circle across the $A$–$B$ plane, and roll leaves that circle where it is.

vertex-first splits the polar angles into 0° and 109.5°. face-first piles three of them on the plateau near 70°. Watch whether the circle travels between $B=1/3$ and $B=1$.

1980, Passey–Melosh: two fragments, one constant#

Meteoroid fragment modeling starts with the binary model of Passey and Melosh (1980). They assumed two spheres pushed purely sideways and derived a single formula for the final lateral speed. It scales as $V_T \propto V\sqrt{\rho_a/\rho_m}$. The remaining dimensionless constant $C$ scattered widely from one observation to the next.

The other branch is the pancake model (a strengthless mass assumed to flatten and spread). It assumes fragments expand purely radially. Almost no calculation had tested that assumption. The paper asks both questions. What sets one fragment's speed? Is its direction really radial?

Scales for time and velocity#

The paper nondimensionalizes with two scales.

$\rho_{sph}$ is the sphere density, $\rho_\infty$ and $u_\infty$ the free-stream density and velocity, and $r_c$ the circumscribed radius of the cluster. The density ratio is $10^4$, so $v_s = u_\infty/100$. For a 20 km/s meteoroid, $V_T' = 0.35$ means 70 m/s sideways.

The lateral velocity $V_{T,i}$ is the velocity relative to the center of mass, projected onto the radial direction in the $y$–$z$ plane. The collective separation speed $V_T$ is the average over the four spheres. AMROC (an adaptive-mesh Euler solver) handles the flow. DYNA3D (an explicit finite-element code) handles sphere contact and elastic collisions.

B and A instead of pitch and yaw#

Two angles, pitch and yaw, fix the tetrahedron's attitude. The paper argues these carry no physical meaning and defines two new indices.

$\hat{n}_i$ is the outward unit normal of a face through three sphere centers. Flow runs in $+x$. With one vertex in front, $B=1/3$. With a face pointed upstream, $B=1$. Larger $B$ means more spheres share the bow shock.

$l$ denotes the lateral component, $w_i'$ the weights normalized to sum to one, and $r_v$ the radius of the virtual circumscribing sphere. $A$ grows as upstream spheres sit farther off the axis.

One polar angle decides one fragment#

The paper bins every sphere by its initial polar angle $\theta_0$ (the angle from the upstream axis, seen from the center of mass) in 15° groups. The result is clear. Final lateral speed rises almost linearly from 0° to about 70°. It plateaus briefly between 70° and 90° near 0.7. That band is where shock surfing (a sphere riding along another sphere's shock) is most likely. Keeping that shock from smearing across several cells is the reconstruction problem in Fu's 2019 TENO scheme.

Past 90°, more spheres get entrained in a leading sphere's wake and speeds fall. Above 120° they rise again. These spheres hide behind an upstream body until a two-sphere interaction starts. Near 180°, symmetry keeps speeds low.

Direction was checked too. 75% of spheres stayed within 20° of their initial azimuth. For equal-sized fragments, the pancake model's radial-expansion assumption mostly holds.

Code: 40,000 tetrahedron attitudes#

$B$, $A$ and the polar angles depend only on attitude. We can count them without any flow solve. The code below loops over 40,000 uniform random rotations and a 36-point pitch–yaw grid like the paper's. Each sphere's speed comes from the sketch curve CURVE. It is transcribed from the paper's prose, not from the Fig. 16(a) data.

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      # four touching spheres: circumradius of the centre tetrahedron
    s = rc * math.sqrt(8) / 3              # radius of the base circle
    X = [[-rc, 0.0, 0.0]]                  # one vertex points upstream (-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]   # outward normal
        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]          # quadratic weight on upstream spheres only
    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]
 
# sketch curve transcribed from the paper's description of Fig. 16(a): initial polar angle θ0 -> final lateral velocity 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                            # radius of the virtual sphere enclosing the cluster
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

What 0.987 says, and what it misses#

Over random attitudes, $B$ spans $[0.336, 1.000]$. That matches the tetrahedral limit $(1/3, 1)$ stated in the paper. Only 0.9% of attitudes have $B<0.4$. A vertex pointing exactly upstream is rare.

The maximum of $A$ came out at 0.364. The paper says about 0.32. The likely cause is our choice of $r_v$: center circumradius plus one sphere radius. A small difference in definition easily produces a gap this size.

The cluster-mean speed from the polar-angle curve correlates with $B$ at 0.987. With $A$, only 0.511. $B$ is really the distribution of four polar angles folded into one number. With a face upstream, three spheres sit at 70.5°, right in the surfing band. That is why the paper's Fig. 18(a) shows collective speed climbing from 0 to 0.5 with $B$.

Read the bin means downward and $V'$ climbs with $B$ in one direction. The bins give 0.379, 0.422, 0.456, 0.500, 0.554, and 0.600. The end-to-end gap is 0.221. That monotone is why the correlation is 0.987. The attitude average, 0.504, sits just above the 0.500 of the 0.70-0.80 bin. The random sample is weighted past the sharp vertex, a little blunter than the middle.

The counts are not even. $B \lt 0.5$ has 2466 attitudes and $B \gt 0.9$ has 8147. Divided by 40,000 that is 0.062 and 0.204. The printed $P(B \gt 0.9)=0.204$ matches, and $P(B \lt 0.4)=0.009$ is the thin tail of a vertex almost dead upstream. Mean $A$ by bin runs 0.098, 0.184, 0.193, 0.224, 0.255, 0.243. It peaks in 0.80-0.90 and drops in the last bin, so the correlation of $A$ with $V'$ stops at 0.511. Closer to a face-on attitude, the lateral offset shrinks again.

The $(0°,0°)$ row shows the symmetry that average hides. Polar angles 0.0, 109.5, 109.5, 109.5, with $B=0.333$, $A=0.000$, and mean $V'=0.355$. On CURVE, 0° is 0.10, 90° is 0.70, and 120° is 0.30. At 0.65 of the way from 90° to 120°, 109.5° falls from 0.70 toward 0.30 and lands on 0.44. Then $(0.10+3\times0.44)/4=0.355$, the printed value. The paper's actual run leaves this attitude near zero. The curve bills each of the three spheres at 0.44 and never sees contact canceling the trio. The largest grid $V'$ is 0.603 at pitch 54°, yaw 36°, where $B=0.970$. $B=1$ is still out of reach.

The curve also misses something. At attitude $(0°, 0°)$ it predicts 0.355. In the paper's actual simulation this cluster hardly separated. The three rear spheres sat symmetrically at 109.5°, and contact forces canceled the pressure. Polar angle tells one fragment's fate. It cannot see the collective lock that symmetry creates.

Absolute values are not trustworthy either. The curve's attitude average is 0.504, while the paper's area-weighted mean is about 0.35. That is still clearly above the pair average of about 0.22. With a sketch curve, read only ordering and correlation.

Only 21 of 36 grid points are distinct#

With this post's rotation order (pitch first, then yaw), all six attitudes in the yaw 90° row collapse to the same cluster. $\pm$pitch are mirror images with identical $(B, A)$. So only 21 of 36 points differ. The grid never reaches $B=1$ either; its maximum is 0.970.

The paper's angle convention may differ from ours. The lesson is the same. An attitude-angle grid does not cover physical space evenly. That is why the paper drew the map boundary from random samples and replotted the results in the $A$–$B$ plane. Count the interval the samples actually cover before trusting the map.

With 13 spheres, attitude matters less#

§5 moves to a face-centered-cubic cluster of 13 spheres. One sphere is now fully enclosed. The overall shape of the polar-angle trend survives. Mean $V_T'$ is 0.474 after the primary phase and 0.551 at the end. The asymmetry index drops to order $10^{-7}$ and becomes useless.

The conclusion: attitude effects weaken as the population grows. Thirteen is still small. The authors expect attitude dependence to nearly vanish for larger clusters. In the tetrahedra, the maximum lateral speed of 0.82 came from a lifting triplet, about 2.5 times the mean.

Carrying these numbers into a strewn-field model#

With few fragments, attitude decides the outcome. A distribution conditioned on $B$ beats a single mean. The center of mass itself drifts sideways by about 0.13, which moves where the energy is deposited.

One fragment's speed can be approximated with the polar-angle curve. Symmetric arrangements need a separate contact constraint. The radial-direction assumption mostly holds. Keeping a contact face inside one cell is the reconstruction problem in MUSCL-THINC-BVD, which keeps the smaller jump. Beyond 13 fragments, population size and a single mean may be enough.