[論文レビュー] 2つの破片は反発せず、くっついたまま飛んだ — Mach 20 球ペアの144°平衡角

隕石が大気圏で2つに割れたとしましょう。2つの破片は互いに離れていくのでしょうか。答えが「必ずそうなる」であれば、地上のクレーター群から逆算する作業はもっと楽だったはずです。Whalen・Deiterding・LaurenceがJFM 2026 (vol. 1029, A35) に発表したMach 20の計算は、その答えを揺さぶります。ある整列角では2つの球がくっついたまま飛び続け、しかもペア全体が揚力を生みます。この記事では、その角度がどこから来るのかを30行のNewton近似で再現します。出てきた値は143.6°で、論文の安定区間132°–145.7°の中に収まっていました。

Q1. 破片2個の問題がなぜ今も未解決なのか#

隕石破片の計算には古い2つの流れがあります。破片が数個しかない場合は、各破片を個別に追跡するdiscrete-fragment（離散破片）アプローチを使います。破片が事実上無限にある場合は、塊を液体のように押し広げるdebris-cloud（破片雲）アプローチを使います。PasseyとMeloshが1980年に示した二体モデルが前者の原型です。2つの球が純粋に横方向へ押されるだけと置くと、最終的な横方向速度は次の比例式にまとまります。

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

$V$は突入速度、$\rho_a$は大気密度、$\rho_m$は隕石密度です。残るのは比例定数1つだけですが、地上のクレーター群から逆算した値は0.03から2.28まで開きました。2桁の幅です。

問題はその中間にあります。破片が2個でも数千個でもない中間個数の領域は、ほとんど空白でした。論文はそこを正面から狙います。同じ大きさの球を2個・4個・13個に規則配列し、初期姿勢だけを変えて83ケースを計算しました。流れはAMROCがEuler方程式をembedded boundaryで解き、球はDYNA3Dが接触まで含めて解きます。気体は$\gamma=1.4$の完全気体、流入はMach 20、球と気体の密度比は$10^4$です。格子は密度勾配に沿って自動的に細かくなります。この方式のコスト構造はAMRタギングとリフラクシングで扱ったものと同じです。

下のシミュレーションで実際に操作してみましょう。

整列角$\theta_0$を好きな位置に置いてreleaseを押すだけです。曲がった矢印が緑ならペアが開く向き、赤なら閉じる向きを表します。右側の曲線が0を横切る一点へ、あらゆる初期角度が吸い込まれていきます。

Q2. 接触したままのペアはなぜ特定の一角度へ向かうのか#

整列角$\theta$は、自由流の向きと「二次球 → 一次球」を結ぶ線分のなす角として定義します。$\theta = 180^\circ$なら後ろの球が前の球の後流にぴったり収まります。$\theta = 90^\circ$なら横並びです。

2つの球が接触を保っている間は、1つの剛体とみなせます。このとき接触力は内力なので、ペア全体のモーメントからは消えます。残るのは空力モーメントだけです。

$S_i$は各球の表面、$\hat{n}$は外向き法線、$\mathbf{x}_{\rm com}$はペアの質量中心です。平衡条件はこの値が0になることだけで、安定性はその近傍での符号が決めます。

まず$\theta = 180^\circ$から見ます。上流の球は自由流を正面から受け、下流の球は完全に隠れて抗力がほぼ0です。2つの球の抗力差が質量中心から外れた位置に働くので、角度が少しでもずれると、そのずれを拡大する向きにモーメントが生じます。不安定平衡です。論文も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$における淀み点の値です。ここに遮蔽判定を1行だけ足します。ある面がもう一方の球の陰に隠れていれば、圧力を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):
    """単位球面上のほぼ一様な点群と、パネル1枚の面積。"""
    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近似は衝撃波同士の干渉を知りません。論文の図で2つの球の内側の面に現れる高圧の斑点は、2つの弓形衝撃波が出会って作られたものですが、このモデルにはそうした場所がありません。斜め衝撃波列で見たように、衝撃波が重なる場所の圧力は単純な和より大きくなります。つまりこのコードが当てたのは大きさではなく、符号と位置です。

表の中で符号が変わる場所も読んでおく価値があります。$C_L$は150°で0.215と最大になりますが、$C_M$はそれより手前の143.6°で0を通過します。揚力が最大になる姿勢と、自ら維持される姿勢は一致しません。

Q4. 離れた後で横に押す力はどこから来るのか#

接触が切れると話は変わります。2つの球が最初から横並びだった場合(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倍になります。このスケールで無次元化すると、相対運動方程式には係数が1つだけ残ります。

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

$3/8$は球の体積と断面積の比から出る値で、密度比はスケールの中にすべて吸収されています。下のアニメーションはこの方程式をそのまま積分します。

release angleを90°から上へ引っ張ってみましょう。90°付近では実線と点線が重なっています。共通衝撃波が仕事をすべてこなし、サーフィンが付け加えるものがないという意味です。135°を超えると2つの曲線が分かれ、球の縁はピンクの衝撃波線をくわえている間だけ明るくなります。

接触区間の挙動もこのモデルに入っています。165°で放すと2つの球は3–4$\tau_s$の間くっついたまま転がり、整列角が120°付近まで下がってからようやく離れます。論文が報告した接触継続4–6$\tau_s$、接触解除角130°付近とおおむね一致します。後流で抗力が急減する仕方は、抗力危機で見た単一球の抗力曲線とは別の問題です。ここでは剥離点ではなく遮蔽が原因です。

Q5. では180°でなぜ何も起きないのか#

正確に整列したタンデムは対称です。横方向の力も0、モーメントも0です。計算は8$\tau_s$の間、2つの球がくっついたまま姿勢をまったく変えないまま終わります。論文の観測も同じです。

しかしこれは鉛筆を立てておくのと同じ平衡です。上流球の抗力が下流球より大きいので(上の表で0.920対0.000)、少し傾いただけでその傾きを拡大するモーメントが付きます。実際の計算では172.5°で放しても、初めはほとんど動かず次第に速くなる「転がる」運動が現れ、最終的に接触が切れます。

ここで実務的に重要な値が1つ出てきます。接触を保ったまま揚力を出したペアは、質量中心そのものが横へ押されます。論文は150°・157.5°・172.5°で放したペアの質量中心の横方向速度がそれぞれ0.42・0.32・0.28で、個々の破片の相対速度より大きいと記しています。破片1つ1つの反発だけを数えると、この成分は丸ごと抜け落ちます。

13個まで行って残るもの#

論文は4個の四面体配列38通り、13個の面心立方配列34通りで同じ調査を続けます。個々の球の最終的な横方向速度は初期の極角1つでかなりよく整理でき、この傾向は4個でも13個でも同様に現れます。ただしクラスター全体の鈍さがバルク挙動に及ぼす影響は、個数が増えるほど薄れていきます。破片が多くなるほど分離が均質になるということであり、debris-cloudアプローチがどこから使えるようになるのかのヒントになります。

まとめると、破片ペアの分離を支配するのは反発力の大きさではありません。接触がいつ切れるかであり、それを決めるのはモーメントの符号です。そしてその符号は格子なしでもわかります。表面圧力をNewton近似で置き、遮蔽さえきちんと数えれば、2つの球が自ら向かう角度が143.6°だという答えが30行の中で出てきます。