心臓弁を格子に刺さずに — Peskinの埋め込み境界法と離散デルタ関数

補間と拡散が同じ離散デルタを使わないと、マーカーの力は流体の運動量として戻ってきません。

1972年、Charles Peskin はコロンビア大学医学部で心臓弁の格子に行き詰まりました。弁は拍動のたびに曲がり、膨らみ、閉じます。その動きに合わせて毎ステップ格子を切り直すのは持てません。格子はそのままにして、物体は点として流体へ力を流します。

マーカーを格子の間へ動かすスライダーを送ると、4点の重みが敷き直され、和は 1 の近くに残ります。

Marker position x_L = 10.30 h

Each bar is the weight δ_h(x_i − x_L) the Lagrangian marker hands to grid node i. The sum stays close to 1 by design — that is the zeroth moment condition.

幅を 2 点ハットから 6 点 Peskin へ広げると、力はより広く散らばります。マーカーが格子線をまたぐとき、和が 1 から跳ねないかを見ます。

1972年、Peskinが向き合った心臓弁#

心臓シミュレーションで最も厄介なのは格子でした。弁が動くと ALE(Arbitrary Lagrangian–Eulerian)で格子も追従させる必要がありますが、動きが大きく速すぎます。Peskin の発想は二つの座標系を完全に分離することでした。流体は固定した直交格子(Eulerian)で解き、物体は動く点群(Lagrangian)で表す。二つの座標系の橋渡しはディラックのデルタ関数 $\delta$ が行います。

一行で書くと、

\mathbf{f}(\mathbf{x}, t) = \int_\Gamma \mathbf{F}(s, t)\, \delta(\mathbf{x} - \mathbf{X}(s, t))\, ds

ここで $\mathbf{x}$ は流体格子点、 $\mathbf{X}(s,t)$ はラグランジアンマーカー、 $\mathbf{F}$ はマーカーが流体に与える力です。弁の形も動きも、すべて $\mathbf{X}(s,t)$ という一つの関数に入ります。

メッシュを切り直さない

問題はコンピュータです。コンピュータはディラックのデルタを直接扱えません。格子間隔 $h$ より狭い領域で無限大に立ち上がる関数は離散化できません。Peskin が導入したのは離散デルタ関数 $\delta_h$ です。ディラックを幅 $h$ の滑らかな山に置き換えたものです。

二つの条件を要求します。

全格子で重みの和がちょうど 1(0 次モーメント)
重みと座標の積の和がマーカー座標に一致(1 次モーメント)

和が 1 は加える力が保存することを意味します。1 次モーメントはトルクがずれないことの保証です。この二条件だけで 4 点カーネルが一意に決まります。

離散デルタ — 一点を格子に広げる

2 点ハットから 6 点 Peskin カーネルまで、サポート幅が広くなるほど重みは滑らかに散らばります。滑らかなほど、マーカーが格子線をまたぐ瞬間に力が飛び跳ねません。これが IBM の時間振動が小さくなる理由です。4 点 Peskin カーネルが最もよく使われるのは、上記の二条件に加えて「重みの二乗和がマーカー位置に依らない」という第三の条件まで満たすからです。一度広げた力が滑るなら、同じカーネルで何度も強制し直します。その繰り返しは格子が壁を知らない側のデルタカーネルにあります。

閉形式では、

\delta_h(r) = \begin{cases} \dfrac{1}{8}\left(3 - 2|r| + \sqrt{1 + 4|r| - 4r^2}\right), & 0 \le |r| < 1 \\ \dfrac{1}{8}\left(5 - 2|r| - \sqrt{-7 + 12|r| - 4r^2}\right), & 1 \le |r| < 2 \\ 0, & |r| \ge 2 \end{cases}

ただし $r = (x_i - x_L)/h$ 、 $x_i$ は格子点、 $x_L$ はマーカー位置です。

直接強制 — 二つの矢印の往復

二次元に移ります。円柱表面にマーカーを並べ、それぞれが no-slip 条件を流体に課すとします。各時間ステップで以下の四段階を繰り返します。

