[论文评述] 把一个球顶在前面的四面体没有散开

单个碎片的横向速度几乎由它在簇中的初始极角决定，整个簇的分离速度则由把这些极角折叠成一个数的钝度 B 决定。

把一个球顶在前面的四面体没有散开

把四个同样大小的铁球粘成四面体，放进 Mach 20 气流。让一个顶点正对来流，然后松手。直觉告诉我们，它会像爆开一样散开。在 Whalen–Deiterding–Laurence(JFM 2026)的计算中，这个簇几乎没有散开。把后方三个球向外推的压力，恰好被接触力抵消了。

本文阅读该论文中关于四球四面体的一节(§4)。单个碎片的横向速度由什么决定？整个簇能用哪一个数来概括？我们用 Python 重新计算论文定义的钝度指数 $B$ 和非对称指数 $A$ 。并在4万个随机姿态上检查两个指数的范围和相关性。

两个球在 144° 附近贴在一起飞行并产生升力的故事，在上一篇评述中讲过。这次球增加到四个。

pitch 和 yaw 滑块移动 $A$ – $B$ 平面上的白圈，roll 则把白圈留在原地。

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 把极角分成 0° 和 109.5°，face-first 让三个点落在 70° 附近的平台上。看白圈是否在 $B=1/3$ 和 $B=1$ 之间来回。

1980 年 Passey–Melosh：两个碎片，一个常数#

流星体碎片模型的起点是 Passey 和 Melosh(1980)的二体模型。他们假设两个球只沿横向被推开，并用一个公式给出最终横向速度。其量级为 $V_T \propto V\sqrt{\rho_a/\rho_m}$ 。剩下的无量纲常数 $C$ 在不同观测之间差别很大。

另一条路线是薄饼模型(pancake model，把无强度的团块看作被压扁并向外铺开)。该模型假设碎片纯粹沿径向扩张。几乎没有计算检验过这个假设。论文同时提出两个问题。单个碎片的速度由什么决定？它的方向真的是径向吗？

时间与速度的尺度

论文用两个尺度对结果做无量纲化。

\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}$ 是球的密度， $\rho_\infty$ 和 $u_\infty$ 是来流密度和速度， $r_c$ 是簇的外接半径。密度比为 $10^4$ ，所以 $v_s = u_\infty/100$ 。对于 20 km/s 的流星体， $V_T' = 0.35$ 相当于横向 70 m/s。

横向速度 $V_{T,i}$ 是相对质心的速度在 $y$ – $z$ 平面径向上的投影。簇的集体分离速度 $V_T$ 是四个球的平均值。流场由 AMROC(自适应网格 Euler 求解器)计算，球的接触与弹性碰撞由 DYNA3D(显式有限元程序)处理。

用 B 和 A 代替 pitch 和 yaw#

四面体的姿态由 pitch 和 yaw 两个角确定。论文认为这两个角没有物理意义，于是定义了两个新指数。

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

$\hat{n}_i$ 是经过三个球心的外表面的单位外法线。流动沿 $+x$ 方向。一个顶点在前时 $B=1/3$ ，一个面正对来流时 $B=1$ 。 $B$ 越大，越多的球共同形成弓形激波(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$ 表示横向分量， $w_i'$ 是归一化后和为 1 的权重， $r_v$ 是虚拟外接球半径。上游球偏离轴线越远， $A$ 越大。

一个极角决定一个碎片

论文按初始极角 $\theta_0$ (从质心看与上游轴的夹角)把所有球分成 15° 一组。结果很清楚。最终横向速度从 0° 到约 70° 几乎线性上升。在 70°–90° 经过一段短暂平台，接近 0.7。这一区间最容易出现激波冲浪(shock surfing，沿前方球的激波滑行的轨迹)。不让这道激波糊成好几格的重构，是Fu(2019)的 TENO单独处理的问题。

超过 90° 后，被前方球尾流卷入的情况增多，速度下降。120° 以上又重新上升。这些球先躲在前方球身后，直到两球相互作用开始。180° 附近由于对称性，速度较低。

方向也做了检查。75% 的球方位角变化不超过初始值 20°。这说明对于等大碎片，薄饼模型的径向扩张假设大体成立。

代码：4万个四面体姿态#

$B$ 、 $A$ 和极角只取决于姿态，不需要求解流场就能直接统计。下面的代码遍历4万个均匀随机旋转，以及仿照论文的 36 点 pitch–yaw 网格。每个球的速度由示意曲线 CURVE 给出。这条曲线不是图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      # 四个相切的球: 中心四面体的外接半径
    s = rc * math.sqrt(8) / 3              # 底面圆的半径
    X = [[-rc, 0.0, 0.0]]                  # 一个顶点朝向上游(-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]   # 外法线
        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]          # 只对上游的球施加二次权重
    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]
 
# 按论文正文对图16(a)的描述转写的示意曲线: 初始极角 θ0 -> 最终横向速度 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                            # 包围簇的虚拟球半径
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

0.987 说明了什么，又漏掉了什么#

