[论文评述] 两块碎片没有互相推开，而是贴在一起飞行

决定碎片对分离速度的不是斥力的大小，而是维持接触的力矩的符号。

假设一颗陨石在大气层中裂成了两块。这两块会彼此远离吗？如果答案永远是"会"，那么从地面的陨石坑群反推入射条件就轻松多了。Whalen、Deiterding 与 Laurence 发表在 JFM 2026（vol. 1029, A35）的 Mach 20 计算动摇了这个答案。在某些对齐角下，两个球会贴在一起同行，甚至整个球对还会产生升力。本文用 30 行的 Newton 近似复现这个角度是从哪里来的。算出来是 143.6°，落在论文给出的稳定区间 132°–145.7° 之内。

Q1. 只有两块碎片的问题为什么至今仍未解决#

陨石碎片的计算长期存在两条路线。碎块只有几个时，用逐块追踪的 discrete-fragment（离散碎片）方法。碎块多到近乎无穷时，用把整团当作液体铺开的 debris-cloud（碎片云）方法。Passey 与 Melosh 在 1980 年提出的双体模型是前者的原型。若假设两球只沿横向被推开，最终横向速度可整理成下面的比例关系。

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

$V$ 是进入速度， $\rho_a$ 是大气密度， $\rho_m$ 是陨石密度。剩下的只有一个比例常数，而从地面陨石坑群反推出的值从 0.03 一直散到 2.28。整整两个数量级。

麻烦出在中间地带。碎块既不是 2 个也不是几千个的中等数量区间几乎无人涉足。这篇论文正面瞄准了这一段。作者把同尺寸的球按 2 个、4 个、13 个规则排列，只改变初始姿态，算了 83 个算例。流场由 AMROC 用 embedded boundary 求解 Euler 方程，球体由 DYNA3D 连同接触一起求解。气体取 $\gamma=1.4$ 的完全气体，来流 Mach 20，球与气体的密度比为 $10^4$ 。网格沿密度梯度自动加密。这套做法的成本结构与 AMR 打标与重通量中讨论的一致。

在下面的模拟里亲手操作一下。

theta_0170°contact damping0.08

theta 170.0°C_M 0.0000L/D 0.000C_D 0.000omega 0.00°/tequilibrium 0.0°

Set theta_0 anywhere and press release. Watch the curved torque arrow: it is green wherever the pair is being opened and red wherever it is being closed, and it vanishes only on the dashed line. Drop damping to zero and the pair never settles — it swings through the equilibrium angle forever. Start at exactly 180° and nothing happens at all, until you nudge it.

把对齐角 $\theta_0$ 放在任意位置，按下 release 即可。弯箭头为绿色表示球对张开的方向，红色表示闭合的方向。右侧曲线与 0 相交的那一点，会把所有初始角度都吸过去。

Q2. 贴在一起的球对为什么偏偏走向某一个角度#

对齐角 $\theta$ 定义为自由来流方向与"次级球 → 主球"连线之间的夹角。 $\theta = 180^\circ$ 时后球正好坐进前球的尾迹里。 $\theta = 90^\circ$ 时两球并排。

只要两球保持接触，就可以当成一个刚体来看。此时接触力属于内力，在球对整体的力矩中相互抵消。剩下的只有气动力矩。

\mathbf{M}_{\rm com} = \sum_{i=1}^{2}\oint_{S_i}(\mathbf{x}-\mathbf{x}_{\rm com})\times\bigl(-C_p\,\hat{n}\bigr)\,\mathrm{d}A

$S_i$ 是各球的表面， $\hat{n}$ 是外法线， $\mathbf{x}_{\rm com}$ 是球对的质心。平衡条件只要求这个量为 0，而稳定性由它附近的符号决定。

先看 $\theta = 180^\circ$ 。上游球正面迎着自由来流，下游球被完全遮挡，阻力几乎为 0。两球的阻力差作用在偏离质心的位置上，因此角度只要稍有偏转，就会产生放大这一偏转的力矩。这是不稳定平衡。论文也写到 180° 下未观测到姿态变化，并补充说由于上游球阻力更大，这个状态应当是不稳定的。

另一端的 $\theta = 90^\circ$ 因为对称，力矩恰好为 0。那么在两者之间必然存在一个符号翻转的位置。那个位置就是稳定平衡角。

Q3. 仅凭 Newton 近似能不能算出这个角度#

在高超声速下，压力分布可以用 Newton 撞击理论来近似。若认为流动微团撞上表面后完全损失法向动量，压力系数就是下式。

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

$\hat{u}$ 是自由来流单位矢量， $\hat{n}$ 是表面外法线，只在 $\hat{u}\cdot\hat{n}<0$ 的迎风面上有效。 $C_{p,\max}$ 是 $\gamma=1.4$ 、 $M\to\infty$ 时的驻点值。在此基础上只再加一行遮挡判定：某个面元若藏在另一个球后面，压力就置 0。

把球面用 Fibonacci 格点切成 2.4 万个面元，从 90° 到 180° 以 1° 步长扫描。

import math
 
CP_MAX = 1.8394          # 修正 Newton 理论，gamma=1.4，M -> 无穷大
R = 1.0                  # 球半径
N_PANEL = 24000
 
 
def fib_sphere(n):
    """单位球面上近似均匀分布的点 + 单个面元的面积。"""
    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):
    """次级（下游）球位于原点，主球沿 n_hat 方向偏移 2R。
    theta = 自由来流 x_hat 与 '次级 -> 主球' 连线之间的夹角。"""
    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):
    """该点是否被位于 c 的球遮挡而看不到自由来流？
    来流沿 +x 流动，因此沿 -x 回溯做圆柱判定。"""
    if cx >= px:
        return False
    return (py - cy) ** 2 + (pz - cz) ** 2 < R * R
 
 