補間(E → L):格子速度 $\mathbf{u}$ をマーカー位置に引き寄せる。

\mathbf{U}_l = \sum_{i,j} \mathbf{u}_{ij}\, \delta_h(\mathbf{x}_{ij} - \mathbf{X}_l)\, h^2

力の計算:マーカーで望む速度 $\mathbf{U}_l^{des}$ との差を時間スケールで割る。

\mathbf{F}_l = \frac{\mathbf{U}_l^{des} - \mathbf{U}_l}{\Delta t}

拡散(L → E):その力を格子に逆向きに撒く。

\mathbf{f}_{ij} = \sum_l \mathbf{F}_l\, \delta_h(\mathbf{x}_{ij} - \mathbf{X}_l)\, \Delta s_l

流体更新: $\mathbf{u}^{n+1} = \mathbf{u}^n + \Delta t\,(\dots + \mathbf{f}/\rho)$ 。

補間と拡散は同じカーネル $\delta_h$ を使います。この対称性で、マーカーの力が流体の運動量として戻ります。

Markers around the circle: 24

Orange dots are Lagrangian markers traced along the circle. The cyan halo is the total spreading weight Σ_l δ_h(x − X_l) on each Eulerian cell. Increase markers until the band is uniform — that is the rule of thumb Δs < h.

マーカーが少なすぎると円周に隙間ができます。経験則は $\Delta s_l < h$ 、マーカー間隔が格子より細かいこと。64 個まで上げるとバンドが均一になります。

100 行 IBM — 円柱まわり一ステップ#

import numpy as np
 
def peskin4(r):
    # 4点正則化デルタ (Peskin 2002)
    a = np.abs(r)
    w = np.zeros_like(r)
    m1 = a < 1
    m2 = (a >= 1) & (a < 2)
    w[m1] = (3 - 2*a[m1] + np.sqrt(1 + 4*a[m1] - 4*a[m1]**2)) / 8
    w[m2] = (5 - 2*a[m2] - np.sqrt(-7 + 12*a[m2] - 4*a[m2]**2)) / 8
    return w
 
def spread_force(F_l, X_l, ds, shape, h):
    # ラグランジアン力 F_l -> オイラー格子力 f_ij
    nx, ny = shape
    f = np.zeros((nx, ny, 2))
    for l in range(len(X_l)):
        i0, j0 = int(X_l[l, 0] / h), int(X_l[l, 1] / h)
        for di in range(-2, 3):
            for dj in range(-2, 3):
                ii, jj = i0 + di, j0 + dj
                if not (0 <= ii < nx and 0 <= jj < ny):
                    continue
                wx = peskin4(np.array([ii - X_l[l, 0] / h]))[0]
                wy = peskin4(np.array([jj - X_l[l, 1] / h]))[0]
                f[ii, jj] += F_l[l] * wx * wy * ds[l]
    return f
 
def interp_velocity(u, X_l, h):
    # 格子速度 u_ij -> マーカー速度 U_l
    nx, ny, _ = u.shape
    U = np.zeros((len(X_l), 2))
    for l in range(len(X_l)):
        i0, j0 = int(X_l[l, 0] / h), int(X_l[l, 1] / h)
        for di in range(-2, 3):
            for dj in range(-2, 3):
                ii, jj = i0 + di, j0 + dj
                if not (0 <= ii < nx and 0 <= jj < ny):
                    continue
                wx = peskin4(np.array([ii - X_l[l, 0] / h]))[0]
                wy = peskin4(np.array([jj - X_l[l, 1] / h]))[0]
                U[l] += u[ii, jj] * wx * wy * h**2
    return U
 
# direct-forcing IBM の1ステップ
N, h, dt = 32, 1.0, 0.05
nL = 24
theta = np.linspace(0, 2*np.pi, nL, endpoint=False)
R = 8.0
X_l = np.stack([N/2 + R*np.cos(theta), N/2 + R*np.sin(theta)], axis=1)
ds = np.full(nL, 2*np.pi*R / nL)
 
