[論文レビュー] 球を一つ前に出した四面体は散らばらなかった

破片一つの横速度はクラスタ内の初期極角でほぼ決まり、クラスタ全体の分離速度はその極角を一つの数にたたんだ鈍さ B で決まります。

球を一つ前に出した四面体は散らばらなかった

同じ大きさの鉄球4個を四面体に組み、Mach 20 の気流に置きます。頂点を一つ正面に向けて手を離します。直感では、すぐにはじけるように散るはずです。Whalen–Deiterding–Laurence(JFM 2026)の計算では、このクラスタはほとんど散りませんでした。後ろの3球を外へ押す圧力を、接触力がちょうど打ち消したためです。

本記事ではこの論文の4球四面体の節(§4)を読みます。破片一つの横速度を何が決めるのか、クラスタ全体はどんな一つの数で要約できるのかを扱います。論文が定義した鈍さ指数 $B$ と非対称指数 $A$ を Python で再計算します。ランダムな姿勢4万個で、二つの指数の範囲と相関を確かめます。

球2個が 144° 付近でくっついたまま揚力を出す話は前回のレビューで扱いました。今回は球が4個に増えた場合です。

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 では3点が 70° 付近の平坦部に集まります。白い円が $B=1/3$ と $B=1$ の間を動くかを見てください。

1980年 Passey–Melosh: 破片2個と定数1個#

流星体の破片モデルの出発点は、Passey と Melosh(1980)の二体モデルです。2球が真横にだけ押し出されると考え、最終横速度を一つの式で与えました。大きさは $V_T \propto V\sqrt{\rho_a/\rho_m}$ でまとまります。残った無次元定数 $C$ は観測ごとに大きくばらつきました。

もう一つの系統はパンケーキモデル(強度のない塊が平たく広がると見るモデル)です。このモデルは、破片が純粋に半径方向へ広がると仮定します。その仮定を検証した計算はほとんどありませんでした。論文は二つの問いを同時に立てます。破片一つの速度は何で決まるか。その向きは本当に半径方向か。

時間と速度の尺度

論文は結果を二つの尺度で無次元化します。

\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$ は4球の平均です。流れは AMROC(適応格子 Euler ソルバー)が、球の接触と弾性衝突は DYNA3D(陽解法有限要素コード)が解きます。

pitch・yaw の代わりに B と A#

四面体の姿勢は pitch と yaw の二つの角で決まります。論文はこれらの角に物理的な意味はないと考え、二つの指数を新たに定義します。

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

$\hat{n}_i$ は3つの球中心を通る外側の面の単位法線です。流れは $+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° 以上では再び上がります。前の球の陰に隠れていて、2球の相互作用が始まる球たちです。180° 付近は対称性のため低い値です。

向きも確認されています。75% の球は、方位角が初期値から 20° 以内にとどまりました。同じ大きさの破片なら、パンケーキモデルの半径方向膨張の仮定はおおむね成り立つということです。

コード: 四面体の姿勢4万個#

$B$ 、 $A$ 、極角は姿勢だけで決まります。流れを解かずに直接数えられます。以下のコードは一様ランダムな回転4万個と、論文式の pitch–yaw 格子36個を回します。各球の速度は、上の記述を写した概略曲線 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      # 接触した球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]          # 上流側の球だけに2次の重み
    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$ は結局、4球の極角分布を一つの数にたたんだものです。面が正面を向くと3球が 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 に置きます。曲線は3球をそれぞれ 0.44 と数え、接触がその三つを束ねて打ち消すことは数えません。格子で最大の $V'$ は pitch 54°、yaw 36° の 0.603 で、そのとき $B=0.970$ です。 $B=1$ にはまだ届きません。

曲線が見落とす場所もあります。 $(0°, 0°)$ の姿勢での予測値は 0.355 です。論文の実際の計算では、このクラスタはほとんど分離しませんでした。後ろの3球がすべて 109.5° に対称に並び、接触力が圧力を打ち消しました。極角は破片一つの運命を語ります。対称性が生む集団の拘束は語れません。

絶対値も信用しにくいです。曲線の姿勢平均は 0.504 ですが、論文の面積加重平均は約 0.35 です。それでも2球ペアの平均約 0.22 よりは明らかに大きいです。概略曲線からは、順序と相関だけを読むのが安全です。

格子36点のうち異なるのは21点だけ#

本記事の回転順序(pitch が先、yaw が後)では、yaw 90° の行の6姿勢がすべて同じクラスタになります。 $\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 は、揚力を出す3球の組から出ました。平均の約 2.5 倍です。

落下分布モデルにこの数字を持ち込むとき

破片が少ないと姿勢が結果を分けます。このときは平均一つより、 $B$ に応じた分布を使うほうが良いです。重心そのものも約 0.13 だけ横にずれます。エネルギーを落とす位置が変わるということです。

破片一つの速度は極角曲線で近似できます。対称な配置では接触拘束が別に必要です。向きについては半径方向の仮定がおおむね持ちこたえます。接触面を1セルに立てる再構成は、より小さい振幅を選ぶ MUSCL-THINC-BVDの問題です。13個を超えるクラスタなら、姿勢の代わりに個数と平均一つで十分かもしれません。

[論文レビュー] 球を一つ前に出した四面体は散らばらなかった — Mach 20 破片クラスタの鈍さ指数