A Heart Valve Without a Body-Fitted Mesh — Peskin's Immersed Boundary Method and the Discrete Delta

In 1972 Charles Peskin got stuck on the mesh around a heart valve at Columbia. The valve flexes, bulges, and shuts on every beat. Remeshing that motion every step was not going to happen. Leave the grid alone. Let the body be points that leak force into the fluid.

Drag the marker between grid nodes and the four-point weights rearrange. Their sum stays near 1.

Widening the kernel from a 2-point hat to the 6-point Peskin spreads the force. Watch the sum when the marker crosses a grid line: it should not jump off 1.

1972, Peskin in front of a heart valve#

The worst part of heart simulation was the mesh. When a valve moves, an ALE (Arbitrary Lagrangian–Eulerian) approach forces the grid to follow, but the motion is too large and too fast. Peskin's idea was to decouple the two frames completely. Solve the fluid on a fixed Cartesian grid (Eulerian) and represent the body as a moving cloud of points (Lagrangian). The bridge between the frames is the Dirac delta function $\delta$.

One line says it:

Here $\mathbf{x}$ is a fluid grid point, $\mathbf{X}(s,t)$ is a Lagrangian marker, and $\mathbf{F}$ is the force the marker hands to the fluid. The valve's whole shape and motion live inside a single function $\mathbf{X}(s,t)$.

Don't remesh#

The trouble is the computer. A computer cannot handle a Dirac delta directly. A function that spikes to infinity inside a region narrower than $h$ refuses to be discretized. Peskin introduced the discrete delta $\delta_h$ — the Dirac delta replaced by a smooth bump of width $h$.

It has to satisfy two conditions.

• The weights sum to 1 over all grid nodes (zeroth moment)
• The weighted sum of coordinates equals the marker coordinate (first moment)

Sum equals 1 means the applied force is conserved. The first moment guarantees torque doesn't drift. Just those two pin down the 4-point kernel uniquely.

The discrete delta — spreading one point onto the grid#

From a 2-point hat to a 6-point Peskin kernel, the wider the support, the more smoothly the weights are smeared. A smoother kernel means the force does not jump when the marker slides across a grid line. That is why IBM produces fewer temporal oscillations than a sharp interpolation would. The 4-point Peskin kernel is the workhorse — it satisfies a third condition on top of the two above: the sum of squared weights is independent of marker position. If one pass still slips, the same kernel is forced again. That repeat lives in the delta kernel on a grid that does not know the wall.

In closed form:

with $r = (x_i - x_L)/h$, where $x_i$ is a grid node and $x_L$ is the marker.

Direct forcing — a round trip between two arrows#

Now move to 2D. Place markers along a circular cylinder and demand that each marker enforce no-slip. Every timestep cycles through four steps.

1. Interpolate (E → L): pull grid velocity $\mathbf{u}$ to the marker location.

2. Compute force: divide the velocity mismatch by a timescale.

3. Spread (L → E): push that force back onto the grid.

4. Update fluid: $\mathbf{u}^{n+1} = \mathbf{u}^n + \Delta t\,(\dots + \mathbf{f}/\rho)$.

Interpolation and spreading use the same kernel $\delta_h$. That symmetry returns the marker force as fluid momentum.

If you place too few markers, gaps open along the rim. The rule of thumb is $\Delta s_l < h$ — the marker spacing must be tighter than the grid. Crank it up to 64 and the band becomes uniform.

A hundred lines of IBM — one step around a cylinder#

import numpy as np
 
def peskin4(r):
    # 4-point regularized 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 force F_l -> Eulerian grid force 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):
    # grid velocity u_ij -> marker velocity 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
 
# one step of 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                      # uniform inflow
U_des = np.zeros((nL, 2))             # stationary cylinder = 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}")

After a single step, the interior velocity already collapses toward zero. A real simulation needs a pressure Poisson solve and viscous terms, but the core of IBM — interpolate, compute, spread — sits in those thirty lines.

The price of a blurred wall#

Peskin's original IBM has a cost. Smear $\delta_h$ once and the boundary becomes a diffuse interface of width $h$. The marker position is exact, but the wall the fluid sees is slightly fuzzy. At low Re, no harm done — the boundary layer is thick enough to swallow the smear. At high Re the boundary layer shrinks below one grid cell, and a fuzzy wall ruins boundary-layer resolution.

Pressure oscillates too. As markers cross grid lines the forcing changes smoothly, yet the pressure Poisson amplifies even tiny ripples. Implicit IBM tames the oscillations for rigid bodies, at the cost of one extra linear solve per timestep.

Why sharp-interface variants appeared#

These limitations pushed the field, from the late 1990s onward, toward discrete-forcing IBM variants. Mittal and Iaccarino's ghost-cell IBM, Tseng and Ferziger's cut-cell, Tucker's sharp-interface schemes. They share one idea — drop $\delta_h$ and instead solve a local interpolation problem in the boundary cells so that the wall condition is imposed exactly.

Flexible membranes like heart valves stay in the Peskin family; rigid bodies in high-Re external flow lean ghost-cell. They are not rivals, they share the workload.

On a node the weights are 1/2 and 1/4, and 24 markers already exceed 2h#

Sit the marker on a node and the only offsets left are $r = 0$, $\pm 1$, $\pm 2$. At $r = 0$ the first piece is $(3 + \sqrt{1})/8 = 1/2$. At $|r| = 1$ the two pieces meet, both $(3 - 2 + \sqrt{1})/8 = 1/4$. At $|r| = 2$ the weight is 0. The row is $0$, $1/4$, $1/2$, $1/4$, $0$, and the sum is 1. The zeroth moment closes on those five nodes.

The listing sets nL = 24, R = 8, h = 1. The circumference is $2\pi R = 16\pi \approx 50.27$. The spacing $2\pi R/24 \approx 2.09$ is more than twice the grid. The rule $\Delta s_l < h$ needs more than 51 markers. At 64, $50.27/64 \approx 0.79$, and only then is the spacing under h. Twenty-four markers show one interpolate-and-spread step. They are not a closed wall.

The listing's 24 markers already break the spacing rule#

A ring spaced at about $2.09h$ leaves empty arcs between nodes. The note that 64 markers even out the band is not describing the same wall as the listing's 24. Read the one-step print as a closed no-slip and the holes from the coarse markers get missed.

A valve is a thin membrane before it is a chain of points. Thin shells that go stiff in shear for the wrong reason are the subject of shear locking in shell elements. Matching the fluid weights to 1 does not help if the membrane itself locks and the valve angle goes wrong for that other reason.