[논문 리뷰] 앞에 구 하나를 내민 사면체는 흩어지지 않았다

파편 하나의 횡속도는 클러스터 안의 초기 극각이 거의 정하고, 클러스터 전체의 분리 속도는 그 극각들을 한 숫자로 접은 뭉툭함 B가 정한다.

앞에 구 하나를 내민 사면체는 흩어지지 않았다#

같은 크기의 철구 네 개를 사면체로 붙여 Mach 20 기류에 놓는다. 꼭지점 하나가 정면을 향하게 두고 손을 놓는다. 직관은 "곧 터지듯 흩어진다"고 말한다. Whalen–Deiterding–Laurence(JFM 2026)의 계산에서 이 클러스터는 거의 흩어지지 않았다. 뒤쪽 세 구를 바깥으로 미는 압력을 접촉력이 정확히 상쇄했기 때문이다.

이 글은 그 논문의 4구 사면체 절(§4)을 읽는다. 파편 하나의 횡속도를 무엇이 정하는지, 클러스터 전체는 어떤 숫자 하나로 요약되는지를 다룬다. 논문이 정의한 뭉툭함 지수 $B$ 와 비대칭 지수 $A$ 를 Python으로 다시 계산한다. 무작위 자세 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$ 이다. $V_T' = 0.35$ 는 20 km/s 유성체라면 횡방향 70 m/s다.

횡속도 $V_{T,i}$ 는 질량중심 기준 상대속도를 $y$ – $z$ 평면의 반지름 방향에 투영한 값이다. 클러스터의 집단 분리 속도 $V_T$ 는 네 구의 평균이다. 유동은 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$ 는 세 구 중심을 지나는 바깥 면의 단위 법선이다. 유동은 $+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만 개와 논문식 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$ 는 결국 네 구의 극각 분포를 한 숫자로 접은 것이다. 면이 정면을 보면 세 구가 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'$ 는 1차 분리 후 0.474, 최종 0.551이다. 비대칭 지수는 $10^{-7}$ 수준으로 사라져 쓸모가 없다.

결론은 자세 효과가 개체 수와 함께 약해진다는 것이다. 13개도 작은 수다. 더 많은 파편이면 자세 의존성은 거의 사라질 것으로 저자들은 본다. 사면체에서 최대 횡속도 0.82는 양력을 내는 세 구 묶음에서 나왔다. 평균의 약 2.5배다.

낙하 분포 모델에 이 숫자를 옮길 때#

파편 수가 적으면 자세가 결과를 가른다. 이때는 평균 하나보다 $B$ 에 따른 분포를 쓰는 편이 낫다. 질량중심 자체도 약 0.13만큼 옆으로 밀린다. 에너지 투하 위치가 바뀐다는 뜻이다.

파편 하나의 속도는 극각 곡선으로 근사할 수 있다. 대칭 배치에서는 접촉 구속이 따로 필요하다. 방향은 반지름 방향 가정이 대체로 버틴다. 접촉면을 한 칸에 세우는 재구성은 더 작은 진폭을 고르는 MUSCL-THINC-BVD의 문제다. 13개를 넘는 클러스터라면 자세 대신 개체 수와 평균 하나로 충분할 수 있다.

[논문 리뷰] 앞에 구 하나를 내민 사면체는 흩어지지 않았다 — Mach 20 파편 클러스터의 뭉툭함 지수