[논문 리뷰] 인장에서 0.333, 압축에서 0 — SPH 입자가 짝을 짓는 자리

SPH 인장 불안정을 정하는 것은 응력의 크기가 아니라 응력 부호와 커널 2차 도함수의 곱이다. 커널을 초기 배치에 고정하면 성장률이 0이 된다.

늘였을 때만 입자가 뭉쳤다#

같은 재료를 같은 코드로 두 번 계산했다. 한 번은 눌렀고 한 번은 당겼다. 누른 계산은 끝까지 돌았다. 당긴 계산은 입자가 둘씩 달라붙고 그 사이에 구멍이 생기며 죽었다. 이 글은 그 비대칭이 어디서 오는지를 고유값 하나로 재고, 해법으로 제안된 것들이 각각 무엇을 대가로 내는지까지 정리한다. 읽고 나면 SPH 코드의 인장 계산이 터질 때 응력 크기 대신 무엇을 먼저 볼지 알게 된다.

소재는 Xu 등이 Acta Mechanica Sinica(2023)에 실은 SPH 리뷰다. 유체·고체·생체역학을 훑는 글인데, 고체 파트가 인장 불안정을 "도전 과제"로 따로 떼어 놓고 응력점법, 보존 스무딩, 인공응력, 그리고 Total Lagrangian SPH를 나란히 올려 둔다. 왜 해법이 네 개나 필요한지, 직접 재봤다.

커널이 이웃을 다시 세는 순간 생기는 되먹임#

SPH의 운동방정식은 이렇게 생겼다. 응력 $\sigma$ 는 인장을 양수로 잡는다.

\frac{dv_i}{dt}=\sum_j m_j\left(\frac{\sigma_i}{\rho_i^{2}}+\frac{\sigma_j}{\rho_j^{2}}\right) \frac{\partial W_{ij}}{\partial x_i}

$m_j$ 는 입자 질량, $\rho_i$ 는 밀도, $W_{ij}=W(|x_i-x_j|/h)$ 는 매끄럼 길이 $h$ 를 쓴 커널이다. 밀도는 격자에서 오지 않고 이웃을 세서 나온다.

\rho_i=\sum_j m_j\,W(|x_i-x_j|/h)

여기에 함정이 하나 숨어 있다. 입자가 움직이면 $W_{ij}$ 도, $\partial W_{ij}/\partial x_i$ 도 같이 바뀐다. 즉 변형이 힘을 만들고, 그 힘이 다시 힘을 계산하는 도구를 바꾼다. 유한요소법에서는 형상함수가 절점 배치에 붙어 있어 이런 되먹임이 없다. SPH는 커널을 매 스텝 현재 위치에서 다시 짓기 때문에 있다.

sigma0 / E 슬라이더를 양수로 밀면 파장 $2\Delta x$ 씨앗이 자라고, 음수로 두면 로그 축이 눕는다.

sigma0 / E = +0.05 h / dx = 1.20  freeze kernel on initial lattice (TL-SPH)

48 particles on a periodic line, seeded with a 2dx zig-zag of amplitude 1e-4 dx and nothing else. Push sigma0 positive and the seed grows until neighbours collide into visible pairs; pull it negative and the same seed just sits there. The trace is logarithmic, so a straight rising line is exponential growth. Tick the checkbox to freeze the kernel on the initial lattice — the line goes flat at every sigma0, step 0 included.

sigma0 / E를 양수로 밀면 파장 $2\Delta x$ 지그재그 씨앗이 지수로 자라 입자가 눈에 보이게 짝을 짓는다. 같은 씨앗을 음수 응력에 두면 로그 축의 선이 평평하게 눕는다. 체크박스로 커널을 초기 배치에 고정하면 응력을 어떻게 놓아도 선이 눕는다는 점까지 확인할 수 있다.

부호 하나로 갈리는 조건#

Swegle이 1995년에 정리한 결론은 응력의 크기가 아니라 부호 조합에 관한 것이다. 균일하게 응력을 받은 입자 배치에 파수 $k$ 의 미소 변위를 얹고 선형화하면, 진동수는 두 항의 차로 나온다.

\rho_0\,\omega^{2}(k)=E\,A(k)-\sigma_0\,B(k)

