[논문 리뷰] 두 파편은 밀어내지 않고 붙어서 날았다 — Mach 20 구 쌍의 144° 평형각

운석이 대기권에서 두 조각으로 깨졌다고 하자. 두 조각은 서로 멀어질까? 답이 "항상 그렇다"였다면 지상의 크레이터 밭을 역산하는 일은 훨씬 쉬웠을 것이다. Whalen·Deiterding·Laurence가 JFM 2026(vol. 1029, A35)에 실은 Mach 20 계산은 그 답을 흔든다. 어떤 정렬각에서는 두 구가 붙은 채로 함께 날아가고, 심지어 쌍 전체가 양력을 낸다. 이 글은 그 각도가 어디서 오는지를 30줄짜리 Newton 근사로 재현한다. 나온 값은 143.6°였고, 논문의 안정 구간 132°–145.7° 안에 들어 있었다.

θ_0를 옮기고 release를 누르면 굽은 화살표가 초록(벌어짐)과 빨강(닫힘)으로 갈린다.

오른쪽 모멘트 곡선이 0을 지나는 각도로 초기값이 모인다. 그 교차가 143° 부근인지 먼저 본다.

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이 Euler 방정식을 embedded boundary로 풀고, 구는 DYNA3D가 접촉까지 포함해 푼다. 기체는 $\gamma=1.4$ 완전기체, 유입 Mach 20, 구와 기체의 밀도비는 $10^4$이다. 격자는 밀도 구배를 따라 자동으로 세밀해진다. 세밀 패치가 활꼴 충격파를 쫓으므로, 패치 부피가 곧 비용이다.

83번 계산이 쫓는 것도 그 영점 하나다.

Q2. 붙어 있는 쌍은 왜 하필 한 각도로 가는가#

정렬각 $\theta$는 자유류 방향과 "이차 구 → 일차 구" 선분 사이의 각으로 정의한다. $\theta = 180^\circ$면 뒤쪽 구가 앞쪽 구의 후류에 정확히 들어앉는다. $\theta = 90^\circ$면 나란히 선다.

두 구가 접촉을 유지하는 동안은 하나의 강체로 봐도 된다. 이때 접촉력은 내력이라 쌍 전체의 모멘트에서 사라진다. 남는 것은 공력 모멘트뿐이다.

$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

release angle을 90°에서 올리면 실선 횡속도와 점선이 갈라지며 서핑 구간이 생긴다.

90° 근처에서는 두 곡선이 겹친다. 135°를 넘기면 구 테두리가 분홍 충격파를 무는 동안만 밝아진다.

143.61°가 나왔다. 논문이 계수·모멘트 이력을 뒤져 찾아낸 안정 구간은 132°–145.7°이고, 접촉 쌍의 양항비 최대점은 141.5°다. 격자도 없고 충격파도 없는 30줄이 그 안에 들어갔다.

양항비는 0.119로, 논문의 안정 구간 평균 0.197과 최대 0.22보다 낮다. 차이가 나는 이유는 분명하다. Newton 근사는 충격파-충격파 간섭을 모른다. 논문 그림에서 두 구 사이 안쪽 면에 생기는 고압 반점은 두 활꼴 충격파가 만나 만든 것인데, 이 모형에는 그런 자리가 없다. 충격파가 겹치는 자리의 압력은 단순 합보다 크다. 그 중첩을 몇 칸 안에 붙잡는 일은 TENO-THINC 재구성의 몫이다. 그러니 이 코드가 맞힌 것은 크기가 아니라 부호와 위치다.

표에서 부호가 바뀌는 곳도 읽어둘 만하다. $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$은 구의 부피와 단면적 비에서 나온 값이고, 밀도비는 척도 안으로 전부 흡수됐다. 위 적분의 90° 근처는 실선과 점선이 겹친다. 공통 충격파가 일을 다 하고, 서핑이 보탤 횡속도는 없다.

접촉 구간의 거동도 이 모형에 들어 있다. 165°에서 놓으면 두 구는 3–4$\tau_s$ 동안 붙어서 구르다가 정렬각이 120° 근처로 내려온 뒤에야 떨어진다. 논문이 보고한 접촉 지속 4–6$\tau_s$, 접촉 해제 각 130° 부근과 대체로 맞는다. 항력이 후류에서 급감하는 방식은 단일 구의 항력 곡선과는 다른 문제다. 180°에서 하류 항력은 0.000, 상류는 0.920이다. 여기서는 박리점이 아니라 가림이 원인이다.

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 접근이 어디서부터 쓸 만해지는지에 대한 힌트다.

분리 속도를 상수 하나로 박아 두는 선택은 표에서 바로 깨진다. 143.61°에서 $C_D = 1.580$, $L/D = 0.119$다. 양력 최대는 150°의 $C_L = 0.215$인데, 그 행의 $C_M$은 이미 $-0.0020$이다. 135° 행은 $C_D = 1.714$, $C_L = 0.128$, $C_M = 0.0007$로 아직 양이다. 1° 스캔의 영점은 그 사이 143.61°에 있고, 기울기는 $-2.09\times 10^{-4}$/deg라 잡음 바닥 $5.0\times 10^{-5}$보다 가파르다. 90°와 105°의 $C_M$은 0.0000으로 찍혀 교차로 세지 않는다. 복원 부호가 버티는 시작은 111.0°다.

붙잡아야 할 각과 양력이 큰 각은 6° 넘게 어긋난다. 180° 탠덤은 상류 $C_D = 0.920$, 하류 $C_D = 0.000$이다. 150°·157.5°·172.5°에서 놓은 쌍의 질량중심 횡속도는 0.42·0.32·0.28로, 파편 상대속도보다 크다. 양항비 0.119는 논문 평균 0.197, 최대 0.22보다 낮다. 크기까지 맞추려면 활꼴 충격파가 만나는 고압을 되돌려야 하고, 부호만 필요하면 $C_{p,\max} = 1.8394$와 가림만으로 132°–145.7° 안에 들어온다.

90° 행은 $C_D = 1.839$, $C_L = 0.000$, 두 구가 각각 0.920이다. 120°는 $C_D = 1.819$, $C_L = 0.037$, $C_M = 0.0003$으로 바닥 $5.0\times 10^{-5}$ 근처라 영점으로 안 센다. 165°는 $C_D = 1.112$, $C_L = 0.170$, $C_M = -0.0122$, $L/D = 0.152$다. 안정각의 $L/D = 0.119$보다 165°의 양항비가 더 큰데, 모멘트 부호는 이미 음이라 그 자세는 머물지 못한다. $C_{p,\max} = 1.8394$를 90°의 $C_D = 1.839$가 거의 그대로 받아 간다. 가림이 없는 옆나란 자세의 항력 합이 정체 압력계수와 같다는 뜻이다.