u = np.zeros((N, N, 2))
u[:, :, 0] = 1.0                      # 一様流入速度
U_des = np.zeros((nL, 2))             # 静止円柱 = no-slip
U_l = interp_velocity(u, X_l, h)
F_l = (U_des - U_l) / dt
f = spread_force(F_l, X_l, ds, (N, N), h)
u_new = u + dt * f / 1.0
print(f"max |u| inside cylinder: {np.abs(u_new[12:20, 12:20]).max():.4f}")

一ステップだけで円柱内部の速度がほぼゼロまで落ちます。実際のシミュレーションでは圧力ポアソンと粘性項が要りますが、IBM の核心 — 補間・計算・拡散の三拍子 — はこの 30 行に入っています。

境界がぼやける代償

Peskin のオリジナル IBM には対価があります。 $\delta_h$ で一度伸ばすと境界は幅 $h$ の拡散界面になります。マーカー位置は正確でも、流体が見る壁はわずかにぼやけます。低 Re なら問題ありません。境界層が厚いので埋もれます。しかし Re が大きくなると境界層が格子一マスより薄くなり、ぼやけた壁の上では BL を解けません。

圧力にも振動が出ます。マーカーが格子線をまたぐと強制力は滑らかに変わりますが、圧力ポアソンは微小なリップルを増幅します。剛体の場合 implicit IBM で振動は抑えられますが、毎時刻線形システムを一回余分に解く代価が付きます。

Sharp-interface 系が登場した理由#

この制限から、1990 年代後半以降、不連続強制(discrete forcing)IBM 系が登場しました。Mittal・Iaccarino の ghost-cell IBM、Tseng・Ferziger の cut-cell、Tucker の sharp-interface 系。共通の発想は一つです。 $\delta_h$ を捨て、境界近傍のセルで直接補間問題を解いて境界条件を正確に課す。

手法	長所	短所
Peskin 連続強制	単純、保存的、変形体も扱える	境界拡散、圧力振動
Ghost-cell IBM	境界鮮明、高次化可能	コード複雑、変形体困難
Cut-cell	幾何が正確	小さなセルの安定性問題
MLS IBM	非一様格子と相性が良い	重み計算が高コスト

心臓弁のようなしなやかな膜は今も Peskin 系が自然で、剛体運動 + 高 Re 外部流は ghost-cell が優勢です。両者は競争ではなく分業です。

節点上の重みは 1/2 と 1/4 で、マーカー 24 個は 2h を超える#

マーカーが格子点にちょうど乗ると、残る距離は $r = 0$ , $\pm 1$ , $\pm 2$ だけです。 $r = 0$ では前半の式が $(3 + \sqrt{1})/8 = 1/2$ です。 $|r| = 1$ では二つの式がつながり、どちらも $(3 - 2 + \sqrt{1})/8 = 1/4$ です。 $|r| = 2$ は 0 です。並べると $0$ , $1/4$ , $1/2$ , $1/4$ , $0$ で、和は 1 です。0 次モーメントはこの五つの節点で閉じます。

コードの数字は nL = 24, R = 8, h = 1 です。円周は $2\pi R = 16\pi \approx 50.27$ です。間隔 $2\pi R/24 \approx 2.09$ は格子幅の 2 倍を超えます。経験則 $\Delta s_l < h$ が成り立つにはマーカーが 51 個を超える必要があります。64 個なら $50.27/64 \approx 0.79$ で、そこで初めて h より狭くなります。24 個のリストは補間と拡散を一歩見せるだけで、閉じた壁ではありません。

リストのマーカー 24 個は間隔の規則をすでに超えている#

間隔が約 $2.09h$ の円周には、格子の間に空いた弧が残ります。64 個で帯が均一になるという話と、リストの 24 個は同じ壁ではありません。一歩の出力だけを見て no-slip が閉じたと読むと、粗いマーカーが残した穴を見落します。

弁は点の鎖である前に薄い膜です。薄く切るほどせん断が偽って固まる側はシェル要素のせん断ロッキングが別に示しています。流体側の重み和を 1 に合わせても、膜がロックすれば弁の角度は別の理由でずれます。