$A(k)$ 는 탄성 강성에서 오는 항으로 항상 0 이상이다. $B(k)$ 는 커널을 다시 짓는 데서 오는 항이고, 이웃마다 $W''(r_{ij})$ 를 모아 놓은 합이다. $\omega^2$ 가 음수가 되면 $\omega$ 가 순허수가 되고 변위는 지수로 자란다. 조건을 이웃 하나 단위로 줄여 쓰면 유명한 형태가 된다.

\sigma\,W''(r_{ij})\gt 0 \;\Longrightarrow\; \text{unstable}

여기서 비대칭이 설명된다. 3차 B-스플라인 커널의 $W''$ 는 $r/h$ 가 약 $2/3$ 을 넘으면 부호를 바꾼다. 인장( $\sigma\gt 0$ )에서는 먼 이웃이 위험하고, 압축( $\sigma\lt 0$ )에서는 가까운 이웃이 위험하다. 실제 배치에서 이웃은 대부분 $r/h\gt 2/3$ 쪽에 앉으므로, 위험은 인장에 몰린다. 3차 B-스플라인의 $W''$ 가 $r/h \approx 2/3$ 에서 부호를 바꾸는 자리는 커널과 인공점성부터 적는 SPH 입자법에서 커널을 고른 다음에 오는 문제다.

성장률을 고유값으로 직접 재보기#

부호 조건은 이웃 하나 기준이다. 실제 코드에서 알아야 할 것은 배치 전체에서 가장 빨리 자라는 모드의 성장률이다. 시간적분으로 재면 격자 잡음과 시간오차가 섞이니, 가속도 연산자의 강성행렬을 수치미분으로 뽑아 고유값을 보는 편이 깨끗하다. 주기 경계의 1차원 입자열, $E=\rho_0=1$ , $\Delta x=1$ 로 두면 성장률은 곧 $c/\Delta x$ 단위가 된다.

import numpy as np
 
def spline_w(q):
    w = np.zeros_like(q)
    a, b = q < 1.0, (q >= 1.0) & (q < 2.0)
    w[a] = 1.0 - 1.5 * q[a] ** 2 + 0.75 * q[a] ** 3
    w[b] = 0.25 * (2.0 - q[b]) ** 3
    return (2.0 / 3.0) * w
 
def spline_dw(q):
    d = np.zeros_like(q)
    a, b = q < 1.0, (q >= 1.0) & (q < 2.0)
    d[a] = -3.0 * q[a] + 2.25 * q[a] ** 2
    d[b] = -0.75 * (2.0 - q[b]) ** 2
    return (2.0 / 3.0) * d
 
def pair_tables(x, h, L):
    d = x[:, None] - x[None, :]
    d -= L * np.round(d / L)                 # 주기 경계에서의 최단 간격
    q = np.abs(d) / h
    dW = spline_dw(q) / h ** 2 * np.sign(d)  # dW/dx_i
    np.fill_diagonal(dW, 0.0)
    return spline_w(q) / h, dW
 
def eulerian_accel(u, x0, h, L, m, rho0, E, sig0):
    """고전 SPH: 커널을 변형된 현재 위치에서 다시 짓는다."""
    W, dW = pair_tables(x0 + u, h, L)
    rho = m * W.sum(axis=1)
    sig = sig0 + E * (rho0 / rho - 1.0)
    return (m * ((sig / rho ** 2)[:, None] + (sig / rho ** 2)[None, :]) * dW).sum(axis=1)
 
def lagrangian_accel(u, x0, h, L, m, rho0, E, sig0):
    """Total Lagrangian SPH: 커널을 초기 배치에 고정한다."""
    _, dW0 = pair_tables(x0, h, L)
    F = 1.0 + ((m / rho0) * (u[None, :] - u[:, None]) * dW0).sum(axis=1)
    P = sig0 + E * (F - 1.0)
    return (m / rho0 ** 2 * (P[:, None] + P[None, :]) * dW0).sum(axis=1)
 
