FLIP 流体模拟算法深度解析：从原理到实现

摘要：FLIP（Fluid-Implicit Particle）算法是计算机图形学中流体模拟的经典方法，结合了拉格朗日粒子法和欧拉网格法的优势。本文将系统阐述 FLIP 算法的数学原理、核心步骤、代码实现，以及与传统方法的对比分析。

---

目录

1. 引言：流体模拟的两条路线
2. FLIP 的核心思想
3. 数学基础
4. 算法流程详解
5. 关键步骤深入
6. 完整代码实现
7. PIC、FLIP 与 APIC 的对比
8. 性能优化与工程实践
9. 前沿发展与总结

---

1. 引言：流体模拟的两条路线

在计算机图形学中，模拟流体（水、烟雾、火焰等）主要有两种视角：

1.1 拉格朗日视角（Lagrangian）—— 粒子法

以 SPH（Smoothed Particle Hydrodynamics） 为代表，将流体离散为大量粒子，每个粒子携带质量、速度、位置等物理量，随流体运动。

优点：质量守恒天然保证，细节保留好
缺点：需要大量粒子，邻居搜索开销大，不可压缩性难以严格保证

1.2 欧拉视角（Eulerian）—— 网格法

将空间划分为固定网格，在每个网格单元上求解 Navier-Stokes 方程。

优点：压力求解精确，不可压缩性好
缺点：数值耗散严重，细节容易被抹平

1.3 FLIP 的诞生

FLIP（Fluid-Implicit Particle） 由 Brackbill 和 Ruppel 于 1986 年首次提出，后由 Zhu 和 Bridson（2005）引入图形学领域。它巧妙地融合了两种视角：

• 粒子 负责携带速度和质量（拉格朗日优势）
• 网格 负责求解压力和确保不可压缩性（欧拉优势）

这种"取其精华"的设计，使 FLIP 成为现代流体求解器（如 Houdini 的 Flip Solver、Blender 的 Mantaflow）的核心算法。

---

2. FLIP 的核心思想

2.1 直观理解

想象你有一杯水：

• 纯网格法：你画了一个固定网格，水只能在这个网格上流动，每次移动都在网格间插值，就像把水倒进一格一格的冰格里——水会慢慢失去细节。
• 纯粒子法：你用几百万个小球代表水滴，它们自由运动，但你需要复杂的计算才能让它们保持"水的体积不变"。
• FLIP：你让粒子自由运动，保留它们的细节；但每隔一段时间，你把粒子的速度"告诉"网格，让网格帮你算一下压力，确保水不会压缩或膨胀，然后把修正后的速度传回粒子。

2.2 关键洞察

FLIP 的核心洞察在于：只把"变化量"从网格传输回粒子，而不是整个速度场。

这是它与前身 PIC（Particle-in-Cell） 的根本区别：

PIC:  v_particle_new = v_grid        （完全覆盖，数值耗散）
FLIP: v_particle_new = v_particle + (v_grid_new - v_grid_old)   （只传增量，保留细节）

2.3 工作流程概览

┌─────────────────────────────────────────────────────┐
│                   FLIP 算法主循环                     │
├─────────────────────────────────────────────────────┤
│                                                       │
│  1. 粒子 → 网格 (P2G): 将粒子速度/质量转移到网格      │
│                         ↓                             │
│  2. 网格上求解: 计算压力，投影速度使其无散度          │
│                         ↓                             │
│  3. 网格 → 粒子 (G2P): 将网格的速度增量传回粒子       │
│                         ↓                             │
│  4. 粒子平流: 用更新后的速度移动粒子                  │
│                         ↓                             │
│  5. 回到步骤 1（下一帧）                              │
│                                                       │
└─────────────────────────────────────────────────────┘

---

3. 数学基础

3.1 Navier-Stokes 方程

不可压缩流体的运动由 Navier-Stokes 方程描述：

$$\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot \nabla )\mathbf{u}=-\frac{1}{\rho }\nabla p+\nu {\nabla }^2\mathbf{u}+\mathbf{f}$$

$$\nabla \cdot \mathbf{u}=0$$

其中：

• $\mathbf{u}$ — 速度场
• $p$ — 压强
• $\rho$ — 密度
• $\nu$ — 运动粘性系数
• $\mathbf{f}$ — 外力（如重力）

3.2 算子分裂（Operator Splitting）

我们不直接求解完整的 Navier-Stokes 方程，而是将其拆分为几个独立步骤依次处理：

步骤 1：外力项

$$\frac{\partial \mathbf{u}}{\partial t}=\mathbf{f}$$

离散形式：${\mathbf{u}}^{\ast }={\mathbf{u}}^n+\Delta t\cdot \mathbf{f}$

步骤 2：对流项

$$\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot \nabla )\mathbf{u}=0$$

这正是粒子做的事情——以自身速度运动，天然处理了对流。

步骤 3：压力投影（核心步骤）

$$\frac{\partial \mathbf{u}}{\partial t}=-\frac{1}{\rho }\nabla p\text{满足}\nabla \cdot {\mathbf{u}}^{n+1}=0$$

这导出著名的 压力泊松方程：

$$\nabla \cdot \left(\frac{1}{\rho }\nabla p\right)=\frac{1}{\Delta t}\nabla \cdot {\mathbf{u}}^{\ast }$$

