流固耦合求解器

耦合策略

PhysEngine采用分区型耦合策略，SPH和FEM求解器独立运行，通过专用界面交换数据。

graph LR
    subgraph SPH求解器
        A1[粒子追踪] --> A2[压力计算]
        A2 --> A3[密度更新]
    end

    subgraph 界面数据交换
        B1[压力映射] <--> B2[位移/速度映射]
    end

    subgraph FEM求解器
        C1[动力学求解] --> C2[应力计算]
        C2 --> C3[变形更新]
    end

    A3 --> B1
    B2 --> C1

耦合界面

界面条件

  1. 位移协调：$\mathbf{d}{FSI} = \mathbf{d}$
  2. 速度协调：$\mathbf{v}{FSI} = \mathbf{v}$
  3. 力平衡：$\mathbf{f}{fluid} = \mathbf{f}$

压力-速度映射

从SPH粒子压力到FEM节点力的映射：

$$f_i = \sum_b A_{ib} P_b V_b$$

其中 $A_{ib}$ 是从粒子 $b$ 到节点 $i$ 的映射权重。

耦合算法

主从耦合流程

┌──────────────────────────────────────────────────────────────┐
│                      时间步 n → n+1                           │
├──────────────────────────────────────────────────────────────┤
│  1. 界面搜索：                                                 │
│     - 识别SPH粒子与FEM边界的最近点                            │
│     - 建立映射关系                                            │
├──────────────────────────────────────────────────────────────┤
│  2. SPH求解（子步）：                                         │
│     - 从FEM获取边界位移/速度（Dirichlet边界）                │
│     - 执行一个SPH时间步                                       │
│     - 计算界面压力                                            │
├──────────────────────────────────────────────────────────────┤
│  3. 界面映射：                                                │
│     - 将SPH压力映射到FEM边界（Neumann边界）                  │
│     - 将FEM位移/速度映射到SPH边界                              │
├──────────────────────────────────────────────────────────────┤
│  4. FEM求解（子步）：                                         │
│     - 应用流体压力载荷                                        │
│     - 执行FEM时间步                                           │
│     - 输出结构响应                                            │
└──────────────────────────────────────────────────────────────┘

数值稳定性

耦合系统的时间步长由两子系统共同决定：

$$\Delta t = \min\left(\Delta t_{SPH}, \Delta t_{FEM}\right)$$

时间步长估计

子系统 约束条件
SPH CFL条件、声波传播
FEM Newmark稳定性要求、结构特征周期

边界处理

无滑移边界（SPH侧）

$$\mathbf{v}{particle} = \mathbf{v}$$

压力边界（FEM侧）

$$\mathbf{t} = -P \mathbf{n}$$

其中 $\mathbf{n}$ 为界面法向单位向量。