def worst_mode(accel, x0, h, L, m, rho0, E, sig0, eps=1e-7):
    """강성행렬을 수치미분으로 뽑아 가장 빨리 자라는 변위 모드를 찾는다."""
    n = len(x0)
    K = np.zeros((n, n))
    for j in range(n):
        p = np.zeros(n)
        p[j] = eps
        K[:, j] = (accel(p, x0, h, L, m, rho0, E, sig0)
                   - accel(-p, x0, h, L, m, rho0, E, sig0)) / (2.0 * eps)
    lam, V = np.linalg.eig(K)                # u'' = K u 이므로 omega^2 = -lam
    i = int(np.argmax(lam.real))
    rate = np.sqrt(max(lam[i].real, 0.0))    # c/dx 단위의 성장률
    spec = np.abs(np.fft.rfft(V[:, i].real))
    kmode = int(np.argmax(spec[1:])) + 1
    return rate, L / kmode
 
n, dx, rho0, E = 64, 1.0, 1.0, 1.0
x0 = np.arange(n) * dx
L, m = n * dx, rho0 * dx
 
print("sigma0/E   h/dx   Eulerian rate   worst lambda/dx   frozen-kernel rate")
for sig in (+0.05, 0.0, -0.05):
    for hf in (1.0, 1.2, 1.5):
        ge, lam = worst_mode(eulerian_accel, x0, hf * dx, L, m, rho0, E, sig * E)
        gl, _ = worst_mode(lagrangian_accel, x0, hf * dx, L, m, rho0, E, sig * E)
        tag = f"{lam:.2f}" if ge > 1e-4 else "none"       # 자라는 모드가 없으면 파장도 없다
        print(f"{sig:+8.2f}   {hf:4.2f}   {ge:13.6f}   {tag:>15}   {gl:18.6f}")
 
print("\nwidening the kernel at sigma0/E = +0.05")
print(" h/dx   rate   worst lambda/dx")
for hf in (1.0, 1.2, 1.4, 1.5, 1.7, 2.0):
    ge, lam = worst_mode(eulerian_accel, x0, hf * dx, L, m, rho0, E, +0.05 * E)
    print(f"{hf:5.2f}   {ge:6.4f}   {lam:13.2f}")
 
print("\nrate scaling with the pre-stress at h/dx = 1.2")
for sig in (0.01, 0.02, 0.05, 0.10, 0.20):
    ge, _ = worst_mode(eulerian_accel, x0, 1.2 * dx, L, m, rho0, E, sig * E)
    print(f"sigma0/E = {sig:+5.2f}   rate = {ge:6.4f}   rate/sqrt(sigma0/E) = {ge / np.sqrt(sig):6.4f}")

sigma0/E   h/dx   Eulerian rate   worst lambda/dx   frozen-kernel rate
   +0.05   1.00        0.632456              2.00             0.000000
   +0.05   1.20        0.333437              2.00             0.000012
   +0.05   1.50        0.048837              2.29             0.000007
   +0.00   1.00        0.000000              none             0.000000
   +0.00   1.20        0.000000              none             0.000000
   +0.00   1.50        0.000000              none             0.000000
   -0.05   1.00        0.000000              none             0.000000
   -0.05   1.20        0.000004              none             0.000012
   -0.05   1.50        0.000000              none             0.000000
 
widening the kernel at sigma0/E = +0.05
 h/dx   rate   worst lambda/dx
 1.00   0.6325            2.00
 1.20   0.3334            2.00
 1.40   0.1372            2.00
 1.50   0.0488            2.29
 1.70   0.0920            2.56
 2.00   0.1335            3.20
 
rate scaling with the pre-stress at h/dx = 1.2
sigma0/E = +0.01   rate = 0.1375   rate/sqrt(sigma0/E) = 1.3755
sigma0/E = +0.02   rate = 0.2049   rate/sqrt(sigma0/E) = 1.4489
sigma0/E = +0.05   rate = 0.3334   rate/sqrt(sigma0/E) = 1.4912
sigma0/E = +0.10   rate = 0.4759   rate/sqrt(sigma0/E) = 1.5050
sigma0/E = +0.20   rate = 0.6761   rate/sqrt(sigma0/E) = 1.5119

h / dx 슬라이더를 1.0에서 2.0으로 밀면 이웃이 $r/h = 2/3$ 선을 넘고, $r = 2h$ 에서 새 점이 들어온다.

h / dx = 1.20 sigma0 / E = +0.05  show W'