def newtonian_cp(nx):
    """nx = 外法线的 x 分量。nx < 0 表示迎风面。"""
    return CP_MAX * nx * nx if nx < 0.0 else 0.0
 
 
def pair_loads(theta_deg):
    """各球的力系数，以及相对球对质心的力矩。
    参考量：力用 q_inf * pi * R^2，力矩再乘以 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     # 积分噪声。真正的过零比这陡峭得多
 
 
def stable_window(rows):
    """只保留 C_M 零点中斜率超过噪声阈值的那些。"""
    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

结果是 143.61°。论文翻遍系数与力矩历程找出的稳定区间是 132°–145.7°，接触球对的升阻比最大点在 141.5°。没有网格、也没有激波的 30 行代码落进了这个区间。

升阻比只有 0.119，低于论文稳定区间的平均值 0.197 和最大值 0.22。差距的原因很清楚：Newton 近似不知道激波与激波之间的干扰。论文插图中两球之间内侧面上出现的高压斑点，是两道弓形激波相遇造成的，而这个模型里根本没有那个位置。正如在斜激波串中看到的，激波叠加处的压力高于简单相加。所以这段代码算准的不是量值，而是符号和位置。

表里符号翻转的位置也值得留意。 $C_L$ 在 150° 处达到最大的 0.215，而 $C_M$ 却在更早的 143.6° 就过了 0。升力最大的姿态与能自我维持的姿态并不相同。

Q4. 分开之后把它们推向侧向的力从哪来#

一旦接触断开，故事就变了。若两球从一开始就并排站立（90° 附近），共同的弓形激波压在内侧面上，把两者推开。这种相互排斥在间距拉开后很快消失，论文中的无量纲横向速度停在 0.2 上下。

但当后球骑上前球的激波时，情况就不同了。激波正后方压力高，外侧则是自由来流。球体若一半压着这道边界前进，向外推的力就会持续作用。这就是 Laurence 与 Deiterding 在 2011 年命名的 shock surfing（激波冲浪）。论文中这一效应把最终横向速度从 0.2 提到 0.25，并把分离所需时间从 $1.7\tau_s$ 拉长到 $4\tau_s$ 。

时间与速度的参考尺度如下。

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

$r_c$ 是团簇的外接半径。密度比为 $10^4$ ，所以 $\tau_s$ 是通过时间的 100 倍。用这个尺度无量纲化之后，相对运动方程只剩下一个系数。

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

$3/8$ 来自球的体积与截面积之比，密度比则被完全吸收进了尺度里。下面的动画就是直接积分这个方程。

release angle150°surf gain0.25speed1.0x

t' 0.00V'_T 0.000no surf 0.000gap 0.00 rtheta 150°free

Drag the release angle from 90° upward. Near 90° the two curves sit on top of each other — the shared bow shock does all the work and surfing has nothing to add. Past about 135° the amber curve lifts off the dashed one, and the sphere outline glows exactly while it is straddling the pink shock line. At 180° nothing separates at all.

把 release angle 从 90° 往上拖。90° 附近实线与虚线重合在一起，意思是共同激波干完了全部的活，冲浪没什么可加的。越过 135° 之后两条曲线开始分岔，而球的边缘只有在咬住粉色激波线的那段时间里才会变亮。

接触阶段的行为也包含在这个模型里。从 165° 释放时，两球会贴在一起滚动 3–4 $\tau_s$ ，直到对齐角降到 120° 附近才分开。这与论文报告的接触持续 4–6 $\tau_s$ 、脱离接触角在 130° 附近大致吻合。阻力在尾迹中骤降的方式，与阻力危机中看到的单球阻力曲线是两码事。这里的原因不是分离点，而是遮挡。

Q5. 那么 180° 处为什么什么都不发生#

精确对齐的串列是对称的。横向力为 0，力矩也为 0。计算的结果是两球在 8 $\tau_s$ 内始终贴合、姿态毫无变化。论文的观测也一样。

但这就像把一支铅笔竖着立起来的那种平衡。上游球的阻力大于下游球（上表中 0.920 对 0.000），因此只要稍微倾斜，就会附上放大这一倾斜的力矩。实际计算中，即使从 172.5° 释放，初期也几乎不动，随后越来越快，呈现出"滚动"式的运动，最终接触断开。

这里出现了一个在工程上很重要的量。保持接触并产生升力的球对，其质心本身会被推向侧向。论文写到，从 150°、157.5°、172.5° 释放的球对，其质心横向速度分别为 0.42、0.32、0.28，大于单个碎片之间的相对速度。若只统计每块碎片各自的排斥，这个分量会被整块漏掉。

推到 13 个之后剩下什么#

论文继续用 4 球四面体排列的 38 种、13 球面心立方排列的 34 种做同样的考察。单个球的最终横向速度可以用一个初始极角相当好地归并，这一趋势在 4 球和 13 球中表现相近。不过团簇整体钝度对体相行为的影响，会随个数增加而变淡。碎块越多，分离就越均质，这也提示了 debris-cloud 方法从哪里开始变得可用。

归纳起来，支配碎片对分离的不是斥力的大小，而是接触何时断开；决定这一点的是力矩的符号。而这个符号不需要网格也能知道。把表面压力交给 Newton 近似，只要遮挡判得准，两个球自己找到的角度是 143.6° 这个答案，30 行之内就出来了。

[论文评述] 两块碎片没有互相推开，而是贴在一起飞行 — Mach 20 球体对的 144° 平衡角