[論文レビュー] 引張で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でないのは引張の3行だけです。 $h/\Delta x=1.2$ で0.3334、圧縮では0です。 $10^{-5}$ 桁の値は剛体並進モードの丸め誤差なので、0と読んでよいものです。

大事なのは隣の列です。最も速く育つモードの波長が $2.00\,\Delta x$ です。格子で表現できる最短波長、つまり粒子が一つずつ交互に左右へ押されるモードです。だから画面で観察される症状が「粒子対」になります。二つの粒子が互いに近づき、その隣には誰もいない穴が残ります。密度和は穴で $\rho$ を低く読み、するとその位置の応力がさらに強い引張になり、フィードバックが閉じます。粒子ひとつ分の間隔が場全体を変えます。

3番目の表は、成長率が $\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$ を広げれば消えるのか#

2番目の表がこの記事で最も直観に反する部分です。 $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$ は第1Piola–Kirchhoff応力です。位置が $W^0$ の中にないので、変形がカーネルを書き換える経路が切れます。成長率がどの応力でも0として出てくるのがその結果です。離散作用素がどの恒等式をそのまま受け継ぐかで、安定性が分かれます。

対価も明確です。初期の隣接リストを最後まで使うので、亀裂・破砕・大変形で隣接関係そのものが変わる問題は解けません。接触を新しく作ることもできません。ですから固体衝撃解析ではTL定式化を使いながら破壊モデルを別に付け、隣接を更新する時点を手で決めます。

応力点・人工応力・ $\delta$ -SPHがそれぞれ住む場所#

レビューが一緒に載せた残りの処方は、それぞれ違う経路に触ります。

応力点法は変位と応力を別の点に置きます。波長 $2\Delta x$ のモードが応力点で0でないひずみを作るようになるので、先ほどの最悪モードが剛性を得ます。スタガード配置がチェッカーボードモードを抑えるのと同じ原理です。

人工応力は引張領域に短距離でのみ働く反発項を足します。 $\sigma W''\gt 0$ を直接直さず、逆方向の力で覆います。係数の選び方が結果を変えるので、検証の負担が大きくなります。

保存スムージングはひずみ場を毎ステップ平滑化します。短い波長を削るので最悪モードが先に死にますが、物理的にあるべき応力の局所化も一緒に削られます。

$\delta$ -SPHは実は引張不安定の処方ではありません。密度場に拡散項を入れて弱圧縮性流体計算の圧力ノイズを減らす道具です。レビューが二つを同じ場に並べるのは症状が似て見えるからですが、原因は違います。流体側のノイズは状態方程式が密度誤差を増幅して生じます。固体側の対形成は応力の符号と $W''$ の積から来ます。処方を互いに借りてくると症状だけが隠れます。

粒子法を選ぶとき実際に署名するもの

SPHが格子を消した対価は「カーネルが解に依存する」という一行に要約されます。圧縮支配の問題ではこの依存性は無料のように見えます。引張が入ってきた瞬間に請求書が来ます。

コードを選ぶときに問うのは安定化オプションの数ではありません。この計算で引張がどれだけ長く維持されるか、そして隣接関係が変わる必要があるか、その二つだけです。引張が長く隣接がそのままならTL定式化が合います。隣接が変わり続ける必要があるなら更新型カーネルを使い、安定化係数を検証対象に上げなければなりません。成長率0.3334はオプションを何個入れたかを問いません。離散作用素がどの恒等式を受け継ぐかは、カールが経路ごとに0.544違った拘束でも同じ問いでした。