Vertical lines are the lattice neighbours of one particle; the dot sits on W'' at that separation. Red dot = that pair is destabilising, and its arrow swells inward. Slide h/dx and watch neighbours cross the sign change of W'' near r/h = 2/3, and new ones enter at r = 2h. Flip sigma0 negative and every red dot turns blue without the curve moving at all — only the sign of the product changed (frame 0).

수직선은 한 입자의 격자 이웃이고 점은 그 간격의 $W''$ 다. 곡선을 그대로 둔 채 sigma0 / E만 음수로 뒤집으면 빨간 점이 전부 파랗게 바뀐다. 바뀌는 것은 곱의 부호뿐이다.

가장 빨리 자란 것은 파장 $2\Delta x$ 였다#

첫 표의 인장 세 줄만 성장률이 0이 아니다. $h/\Delta x=1.2$ 에서 0.3334, 압축에서는 0이다. $10^{-5}$ 자리의 값들은 강체 평행이동 모드의 반올림이니 0으로 읽으면 된다.

중요한 것은 옆 칸이다. 가장 빨리 자라는 모드의 파장이 $2.00\,\Delta x$ 다. 격자에서 표현할 수 있는 가장 짧은 파장, 즉 입자가 하나씩 번갈아 좌우로 밀리는 모드다. 그래서 화면에서 관찰되는 증상이 "입자쌍"이다. 두 입자가 서로에게 다가가고, 그 옆에는 아무도 없는 구멍이 남는다. 밀도 합산은 구멍에서 $\rho$ 를 낮게 읽고, 그러면 그 자리 응력이 더 강한 인장이 되고, 되먹임이 닫힌다. 입자 하나 간격이 장 전체를 바꾸는 구조다.

세 번째 표는 성장률이 $\sqrt{\sigma_0}$ 에 붙어 있다고 말한다. 마지막 열의 비가 1.38에서 1.51 사이에 머문다. 인장을 4배 키우면 성장률은 2배다. 응력 1%만 걸어도 성장률이 0.1375라는 뜻이기도 하다. 스텝 수를 늘리면 어차피 터진다.

0.6325 / 0.0488 = 13.0이라 $h/\Delta x$ 를 1.0에서 1.5로 늘릴 때의 13배가 맞다. 2.0의 0.1335는 1.5의 0.0488보다 2.7배 크다. 최악 파장도 2.00에서 3.20으로 길어진다. $\sigma_0/E = 0.20$ 의 0.6761을 0.05의 0.3334로 나누면 2.03이다. $\sqrt{0.20/0.05} = 2$ 라, 응력이 4배일 때 성장률이 2배라는 말과 같다. 마지막 열 1.4912와 1.5119는 붙어 있다. $0.1375 / \sqrt{0.01} = 1.375$ 라 표의 1.3755와 같다. 인장 +0.05, $h/\Delta x = 1.2$ 에서 오일러 성장률 0.333437 옆의 고정 커널은 0.000012다. 2만 배를 넘긴다. $h/\Delta x = 1.0$ 에서는 0.632456 대 0.000000이라 고정 커널 열은 영이다. 압축 $-0.05$ , $h/\Delta x = 1.2$ 의 $4\times10^{-6}$ 은 강체 모드 반올림이지 성장이 아니다.

$h$ 를 넓히면 사라지는가#

두 번째 표가 이 글에서 가장 반직관적인 부분이다. $h/\Delta x$ 를 1.0에서 1.5로 늘리면 성장률이 0.6325에서 0.0488로 13배 줄어든다. 여기까지는 예상대로다. 그런데 1.5를 지나면 다시 오른다. 1.7에서 0.0920, 2.0에서 0.1335다.

파장 열이 이유를 알려준다. 1.4까지는 $2.00\,\Delta x$ 였던 최악 모드가 1.5에서 $2.29\,\Delta x$ , 2.0에서 $3.20\,\Delta x$ 로 길어진다. $h$ 를 넓히면 $r=2h$ 안으로 이웃이 더 들어온다. 새로 들어온 이웃의 $W''$ 부호가 다시 불안정 쪽이면, 불안정은 없어지는 게 아니라 더 긴 파장으로 옮겨 앉는다.