3.3 交错网格（MAC Grid）

FLIP 通常使用 MAC（Marker-and-Cell）交错网格，不同物理量存储在不同位置：

      v(i,j+1/2)
         ↑
         |
u(i-1/2,j) —— P(i,j) —— u(i+1/2,j)
         |
         ↓
      v(i,j-1/2)

图：MAC 网格布局（2D）
- 压强 P 存储在单元中心
- 水平速度 u 存储在左右面中心
- 垂直速度 v 存储在上下边中心

这种布局的优势：

• 自然地避免奇偶失联（checkerboard instability）
• 离散散度算子 $\nabla \cdot$ 和梯度算子 $\nabla$ 互为伴随
• 压力泊松方程的系数矩阵对称正定

---

4. 算法流程详解

下面我们以 2D 不可压缩流体为例，逐步详解 FLIP 算法的每一帧计算。

Step 0: 初始化

创建粒子：在流体区域内均匀采样粒子
每个粒子存储：位置 x, 速度 v
创建 MAC 网格：尺寸 Nx × Ny

Step 1: 粒子到网格（P2G, Particle-to-Grid）

将粒子的质量和动量转移到 MAC 网格上。

1.1 质量/权重转移

对每个粒子 $p$，用插值核函数 $w$ 将其质量分配到周围网格点：

$$m_i=\sum _p{m}_p\cdot w({\mathbf{x}}_p-{\mathbf{x}}_i)$$

1.2 速度转移

将粒子的速度也按权重累积到网格：

$${\mathbf{u}}_i=\frac{{\sum }_p{m}_p{\mathbf{v}}_p\cdot w({\mathbf{x}}_p-{\mathbf{x}}_i)}{m_i}$$

插值核函数通常使用二次 B 样条或三次样条：