在随机姿态下， $B$ 覆盖 $[0.336, 1.000]$ 。这与论文给出的四面体范围 $(1/3, 1)$ 一致。 $B<0.4$ 的姿态只占 0.9%。顶点恰好正对来流的姿态很少见。

$A$ 的最大值为 0.364，论文给出约 0.32。原因可能是我们把 $r_v$ 取为中心外接半径加一个球半径。定义细节稍有不同，就很容易产生这样的差距。

只用极角曲线得到的簇平均速度与 $B$ 的相关系数为 0.987，与 $A$ 只有 0.511。 $B$ 本质上是把四个球的极角分布折叠成一个数。一个面正对来流时，三个球位于 70.5°，正好落在冲浪区间。这就是论文图18(a)中集体分离速度随 $B$ 从 0 升到 0.5 的原因。

把分箱平均往下读， $V'$ 几乎只顺着 $B$ 往一个方向升。0.33-0.50 是 0.379，0.50-0.60 是 0.422，0.60-0.70 是 0.456，0.70-0.80 是 0.500，0.80-0.90 是 0.554，0.90-1.00 是 0.600。两端相差 0.221。相关为 0.987，就是因为这条单调。姿态平均 0.504 坐在 0.70-0.80 一档的 0.500 正上方。随机样本的重心不在尖顶点上，而在比中间略钝的一侧。

样本数并不均匀。 $B \lt 0.5$ 有 2466 个， $B \gt 0.9$ 有 8147 个。除以 4 万是 0.062 和 0.204。输出里的 $P(B \gt 0.9)=0.204$ 对得上， $P(B \lt 0.4)=0.009$ 是顶点几乎正对来流的那条尾巴。 $A$ 的分箱平均是 0.098、0.184、0.193、0.224、0.255、0.243，在 0.80-0.90 见顶，最后一档又掉下来。所以 $A$ 和 $V'$ 的相关停在 0.511。面越接近正对，横向偏移重新变小。

$(0°,0°)$ 那一行露出这个平均盖住的对称。极角是 0.0、109.5、109.5、109.5， $B=0.333$ ， $A=0.000$ ，平均 $V'=0.355$ 。CURVE 上 0° 是 0.10，90° 是 0.70，120° 是 0.30。109.5° 落在这两个刻度之间 0.65 的位置，从 0.70 走向 0.30，得到 0.44。 $(0.10+3\times0.44)/4=0.355$ ，和输出一致。论文的实际计算把这个姿态几乎放在 0。曲线把三个球各算成 0.44，并不计入接触把三者捆在一起抵消的事。网格上最大的 $V'$ 是 pitch 54°、yaw 36° 的 0.603，那时 $B=0.970$ 。 $B=1$ 还够不着。

曲线也有漏掉的地方。在 $(0°, 0°)$ 姿态下，它预测 0.355。而论文的实际计算中，这个簇几乎没有分离。后方三个球都对称地位于 109.5°，接触力抵消了压力。极角能说明单个碎片的命运，却看不到对称性造成的集体约束。

绝对值同样不可靠。曲线的姿态平均为 0.504，而论文的面积加权平均约为 0.35。不过这仍明显高于两球对的平均值约 0.22。用示意曲线时，只读顺序和相关性更稳妥。

36 个网格点中只有 21 个不同#

在本文的旋转顺序下(先 pitch，后 yaw)，yaw 90° 那一行的六个姿态全部变成同一个簇。 $\pm$ pitch 互为镜像，给出相同的 $(B, A)$ 。因此 36 点中只有 21 点不同。网格也到不了 $B=1$ ，最大值是 0.970。

论文的角度约定可能与本文不同，但教训一样。姿态角网格并不能均匀覆盖物理空间。这正是论文用随机样本画出图的边界，并把结果重新铺到 $A$ – $B$ 平面上的原因。应先数清样本实际盖住的区间。

13 个球时姿态没那么重要#

§5 研究面心立方的 13 球簇。此时有一个球被完全包在中心。极角与速度相关性的大致形状得以保留。平均 $V_T'$ 在一次分离后为 0.474，最终为 0.551。非对称指数降到 $10^{-7}$ 量级，失去作用。

结论是，姿态效应随个数增加而减弱。13 个仍然是小数目。作者认为碎片更多时，姿态依赖几乎会消失。在四面体中，最大横向速度 0.82 来自一个产生升力的三球组合，约为平均值的 2.5 倍。

把这些数字带进落点分布模型时

碎片少时，姿态决定结果。此时与其用一个平均值，不如用按 $B$ 条件化的分布。质心本身也会横向偏移约 0.13，这意味着能量沉积的位置会改变。

单个碎片的速度可以用极角曲线近似。对称布置需要单独的接触约束。方向上，径向假设大体成立。把接触面立在一个格子里，是挑选更小振幅的 MUSCL-THINC-BVD的问题。超过 13 个碎片时，用个数和一个平均值或许就够了。

[论文评述] 把一个球顶在前面的四面体没有散开 — Mach 20 碎片簇的钝度指数