이것이 실무에서 곤란한 이유는 분명하다. 더 긴 파장의 불안정은 눈에 잘 띄지 않는다. 입자가 둘씩 붙는 것은 로그만 봐도 안다. 세 칸짜리 뭉침은 "물리적인 국소화"처럼 보인다. $h$ 를 키우는 것은 증상을 치료하는 게 아니라 눈에 안 띄게 만드는 쪽에 가깝다.

커널을 초기 배치에 묶는다#

표의 마지막 열이 해법이다. Total Lagrangian SPH는 커널과 그 도함수를 초기 배치 $X$ 에서 한 번만 계산하고 끝까지 고정한다. 변형은 변형기울기로 들어온다.

F_i=1+\sum_j \frac{m_j}{\rho_0}\,(u_j-u_i)\,\frac{\partial W^{0}_{ij}}{\partial X_i}, \qquad P_i=\sigma_0+E\,(F_i-1)

$u$ 는 변위, $W^0$ 는 초기 배치에서 만든 커널, $P$ 는 1차 Piola–Kirchhoff 응력이다. 위치가 $W^0$ 안에 없으므로 변형이 커널을 바꾸는 경로가 끊긴다. 성장률이 어느 응력에서도 0으로 나오는 것이 그 결과다. 이산 연산자가 어떤 항등식을 그대로 물려받는지가 안정성을 가른다.

대가도 분명하다. 초기 이웃 목록을 끝까지 쓰므로 균열·파편화·대변형으로 이웃 관계 자체가 바뀌는 문제는 못 푼다. 접촉을 새로 만들 수도 없다. 그래서 고체 충격 해석에서는 TL 정식화를 쓰면서 파괴 모델을 따로 붙이고, 이웃을 갱신하는 시점을 손으로 정한다.

응력점·인공응력· $\delta$ -SPH가 각각 사는 자리#

리뷰가 함께 올린 나머지 처방은 각기 다른 경로를 건드린다.

응력점법은 변위와 응력을 다른 점에 둔다. 파장 $2\Delta x$ 모드가 응력점에서 0이 아닌 변형률을 만들게 되므로, 앞에서 본 최악 모드가 강성을 얻는다. 엇갈린 배치가 체커보드 모드를 잡는 것과 같은 원리다.

인공응력은 인장 영역에 짧은 거리에서만 작동하는 반발 항을 더한다. $\sigma W''\gt 0$ 를 직접 고치지 않고 반대 방향 힘으로 덮는다. 계수 선택이 결과를 바꾸므로 검증 부담이 크다.

보존 스무딩은 변형률 장을 매 스텝 평활화한다. 짧은 파장을 깎으니 최악 모드가 먼저 죽지만, 물리적으로 있어야 할 응력 국소화도 같이 깎인다.

$\delta$ -SPH는 사실 인장 불안정 처방이 아니다. 밀도 장에 확산 항을 넣어 약압축성 유체 계산의 압력 잡음을 줄이는 도구다. 리뷰가 둘을 같은 장에 나란히 놓는 것은 증상이 비슷해 보이기 때문인데, 원인은 다르다. 유체 쪽 잡음은 상태방정식이 밀도 오차를 증폭해서 생긴다. 고체 쪽 짝짓기는 응력 부호와 $W''$ 의 곱에서 온다. 처방을 서로 빌려 쓰면 증상만 가려진다.

입자법을 고를 때 실제로 서명하는 것#

SPH가 격자를 지운 대가는 "커널이 해에 의존한다"는 한 줄로 요약된다. 압축 지배 문제에서는 이 의존성이 공짜처럼 보인다. 인장이 들어오는 순간 청구서가 온다.

코드를 고를 때 물어볼 것은 안정화 옵션의 개수가 아니다. 이 계산에서 인장이 얼마나 오래 유지되는지, 그리고 이웃 관계가 바뀌어야 하는지 둘뿐이다. 인장이 길고 이웃이 그대로면 TL 정식화가 맞다. 이웃이 계속 바뀌어야 하면 갱신형 커널을 쓰고 안정화 계수를 검증 대상으로 올려야 한다. 성장률 0.3334는 옵션을 몇 개 켰는지 묻지 않는다. 이산 연산자가 어떤 항등식을 물려받는지는 컬이 경로마다 0.544 달라진 구속에서도 같은 질문이었다.