Una válvula cardíaca sin malla ajustada — el método de frontera inmersa de Peskin y la delta discreta

En 1972 Charles Peskin se quedó atascado con la malla de una válvula cardíaca en Columbia. La válvula se dobla, se hincha y se cierra en cada latido. Volver a mallar ese movimiento en cada paso no era viable. La malla se deja quieta. El cuerpo pasa a ser puntos que filtran fuerza al fluido.

Al arrastrar el marcador entre nodos, los pesos de cuatro puntos se reordenan y la suma se queda cerca de 1.

Al ensanchar el núcleo desde un sombrero de dos puntos hasta el de seis de Peskin, la fuerza se reparte más. Conviene mirar la suma cuando el marcador cruza una línea de malla: no debería saltar fuera de 1.

1972, Peskin frente a una válvula cardíaca#

Lo más espinoso en la simulación cardíaca era la malla. Cuando una válvula se mueve, un esquema ALE (Arbitrary Lagrangian–Eulerian) obliga a la malla a seguirla, pero el movimiento es enorme y veloz. La idea de Peskin fue desacoplar por completo los dos sistemas de coordenadas. El fluido se resuelve sobre una malla cartesiana fija (Euleriana) y el cuerpo se representa como una nube móvil de puntos (Lagrangiana). El puente entre ambos marcos es la delta de Dirac $\delta$.

En una línea:

Aquí $\mathbf{x}$ es un nodo de fluido, $\mathbf{X}(s,t)$ es un marcador Lagrangiano y $\mathbf{F}$ es la fuerza que el marcador entrega al fluido. La forma y el movimiento de la válvula caben dentro de una única función $\mathbf{X}(s,t)$.

No volver a mallar#

El problema es la computadora. Una computadora no puede manejar una delta de Dirac directamente. Una función que sube a infinito en una región más estrecha que el paso $h$ no se discretiza. Peskin introdujo la delta discreta $\delta_h$ — la delta de Dirac reemplazada por una colina suave de ancho $h$.

Debe cumplir dos condiciones.

• La suma de los pesos sobre todos los nodos vale exactamente 1 (momento de orden cero)
• La suma ponderada de coordenadas iguala la coordenada del marcador (momento de orden uno)

Suma igual a 1 significa que la fuerza aplicada se conserva. El primer momento garantiza que el par no se desvía. Solo esas dos condiciones determinan de forma única el núcleo de cuatro puntos.

La delta discreta — un punto que se reparte en la malla#

Desde un sombrero de dos puntos hasta el núcleo de seis puntos de Peskin, cuanto más ancho el soporte, más suavemente se reparten los pesos. Cuanto más suave el núcleo, menos brincos da la fuerza cuando el marcador cruza una línea de malla. Por eso IBM produce menos oscilaciones temporales que una interpolación abrupta. El núcleo de cuatro puntos de Peskin es el caballo de batalla porque cumple, además de las dos condiciones anteriores, una tercera: la suma de pesos al cuadrado no depende de la posición del marcador. Si una pasada todavía resbala, el mismo núcleo se fuerza otra vez. Esa repetición está en el núcleo delta cuando la malla no conoce la pared.

En forma cerrada:

con $r = (x_i - x_L)/h$, donde $x_i$ es un nodo y $x_L$ es el marcador.

Forzado directo — un viaje de ida y vuelta entre dos flechas#

Pasamos a 2D. Se colocan marcadores a lo largo de un cilindro circular y se exige que cada uno imponga la condición de no deslizamiento. En cada paso temporal se ejecutan cuatro acciones.

1. Interpolación (E → L): llevar la velocidad de la malla al marcador.

2. Cálculo de la fuerza: dividir la diferencia entre la velocidad deseada $\mathbf{U}_l^{des}$ y la interpolada por una escala temporal.

3. Reparto (L → E): empujar esa fuerza de vuelta a la malla.

4. Actualización del fluido: $\mathbf{u}^{n+1} = \mathbf{u}^n + \Delta t\,(\dots + \mathbf{f}/\rho)$.

La interpolación y el reparto usan el mismo núcleo $\delta_h$. Esa simetría devuelve la fuerza del marcador como momento del fluido.

