不给心脏瓣膜画贴体网格 — Peskin 浸入边界法与离散 delta 函数

插值和扩散必须共用同一个离散 delta，标记上的力才会回到流体动量里。

1972 年,Charles Peskin 在哥伦比亚大学医学院被心脏瓣膜的网格卡住了。瓣膜每次搏动都在弯、胀、合上。每一步都重划网格,承担不了。网格不动。物体做成一串点,把力漏进流体。

把标记在格子之间拖动,四点权重会重新铺开,和仍停在 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.

支撑从两点帽函数加宽到六点 Peskin 核,力铺得更开。标记跨过格线时,看权重和会不会从 1 上跳开。

1972 年,Peskin 面对心脏瓣膜#

心脏模拟里最棘手的就是网格。瓣膜一动,ALE(Arbitrary Lagrangian–Eulerian)就要求网格跟着动,可这个运动幅度太大、太快了。Peskin 的想法是把两个坐标系彻底解耦。流体在固定的直角网格(Eulerian)上求解,物体用一组移动的点(Lagrangian)表示。两个坐标系之间的桥梁,就是 Dirac delta 函数 $\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)$ 是 Lagrangian 标记点, $\mathbf{F}$ 是标记点施加给流体的力。瓣膜的形状和运动全装在 $\mathbf{X}(s,t)$ 这一个函数里。

不重新划网格

问题在计算机。计算机没法直接处理 Dirac delta。在小于网格间距 $h$ 的区域里冲到无穷的函数没法离散化。Peskin 引入的是离散 delta 函数 $\delta_h$ —— 把 Dirac 换成宽度为 $h$ 的光滑小山。

它要满足两个条件。

所有格点上权重之和精确等于 1(零阶矩)
权重与坐标乘积之和等于标记点坐标(一阶矩)

权重和为 1 表示施加的力守恒。一阶矩保证力矩不漂移。这两个条件就唯一确定了 4 点核。

离散 delta —— 把一点铺到网格上#

从 2 点 hat 到 6 点 Peskin 核,支撑越宽,权重铺得越平滑。光滑度越高,标记点跨过格线时受力越平稳,这就是 IBM 时间振荡更小的原因。4 点 Peskin 核用得最广,因为它除了上面两条还满足第三条:权重平方和与标记点位置无关。一次铺开的力若还在滑,就用同一个核再强迫几次。那套重复写在网格并不认识那面墙时的 delta 核。

闭式表达:

\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$ 是标记点位置。

直接强迫 —— 两个箭头的往返

现在到二维。沿圆柱表面铺标记点,每个标记点都要在流体上强制无滑移。每个时间步循环四步。

插值(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,带就变得均匀了。

一百行 IBM —— 圆柱一步#

import numpy as np
 
def peskin4(r):
    # 4 点正则化 delta (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):
    # Lagrangian 力 F_l -> Eulerian 格点力 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 一步
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}")

一步之后圆柱内部速度就接近零了。真实模拟还要加上压力 Poisson 和黏性项,但 IBM 的核心 —— 插值、计算、扩散三拍 —— 就在这三十几行里。

边界变厚的代价

Peskin 的原始 IBM 有个代价。 $\delta_h$ 一摊开,边界就变成宽度为 $h$ 的弥散界面。标记点位置准,但流体看到的壁稍微模糊。低 Re 下没问题,边界层厚到能把这模糊埋掉。Re 一高,边界层薄到不足一格,模糊壁上就解不准 BL。

压力也会振荡。标记点跨过格线时,强迫力变化是光滑的,但压力 Poisson 会把微小波纹放大。刚体情况用 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。零阶矩在这五个节点上闭合。

代码里的数是 nL = 24, R = 8, h = 1。周长 $2\pi R = 16\pi \approx 50.27$ 。间距 $2\pi R/24 \approx 2.09$ ,超过格子的两倍。经验规则 $\Delta s_l < h$ 要等标记超过 51 个才成立。64 个时 $50.27/64 \approx 0.79$ ,这时才比 h 窄。24 个标记只演示插值和铺力的一步,不是一堵闭合的墙。

列表里的 24 个标记已经越出间距规则#

间距大约 $2.09h$ 的圆周,格子之间会留下空弧。64 个标记让环带变均匀,和列表里的 24 个不是同一面墙。只看一步打印就当成 no-slip 已经合上,就会漏掉稀疏标记留下的孔。

瓣膜在成为一串点之前,先是一层薄膜。切得越薄、剪切越会假性变硬的那一侧,写在壳单元的剪切锁定。流体侧的权重和凑成 1 也没用,膜自己锁住时,瓣膜角度会因另一件事而错。