$$w(r)=\{\begin{array}{ll}\frac{3}{4}-|r{|}^2& 0\le |r|<\frac{1}{2}\\ \frac{1}{2}(\frac{3}{2}-|r|{)}^2& \frac{1}{2}\le |r|<\frac{3}{2}\\ 0& \text{otherwise}\end{array}$$

在 2D 中，$w(\mathbf{r})=w(r_x)\cdot w(r_y)$（可分离核）。

Step 2: 保存旧网格速度

在执行压力投影之前，保存当前的网格速度：

$${\mathbf{u}}_{\text{old}}={\mathbf{u}}_{\text{grid}}$$

这一步至关重要——FLIP 用它来计算速度的增量。

Step 3: 施加外力

将重力等外力加到网格速度上：

$$\mathbf{u}\leftarrow \mathbf{u}+\Delta t\cdot \mathbf{g}$$

其中 $\mathbf{g}=(0,-9.81)$ 表示向下的重力。

Step 4: 求解压力泊松方程

4.1 构建线性系统

对每个流体单元（有粒子存在的单元），离散散度算子：

$$(\nabla \cdot \mathbf{u}{)}_{i,j}=\frac{u_{i+1/2,j}-u_{i-1/2,j}}{\Delta x}+\frac{v_{i,j+1/2}-v_{i,j-1/2}}{\Delta y}$$

离散压力泊松方程（假设均匀网格 $\Delta x=\Delta y=h$）：

$$\frac{4p_{i,j}-p_{i+1,j}-p_{i-1,j}-p_{i,j+1}-p_{i,j-1}}{h^2}=-\frac{\rho }{\Delta t}(\nabla \cdot \mathbf{u}{)}_{i,j}$$

矩阵形式：

$$\mathbf{Ap}=\mathbf{b}$$

其中 $\mathbf{A}$ 是稀疏的拉普拉斯矩阵（五点模板），$\mathbf{b}$ 是负散度。

4.2 边界条件

• 固体边界：Neumann 边界 $\frac{\partial p}{\partial n}=0$（无穿透）
• 自由表面：Dirichlet 边界 $p=0$（大气压强）

4.3 求解

使用迭代法求解：

• 共轭梯度法（Conjugate Gradient）：适合对称正定矩阵
• 预处理共轭梯度法（PCG）：对 Jacobi 预处理效果很好
• 多重网格法（Multigrid）：大规模问题的最优选择

Step 5: 压力投影

利用求解出的压力修正速度，使其无散度：

$$u_{i+1/2,j}\leftarrow u_{i+1/2,j}-\frac{\Delta t}{\rho }\cdot \frac{p_{i+1,j}-p_{i,j}}{h}$$

$$v_{i,j+1/2}\leftarrow v_{i,j+1/2}-\frac{\Delta t}{\rho }\cdot \frac{p_{i,j+1}-p_{i,j}}{h}$$

此步骤后的速度满足 $\nabla \cdot \mathbf{u}=0$。

Step 6: 网格到粒子（G2P, Grid-to-Particle）

这是 FLIP 与 PIC 最关键的区别。

6.1 计算网格速度增量

$$\Delta \mathbf{u}={\mathbf{u}}_{\text{new}}^{\text{grid}}-{\mathbf{u}}_{\text{old}}^{\text{grid}}$$

6.2 将增量插值回粒子

先从网格插值得到该粒子所在位置的增量，然后更新粒子速度：

$$\Delta {\mathbf{v}}_p=\sum _i\Delta {\mathbf{u}}_i\cdot w({\mathbf{x}}_p-{\mathbf{x}}_i)$$

$${\mathbf{v}}_p^{\text{FLIP}}={\mathbf{v}}_p+\Delta {\mathbf{v}}_p$$

6.3 PIC-FLIP 混合（重要！）

纯 FLIP 可能引入噪声和不稳定性。实际中通常采用 PIC-FLIP 混合：

$${\mathbf{v}}_p^{\text{new}}=\alpha \cdot {\mathbf{v}}_p^{\text{FLIP}}+(1-\alpha )\cdot {\mathbf{v}}_p^{\text{PIC}}$$

其中：

• ${\mathbf{v}}_p^{\text{PIC}}$ = 直接从网格插值的速度（稳定但有耗散）
• ${\mathbf{v}}_p^{\text{FLIP}}$ = 粒子原速度 + 网格速度增量（保真但有噪声）
• $\alpha$ 通常取 0.95 ~ 0.99（非常接近纯 FLIP）

Step 7: 粒子平流

用更新后的速度移动粒子：

$${\mathbf{x}}_p^{\text{new}}={\mathbf{x}}_p+\Delta t\cdot {\mathbf{v}}_p^{\text{new}}$$

可使用更高阶的时间积分（如 RK2/RK3）来减少误差。

Step 8: 粒子重采样（可选）

长时间运行后，粒子分布可能不均匀。某些区域粒子过于稀疏，某些区域过于密集。此时需要重采样（reseed）粒子，以保证每个流体单元至少有 4-8 个粒子。

---

5. 关键步骤深入

5.1 插值核函数的实现

import numpy as np

def kernel_quadratic(r):
    """
    二次 B 样条核（常用选择）
    r: 距离，以网格间距为单位
    """
    r = abs(r)
    if r < 0.5:
        return 0.75 - r * r
    elif r < 1.5:
        t = 1.5 - r
        return 0.5 * t * t
    else:
        return 0.0

def kernel_cubic(r):
    """
    三次 B 样条核（更光滑但支撑域更大）
    """
    r = abs(r)
    if r < 1.0:
        return 0.5 * r * r * r - r * r + 2.0 / 3.0
    elif r < 2.0:
        t = 2.0 - r
        return t * t * t / 6.0
    else:
        return 0.0

def interpolate_2d(grid_x, grid_y, px, py, values, kernel):
    """
    2D 双线性样条插值
    grid_x, grid_y: 网格尺寸
    px, py: 粒子坐标（网格单位）
    values: 2D 网格值数组
    kernel: 核函数
    """
    result = 0.0
    # 核函数的支撑域半径
    radius = 1.5  # quadratic; 2.0 for cubic

    i_min = max(0, int(px - radius))
    i_max = min(grid_x - 1, int(px + radius + 1))
    j_min = max(0, int(py - radius))
    j_max = min(grid_y - 1, int(py + radius + 1))

    weight_sum = 0.0

    for i in range(i_min, i_max):
        wx = kernel(px - i)
        for j in range(j_min, j_max):
            wy = kernel(py - j)
            weight = wx * wy
            result += weight * values[i, j]
            weight_sum += weight

    if weight_sum > 1e-10:
        result /= weight_sum
    return result

5.2 粒子到网格（P2G）的完整实现

def particle_to_grid(particles, nx, ny, h, kernel):
    """
    将粒子速度转移到 MAC 网格
    particles: 粒子列表，每个粒子有 (x, y, vx, vy)
    nx, ny: 网格分辨率
    h: 网格间距
    """
    # MAC 网格：u 在 (i+0.5, j), v 在 (i, j+0.5)
    u_grid = np.zeros((nx + 1, ny))
    v_grid = np.zeros((nx, ny + 1))
    u_weight = np.zeros((nx + 1, ny))
    v_weight = np.zeros((nx, ny + 1))

    radius = 1.5  # quadratic kernel

    for p in particles:
        px, py = p.x / h, p.y / h

        # 转移到 u 网格（水平速度面，x 方向偏移 0.5）
        i_min = max(0, int(px - radius))
        i_max = min(nx, int(px + radius + 1))
        j_min = max(0, int(py - 0.5 - radius))
        j_max = min(ny - 1, int(py - 0.5 + radius + 1))

        for i in range(i_min, i_max + 1):
            wx = kernel(px - i)
            for j in range(j_min, j_max + 1):
                wy = kernel(py - (j + 0.5))
                w = wx * wy
                u_grid[i, j] += w * p.vx
                u_weight[i, j] += w

        # 转移到 v 网格（垂直速度面，y 方向偏移 0.5）
        i_min = max(0, int(px - 0.5 - radius))
        i_max = min(nx - 1, int(px - 0.5 + radius + 1))
        j_min = max(0, int(py - radius))
        j_max = min(ny, int(py + radius + 1))

        for i in range(i_min, i_max + 1):
            wx = kernel(px - (i + 0.5))
            for j in range(j_min, j_max + 1):
                wy = kernel(py - j)
                w = wx * wy
                v_grid[i, j] += w * p.vy
                v_weight[i, j] += w

    # 归一化
    mask_u = u_weight > 1e-10
    mask_v = v_weight > 1e-10
    u_grid[mask_u] /= u_weight[mask_u]
    v_grid[mask_v] /= v_weight[mask_v]

    return u_grid, v_grid, u_weight, v_weight

5.3 压力泊松方程的构建与求解

from scipy.sparse import csr_matrix, linalg

def build_pressure_system(u, v, nx, ny, h, dt, rho, fluid_mask):
    """
    构建 Ap = b 稀疏线性系统
    fluid_mask: 标记哪些单元是流体（bool 数组 nx x ny）
    """
    N = nx * ny
    # 建立单元索引映射
    cell_id = np.full((nx, ny), -1, dtype=int)
    idx = 0
    for i in range(nx):
        for j in range(ny):
            if fluid_mask[i, j]:
                cell_id[i, j] = idx
                idx += 1
    n_cells = idx

    # 构建稀疏矩阵
    row = []
    col = []
    data = []
    b = np.zeros(n_cells)

    for i in range(nx):
        for j in range(ny):
            cid = cell_id[i, j]
            if cid < 0:
                continue

            diag = 0.0
            rhs = 0.0

            # 计算散度
            div = (u[i + 1, j] - u[i, j]) / h + (v[i, j + 1] - v[i, j]) / h
            rhs = -(rho / dt) * div

            # 填充拉普拉斯模板
            # 中心
            diag = 0.0

            # 左邻居
            if i > 0 and cell_id[i - 1, j] >= 0:
                row.append(cid); col.append(cell_id[i - 1, j])
                data.append(-1.0 / (h * h))
                diag += 1.0 / (h * h)
            # 右邻居
            if i < nx - 1 and cell_id[i + 1, j] >= 0:
                row.append(cid); col.append(cell_id[i + 1, j])
                data.append(-1.0 / (h * h))
                diag += 1.0 / (h * h)
            # 下邻居
            if j > 0 and cell_id[i, j - 1] >= 0:
                row.append(cid); col.append(cell_id[i, j - 1])
                data.append(-1.0 / (h * h))
                diag += 1.0 / (h * h)
            # 上邻居
            if j < ny - 1 and cell_id[i, j + 1] >= 0:
                row.append(cid); col.append(cell_id[i, j + 1])
                data.append(-1.0 / (h * h))
                diag += 1.0 / (h * h)

            # 对角线
            row.append(cid); col.append(cid)
            data.append(diag)
            b[cid] = rhs

    A = csr_matrix((data, (row, col)), shape=(n_cells, n_cells))
    return A, b, cell_id

def solve_pressure(A, b):
    """
    使用共轭梯度法求解压力
    """
    # 使用 Jacobi 预处理器
    M = A.diagonal()
    M_inv = 1.0 / np.where(M > 1e-10, M, 1.0)
    M_precond = linalg.LinearOperator(
        A.shape, lambda x: M_inv * x
    )
    p, info = linalg.cg(A, b, M=M_precond, atol=1e-6)
    if info != 0:
        print(f"Warning: CG did not converge, info={info}")
    return p

5.4 网格到粒子（G2P）与 PIC-FLIP 混合

def grid_to_particle(particles, u_old, v_old, u_new, v_new,
                     u_weight, v_weight, nx, ny, h, kernel,
                     flip_ratio=0.97):
    """
    网格 → 粒子：FLIP + PIC 混合
    flip_ratio: FLIP 混合比例（越接近 1 越接近纯 FLIP）
    """
    radius = 1.5

    for p in particles:
        px, py = p.x / h, p.y / h

        # ---- PIC 速度：直接从新网格插值 ----
        u_pic = interpolate_mac_u(u_new, px, py, nx, ny, kernel, radius)
        v_pic = interpolate_mac_v(v_new, px, py, nx, ny, kernel, radius)

        # ---- FLIP 增量 ----
        u_old_interp = interpolate_mac_u(u_old, px, py, nx, ny, kernel, radius)
        v_old_interp = interpolate_mac_v(v_old, px, py, nx, ny, kernel, radius)

        delta_u = u_pic - u_old_interp
        delta_v = v_pic - v_old_interp

        v_flip_x = p.vx + delta_u
        v_flip_y = p.vy + delta_v

        # ---- PIC-FLIP 混合 ----
        p.vx = flip_ratio * v_flip_x + (1.0 - flip_ratio) * u_pic
        p.vy = flip_ratio * v_flip_y + (1.0 - flip_ratio) * v_pic


def interpolate_mac_u(u, px, py, nx, ny, kernel, radius):
    """在 MAC u 网格位置 (i, j+0.5) 插值"""
    val = 0.0
    w_sum = 0.0
    i_min = max(0, int(px - radius))
    i_max = min(nx, int(px + radius))
    j_min = max(0, int(py - 0.5 - radius))
    j_max = min(ny - 1, int(py - 0.5 + radius))

    for i in range(i_min, i_max + 1):
        wx = kernel(px - i)
        for j in range(j_min, j_max + 1):
            wy = kernel(py - (j + 0.5))
            w = wx * wy
            val += w * u[i, j]
            w_sum += w

    return val / w_sum if w_sum > 1e-10 else 0.0


def interpolate_mac_v(v, px, py, nx, ny, kernel, radius):
    """在 MAC v 网格位置 (i+0.5, j) 插值"""
    val = 0.0
    w_sum = 0.0
    i_min = max(0, int(px - 0.5 - radius))
    i_max = min(nx - 1, int(px - 0.5 + radius))
    j_min = max(0, int(py - radius))
    j_max = min(ny, int(py + radius))

    for i in range(i_min, i_max + 1):
        wx = kernel(px - (i + 0.5))
        for j in range(j_min, j_max + 1):
            wy = kernel(py - j)
            w = wx * wy
            val += w * v[i, j]
            w_sum += w

    return val / w_sum if w_sum > 1e-10 else 0.0

---

6. 完整代码实现

下面给出一个完整的 2D FLIP 流体求解器的 Python 实现。为清晰起见，代码做了简化，但保留了所有核心要素。

import numpy as np
from scipy.sparse import csr_matrix
from scipy.sparse.linalg import cg, LinearOperator
from dataclasses import dataclass
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation


# ===================== 粒子类型 =====================

@dataclass
class Particle:
    x: float    # x 坐标
    y: float    # y 坐标
    vx: float   # x 方向速度
    vy: float   # y 方向速度


# ===================== 核函数 =====================

def quadratic_kernel(r):
    """二次 B 样条核"""
    r = abs(r)
    if r < 0.5:
        return 0.75 - r * r
    elif r < 1.5:
        t = 1.5 - r
        return 0.5 * t * t
    return 0.0


# ===================== FLIP 求解器 =====================

class FLIPSolver2D:

    def __init__(self, nx, ny, h, dt=0.01, rho=1000.0, gravity=-9.81,
                 flip_ratio=0.97, particles_per_cell=4):
        """
        nx, ny: 网格分辨率
        h: 网格间距 (m)
        dt: 时间步长 (s)
        rho: 密度 (kg/m³)
        gravity: 重力加速度
        flip_ratio: FLIP 混合比例
        particles_per_cell: 每个单元的粒子种子数
        """
        self.nx = nx
        self.ny = ny
        self.h = h
        self.dt = dt
        self.rho = rho
        self.gravity = gravity
        self.flip_ratio = flip_ratio
        self.ppc = particles_per_cell

        # MAC 网格
        self.u = np.zeros((nx + 1, ny))      # x 方向速度
        self.v = np.zeros((nx, ny + 1))      # y 方向速度
        self.u_old = np.zeros((nx + 1, ny))  # 保存旧速度
        self.v_old = np.zeros((nx, ny + 1))
        self.u_weight = np.zeros((nx + 1, ny))
        self.v_weight = np.zeros((nx, ny + 1))

        self.particles = []
        self.kernel = quadratic_kernel
        self.kernel_radius = 1.5

    # ---------- 初始化粒子 ----------
    def init_particles(self, x_range, y_range):
        """在给定矩形区域内均匀采样粒子"""
        x_min, x_max = x_range
        y_min, y_max = y_range

        nx_cells = int((x_max - x_min) / self.h)
        ny_cells = int((y_max - y_min) / self.h)

        for i in range(nx_cells):
            for j in range(ny_cells):
                # 每个单元内随机撒 self.ppc 个粒子
                for _ in range(self.ppc):
                    rx = np.random.uniform(-0.4, 0.4) * self.h
                    ry = np.random.uniform(-0.4, 0.4) * self.h
                    px = x_min + (i + 0.5) * self.h + rx
                    py = y_min + (j + 0.5) * self.h + ry
                    self.particles.append(Particle(px, py, 0.0, 0.0))

    # ---------- P2G ----------
    def particle_to_grid(self):
        # 重置
        self.u.fill(0.0); self.v.fill(0.0)
        self.u_weight.fill(0.0); self.v_weight.fill(0.0)

        h = self.h
        r = self.kernel_radius

        for p in self.particles:
            px, py = p.x / h, p.y / h

            # ---- u 网格 (i, j+0.5) ----
            imin = max(0, int(px - r - 1))
            imax = min(self.nx, int(px + r + 1))
            jmin = max(0, int(py - 0.5 - r - 1))
            jmax = min(self.ny - 1, int(py - 0.5 + r + 1))

            for i in range(imin, imax + 1):
                wx = self.kernel(px - i)
                if wx == 0: continue
                for j in range(jmin, jmax + 1):
                    wy = self.kernel(py - (j + 0.5))
                    w = wx * wy
                    self.u[i, j] += w * p.vx
                    self.u_weight[i, j] += w

            # ---- v 网格 (i+0.5, j) ----
            imin = max(0, int(px - 0.5 - r - 1))
            imax = min(self.nx - 1, int(px - 0.5 + r + 1))
            jmin = max(0, int(py - r - 1))
            jmax = min(self.ny, int(py + r + 1))

            for i in range(imin, imax + 1):
                wx = self.kernel(px - (i + 0.5))
                if wx == 0: continue
                for j in range(jmin, jmax + 1):
                    wy = self.kernel(py - j)
                    w = wx * wy
                    self.v[i, j] += w * p.vy
                    self.v_weight[i, j] += w

        # 归一化
        mask = self.u_weight > 1e-10
        self.u[mask] /= self.u_weight[mask]
        mask = self.v_weight > 1e-10
        self.v[mask] /= self.v_weight[mask]

    # ---------- 标记流体单元 ----------
    def mark_fluid_cells(self):
        """根据粒子位置标记哪些网格单元包含流体"""
        fluid = np.zeros((self.nx, self.ny), dtype=bool)
        for p in self.particles:
            i = int(p.x / self.h)
            j = int(p.y / self.h)
            if 0 <= i < self.nx and 0 <= j < self.ny:
                fluid[i, j] = True
        return fluid

    # ---------- 施加重力 ----------
    def apply_gravity(self):
        # 重力只影响垂直方向
        self.v += self.dt * self.gravity

    # ---------- 施加速度边界条件 ----------
    def apply_boundary_conditions(self):
        """固体边界：法向速度为 0"""
        # 左右墙
        self.u[0, :] = 0.0
        self.u[self.nx, :] = 0.0
        # 底和顶
        self.v[:, 0] = 0.0
        self.v[:, self.ny] = 0.0

    # ---------- 求解压力 ----------
    def solve_pressure(self, fluid):
        nx, ny = self.nx, self.ny
        h = self.h

        # 编号流体单元
        cell_id = np.full((nx, ny), -1, dtype=int)
        idx = 0
        for i in range(nx):
            for j in range(ny):
                if fluid[i, j]:
                    cell_id[i, j] = idx
                    idx += 1
        N = idx
        if N == 0:
            return np.zeros(N), cell_id

        # 构建
        rows, cols, vals = [], [], []
        b = np.zeros(N)

        for i in range(nx):
            for j in range(ny):
                cid = cell_id[i, j]
                if cid < 0:
                    continue

                # 散度
                div = (self.u[i + 1, j] - self.u[i, j]) / h \
                    + (self.v[i, j + 1] - self.v[i, j]) / h
                b[cid] = -(self.rho / self.dt) * div

                diag = 0.0
                inv_h2 = 1.0 / (h * h)

                # 六个邻居方向（2D 为四个邻居）
                for ni, nj in [(i - 1, j), (i + 1, j), (i, j - 1), (i, j + 1)]:
                    if 0 <= ni < nx and 0 <= nj < ny:
                        nid = cell_id[ni, nj]
                        if nid >= 0:
                            rows.append(cid); cols.append(nid)
                            vals.append(-inv_h2)
                            diag += inv_h2

                rows.append(cid); cols.append(cid)
                vals.append(diag)

        A = csr_matrix((vals, (rows, cols)), shape=(N, N))
        # PCG
        M_diag = A.diagonal()
        M_inv = np.where(M_diag > 1e-10, 1.0 / M_diag, 0.0)
        M = LinearOperator((N, N), lambda x: M_inv * x)
        p, info = cg(A, b, M=M, atol=1e-5, maxiter=500)
        return p, cell_id

    # ---------- 压力投影 ----------
    def pressure_projection(self, p, cell_id):
        h = self.h
        dt, rho = self.dt, self.rho
        scale = dt / (rho * h)

        for i in range(self.nx - 1):
            for j in range(self.ny):
                c_left = cell_id[i, j]
                c_right = cell_id[i + 1, j]
                dp = 0.0
                if c_left >= 0 and c_right >= 0:
                    dp = p[c_right] - p[c_left]
                elif c_left >= 0:
                    dp = 0.0 - p[c_left]
                elif c_right >= 0:
                    dp = p[c_right] - 0.0
                self.u[i + 1, j] -= scale * dp

        for i in range(self.nx):
            for j in range(self.ny - 1):
                c_bot = cell_id[i, j]
                c_top = cell_id[i, j + 1]
                dp = 0.0
                if c_bot >= 0 and c_top >= 0:
                    dp = p[c_top] - p[c_bot]
                elif c_bot >= 0:
                    dp = 0.0 - p[c_bot]
                elif c_top >= 0:
                    dp = p[c_top] - 0.0
                self.v[i, j + 1] -= scale * dp

    # ---------- G2P ----------
    def grid_to_particle(self):
        h = self.h
        r = self.kernel_radius

        for p in self.particles:
            px, py = p.x / h, p.y / h

            # PIC 速度
            u_pic = self._interp_u(self.u, px, py, r)
            v_pic = self._interp_v(self.v, px, py, r)

            # 旧网格速度
            u_old = self._interp_u(self.u_old, px, py, r)
            v_old = self._interp_v(self.v_old, px, py, r)

            # FLIP 更新
            v_flip_x = p.vx + (u_pic - u_old)
            v_flip_y = p.vy + (v_pic - v_old)

            # 混合
            alpha = self.flip_ratio
            p.vx = alpha * v_flip_x + (1 - alpha) * u_pic
            p.vy = alpha * v_flip_y + (1 - alpha) * v_pic

    def _interp_u(self, grid, px, py, r):
        val = 0.0; w_sum = 0.0
        imin = max(0, int(px - r - 1))
        imax = min(self.nx, int(px + r + 1))
        jmin = max(0, int(py - 0.5 - r - 1))
        jmax = min(self.ny - 1, int(py - 0.5 + r + 1))
        for i in range(imin, imax + 1):
            wx = self.kernel(px - i)
            for j in range(jmin, jmax + 1):
                w = wx * self.kernel(py - (j + 0.5))
                val += w * grid[i, j]
                w_sum += w
        return val / w_sum if w_sum > 1e-10 else 0.0

    def _interp_v(self, grid, px, py, r):
        val = 0.0; w_sum = 0.0
        imin = max(0, int(px - 0.5 - r - 1))
        imax = min(self.nx - 1, int(px - 0.5 + r + 1))
        jmin = max(0, int(py - r - 1))
        jmax = min(self.ny, int(py + r + 1))
        for i in range(imin, imax + 1):
            wx = self.kernel(px - (i + 0.5))
            for j in range(jmin, jmax + 1):
                w = wx * self.kernel(py - j)
                val += w * grid[i, j]
                w_sum += w
        return val / w_sum if w_sum > 1e-10 else 0.0

    # ---------- 粒子平流 ----------
    def advect_particles(self):
        for p in self.particles:
            p.x += self.dt * p.vx
            p.y += self.dt * p.vy

            # 边界钳制
            p.x = np.clip(p.x, self.h * 0.1, self.nx * self.h - self.h * 0.1)
            p.y = np.clip(p.y, self.h * 0.1, self.ny * self.h - self.h * 0.1)

    # ---------- 删除越界粒子 ----------
    def remove_out_of_bounds(self):
        self.particles = [p for p in self.particles
                          if 0 <= p.x < self.nx * self.h
                          and 0 <= p.y < self.ny * self.h]

    # ---------- 主循环 ----------
    def step(self):
        # 1. P2G
        self.particle_to_grid()

        # 2. 保存旧速度
        self.u_old = self.u.copy()
        self.v_old = self.v.copy()

        # 3. 外力
        self.apply_gravity()

        # 4. 边界条件
        self.apply_boundary_conditions()

        # 5. 标记流体单元 + 求解压力 + 投影
        fluid = self.mark_fluid_cells()
        p, cell_id = self.solve_pressure(fluid)
        self.pressure_projection(p, cell_id)

        # 6. 再次应用边界条件
        self.apply_boundary_conditions()

        # 7. G2P
        self.grid_to_particle()

        # 8. 平流 + 清理
        self.advect_particles()
        self.remove_out_of_bounds()


# ===================== 可视化 =====================

def visualize(flip, ax):
    """绘制粒子位置"""
    ax.clear()
    xs = [p.x for p in flip.particles]
    ys = [p.y for p in flip.particles]
    ax.scatter(xs, ys, s=2, c='dodgerblue')
    ax.set_xlim(0, flip.nx * flip.h)
    ax.set_ylim(0, flip.ny * flip.h)
    ax.set_aspect('equal')
    ax.set_title(f"Particles: {len(flip.particles)}")


if __name__ == '__main__':
    # 参数设置
    NX, NY = 40, 40       # 网格分辨率
    H = 0.025              # 网格间距 (m) — 对应 1m × 1m 域
    DT = 0.005             # 时间步长
    FLIP_RATIO = 0.97

    solver = FLIPSolver2D(NX, NY, H, dt=DT, flip_ratio=FLIP_RATIO,
                          particles_per_cell=4)

    # 初始化一个水块
    solver.init_particles((0.2, 0.6), (0.3, 0.7))

    # 运行动画
    fig, ax = plt.subplots(figsize=(6, 6))

    def update(frame):
        for _ in range(2):  # 每帧多步
            solver.step()
        visualize(solver, ax)
        return ax,

    ani = FuncAnimation(fig, update, frames=200, interval=30, blit=False)
    plt.show()

运行效果说明

运行上述代码，你将看到：

• 蓝色粒子代表流体（水）
• 水在重力作用下下落，撞击底壁
• 飞溅和水花自然形成
• 随着时间推移，粒子维持较好的体积守恒

---

7. PIC、FLIP 与 APIC 的对比

7.1 PIC（Particle-in-Cell）

速度更新：v_new = interp(u_grid_new)

• ✅ 稳定，无条件不会发散
• ❌ 数值耗散严重，动能快速损失，流体看起来像"糖浆"

7.2 FLIP

速度更新：v_new = v_old + (interp(u_grid_new) - interp(u_grid_old))

• ✅ 保留所有粒子细节，飞溅效果自然
• ✅ 几乎无数值耗散
• ❌ 可能引入噪声和不稳定性（纯 FLIP 在有剪切流的场景下）
• ❌ 需要更多的粒子来抑制噪声

7.3 APIC（Affine Particle-in-Cell）

Jiang et al. (2015) 提出的改进版本。在 FLIP 基础上，每个粒子额外存储一个 仿射速度矩阵 C，从而更精确地保持角动量。

v(x) = v_p + C_p (x - x_p)     ← 仿射速度场

特性	PIC	FLIP	APIC
数值耗散	高	极低	极低
稳定性	高	中等	高
角动量守恒	❌	近似	✅
实现复杂度	简单	简单	中等
计算开销	低	低	略高

7.4 为什么实际中用 PIC-FLIP 混合

纯 FLIP（α=1.0）虽然保真度最高，但在长时间模拟中可能出现：

1. 粒子聚类：多个粒子速度渐趋相同，聚成一团
2. 速度场噪声：高频速度扰动不断累积
3. 粒子分离：个别粒子获得不合理的速度飞出

加入少量 PIC（α=0.95~0.99）可以"平滑"掉这些不稳定因素，代价是极其微小的耗散。这是工业界的标准做法。

---

8. 性能优化与工程实践

8.1 稀疏数据结构

对于大规模模拟，大部分网格单元是空的。使用稀疏数据结构：

# 不存储整个网格，只存储活跃区域
from collections import defaultdict

active_u = {}     # (i, j) → velocity
active_v = {}
active_fluid = set()

8.2 GPU 并行化

P2G 和 G2P 步骤天然适合 GPU：

• 每个粒子的插值与累积是独立的
• 使用 CUDA / Taichi / OpenCL 可实现百倍加速

参考实现框架：

• NVIDIA FleX：基于位置的粒子动力学
• Taichi：嵌入 Python 的高性能并行计算语言，特别适合写 FLIP
• mantaflow：Blender 使用的开源流体求解器

8.3 自适应时间步长

CFL 条件限制了 $\Delta t$：

$$\Delta t<\frac{\Delta x}{\max |\mathbf{u}|}$$

即一个时间步内，流体移动的距离不应超过一个网格单元。实践中可在每个时间步检查 CFL 数，动态调整 $\Delta t$。

8.4 水平集辅助表面重建

粒子本身不直接定义表面。常用的表面重建方法：

• Marching Cubes：用水平集函数从粒子重建等值面
• Screen Space Fluids：在屏幕空间做平滑，适合实时应用
• Anisotropic Kernels：各向异性核函数更好地重建薄片和飞溅

8.5 粒子重采样（Reseeding）

def reseed_particles(flip, min_ppc=2, max_ppc=16):
    """确保每个流体单元有合理的粒子数"""
    # 统计每个单元的粒子数
    cell_counts = {}
    for p in flip.particles:
        i, j = int(p.x / flip.h), int(p.y / flip.h)
        cell_counts[(i, j)] = cell_counts.get((i, j), 0) + 1

    # 粒子过少 → 分裂（或新增）
    # 粒子过多 → 随机合并
    for (i, j), count in cell_counts.items():
        if count < min_ppc:
            # 在该单元中心添加粒子
            for _ in range(min_ppc - count):
                cx = (i + 0.5) * flip.h
                cy = (j + 0.5) * flip.h
                flip.particles.append(
                    Particle(cx, cy, 0, 0)
                )

---

9. 前沿发展与总结

9.1 相关研究方向

• Material Point Method（MPM）：FLIP 推广到固体力学，用于模拟雪、沙子、弹性体。迪士尼的《冰雪奇缘》即使用 MPM 模拟雪。
• Narrow Band FLIP：只在自由表面附近使用 FLIP 增量，内部使用 PIC，进一步减少噪声（Ferstl et al., 2016）。
• Learned FLIP：使用神经网络替代压力求解器或插值核函数（Ummenhofer et al., 2020）。
• DualSPHysics + FLIP hybrid：结合 SPH 和 FLIP 用于大规模自由表面流。

9.2 工业应用

应用领域	代表工具
视觉特效（VFX）	Houdini Flip Solver, Maya Bifröst
实时游戏	NVIDIA FleX, Unreal Engine Niagara
开源	Blender (mantaflow), OpenFOAM
科研	Gerris, Basilisk

9.3 总结

FLIP 算法之所以成为流体模拟的主流方法，是因为它巧妙地结合了两种方法的优势：

“粒子负责自由与细节，网格负责约束与物理。”

核心就三点：

1. P2G — 把粒子速度传到网格
2. 压力投影 — 在网格上求解泊松方程，确保不可压缩
3. G2P — 只把速度的变化量传回粒子（这是 FLIP 区别于 PIC 的灵魂所在）

掌握 FLIP，你就掌握了理解 Houdini、Blender 等 DCC 工具中流体求解器的钥匙。希望本文能帮助你建立起对 FLIP 算法的系统认知。

---

参考文献

1. Brackbill, J. U., & Ruppel, H. M. (1986). FLIP: A method for adaptively zoned, particle-in-cell calculations of fluid flows in two dimensions. Journal of Computational Physics, 65(2), 314-343.
2. Zhu, Y., & Bridson, R. (2005). Animating sand as a fluid. ACM Transactions on Graphics (TOG), 24(3), 965-972.
3. Bridson, R. (2015). Fluid Simulation for Computer Graphics (2nd ed.). CRC Press.
4. Jiang, C., Schroeder, C., Selle, A., Teran, J., & Stomakhin, A. (2015). The affine particle-in-cell method. ACM Transactions on Graphics (TOG), 34(4), 1-10.
5. Ferstl, F., Ando, R., Wojtan, C., Westermann, R., & Thuerey, N. (2016). Narrow band FLIP for liquid simulations. Computer Graphics Forum, 35(2), 225-232.
6. Ummenhofer, B., Prantl, L., Thuerey, N., & Koltun, V. (2020). Lagrangian fluid simulation with continuous convolutions. International Conference on Learning Representations (ICLR).

---

作者注：本文所有代码均为教学目的简化实现。生产级代码建议参考 mantaflow 或 Taichi 示例。如有疑问或讨论，欢迎在评论区交流！