Con muy pocos marcadores quedan huecos en el borde. La regla práctica es $\Delta s_l < h$ — el espaciamiento entre marcadores debe ser menor que la malla. Subir a 64 marcadores deja la banda uniforme.

Cien líneas de IBM — un paso alrededor del cilindro#

import numpy as np
 
def peskin4(r):
    # delta regularizada de 4 puntos (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):
    # fuerza Lagrangiana F_l -> fuerza de malla 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):
    # velocidad de malla u_ij -> velocidad de marcador 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
 
# un paso de IBM con forzado directo
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                      # flujo de entrada uniforme
U_des = np.zeros((nL, 2))             # cilindro fijo = no deslizamiento
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}")

Después de un solo paso, la velocidad interior del cilindro ya colapsa hacia cero. Una simulación real necesita la Poisson de presión y los términos viscosos, pero el corazón de IBM — interpolar, calcular, repartir — cabe en estas treinta líneas.

El precio de una pared difusa#

El IBM original de Peskin paga un precio. Al extender $\delta_h$ una vez, la frontera se vuelve una interfaz difusa de ancho $h$. La posición del marcador es exacta, pero la pared que ve el fluido queda algo borrosa. A bajo Re no hay problema, la capa límite es lo bastante gruesa como para absorber la difuminación. A alto Re la capa límite cae por debajo de una celda y sobre esa pared difusa la BL ya no se resuelve.

La presión también oscila. Cuando los marcadores cruzan líneas de malla, el forzado cambia con suavidad, pero la Poisson de presión amplifica esas pequeñas ondulaciones. Para cuerpos rígidos, un IBM implícito controla las oscilaciones, a cambio de un sistema lineal extra por paso.

Por qué aparecieron las variantes sharp-interface#

Estos límites empujaron al campo, desde finales de los 90, hacia las variantes de IBM con forzado discreto. El ghost-cell IBM de Mittal e Iaccarino, el cut-cell de Tseng y Ferziger, los esquemas sharp-interface de Tucker. Todos comparten una idea — dejar $\delta_h$ y resolver un problema local de interpolación en las celdas cercanas a la frontera para imponer la condición exacta.

Las membranas flexibles, como las válvulas cardíacas, siguen siendo terreno natural de la familia Peskin; los cuerpos rígidos en flujo externo a alto Re tienden hacia ghost-cell. No compiten, se reparten el trabajo.

Sobre un nodo los pesos son 1/2 y 1/4, y 24 marcadores ya pasan de 2h#

Si el marcador cae justo en un nodo, las distancias que quedan son $r = 0$, $\pm 1$, $\pm 2$. En $r = 0$ la primera pieza vale $(3 + \sqrt{1})/8 = 1/2$. En $|r| = 1$ las dos piezas se encuentran, ambas en $(3 - 2 + \sqrt{1})/8 = 1/4$. En $|r| = 2$ el peso es 0. La fila es $0$, $1/4$, $1/2$, $1/4$, $0$, y la suma es 1. El momento de orden cero se cierra en esos cinco nodos.

El listado fija nL = 24, R = 8, h = 1. La circunferencia es $2\pi R = 16\pi \approx 50.27$. El espaciado $2\pi R/24 \approx 2.09$ supera el doble de la malla. La regla $\Delta s_l < h$ pide más de 51 marcadores. Con 64, $50.27/64 \approx 0.79$, y solo entonces el espaciado queda bajo h. Veinticuatro marcadores muestran un paso de interpolar y repartir. No son una pared cerrada.

Los 24 marcadores del listado ya rompen la regla de espaciado#

Una circunferencia con separación de unos $2.09h$ deja arcos vacíos entre nodos. Que 64 marcadores uniformicen la banda no describe la misma pared que los 24 del listado. Leer la impresión de un solo paso como un no-slip cerrado pasa por alto los huecos que dejan los marcadores gruesos.

Una válvula es una membrana delgada antes de ser una cadena de puntos. El lado en que un cascarón delgado se endurece en cortante por la razón equivocada está en el bloqueo a cortante de los elementos de cáscara. Igualar a 1 la suma de pesos del fluido no corrige el ángulo de la válvula si la membrana se bloquea por su cuenta.