游戏实时风力模拟分析

发表于 2025-08-24 分类于工业阅读次数：本文字数： 1.3k 阅读时长 ≈ 5 分钟

一、论文核心贡献

创新点
- 提出首个实时Navier-Stokes方程求解器，专为游戏引擎优化（非物理精确解）。
- 核心突破：无条件稳定性（任意时间步长不崩溃），牺牲物理精度换取视觉合理性与速度。
- 完整C代码实现仅需100+行，支持2D/3D网格实时运算（标准PC硬件）。
关键技术
- 分层求解思想：先解密度场（线性），再扩展至速度场（非线性）。
- 半拉格朗日法（advect()函数）：逆向追踪插值避免数值爆炸。
- 质量守恒投影（project()函数）：霍奇分解消除速度场散度，生成漩涡效果。

二、算法流程详解

1. 密度场求解（基础层）

void dens_step(...) {
  add_source(...);   // 注入源（如烟雾）
  SWAP; diffuse(...); // 扩散（稳定高斯-赛德尔迭代） 
  SWAP; advect(...);  // 平流（半拉格朗日法）
}

扩散：
使用20次高斯-赛德尔松弛求解线性系统，避免显式差分不稳定性：
```
x[IX(i,j)] = (x0[...] + a*(相邻四项和)) / (1+4*a); 
```

平流：
逆向追踪网格点来源位置，双线性插值获取密度：

x = i - dt0*u[IX(i,j)];  // 逆向追踪
d[IX(i,j)] = 相邻四点插值; // 双线性插值

2. 速度场求解（进阶层）

void vel_step(...) {
  add_source(...);     // 注入外力
  diffuse(...);        // 粘性扩散
  project(...);        // 质量守恒投影
  advect(...);         // 自平流（速度场跟随自身运动）
  project(...);        // 二次投影增强稳定性
}

关键创新：质量守恒投影
通过霍奇分解将速度场拆解为无散度场（漩涡） + 梯度场（非物理振荡）：

// 步骤1：计算散度场
div[IX(i,j)] = -0.5h*(u右-u左 + v上-v下); 

// 步骤2：解泊松方程求高度场p
p[IX(i,j)] = (div[...] + 相邻四项和)/4; 

// 步骤3：从速度场减去梯度场
u[IX(i,j)] -= 0.5*(p右-p左)/h;

3. 边界处理（set_bnd()）

固体边界：速度场法向分量归零（墙面无穿透）

x[IX(0,i)] = (b==1) ? -x[IX(1,i)] : x[IX(1,i)]; // 水平墙反射

角点特殊处理：避免数值歧义

x[IX(0,0)] = 0.5*(x[IX(1,0)] + x[IX(0,1)]);

三、NS方程原始形式与算法项对应

论文图1给出NS方程的核心形式（简化向量表示）：

右边三项分别对应：

→ **advect()**（自平流项）
- 物理意义：流体自身惯性导致的动量输运（非线性项）
- 算法实现：通过半拉格朗日法逆向追踪速度场（详见§4.3）：
```
x = i - dt0*u[IX(i,j)];  // 逆向追踪质点位置
u_new = 双线性插值(u_old) // 用前一时刻速度场插值
```
→ **diffuse()**（粘性扩散项）
- 物理意义：流体粘性导致的动量扩散（耗散项）
- 算法实现：隐式高斯-赛德尔迭代求解扩散方程（§4.2）：
```
// 解方程：u - dt*ν∇²u = u_old
u[IX(i,j)] = (u_old + a*(相邻速度之和)) / (1+4*a);
```
→ **add_source()**（外力项）
- 物理意义：外部力场（如重力、风扇）的输入
- 算法实现：直接叠加外力到速度场（§4.1）：
```
u[i] += dt * f_ext[i]; 
```

四、压力项的隐藏角色与`project()`

NS方程在不可压缩流体中需补充质量守恒约束：

压力项由此衍生：

物理本质：压力是维持流体不可压缩的约束力（非独立外力）
算法对应 → **project()**（质量守恒投影）
- 通过霍奇分解（图7）将速度场拆解：
  $无散$
- 核心操作（§4.4）：
  1. 计算速度场散度 → 质量亏损量
```
div = -0.5h*(u_right - u_left + v_top - v_bottom); 
```
  2. 解泊松方程→ 压力场
```
p[IX(i,j)] = (div + 相邻p之和)/4; // 高斯-赛德尔迭代
```
  3. 从速度场中减去压力梯度 → 无散速度场
```
u -= 0.5*(p_right - p_left)/h; // 消除质量亏损
```

五、算法流程与NS方程的完整映射

物理过程	算法函数	执行顺序
外力输入	`add_source()`	第一步
粘性扩散	`diffuse()`	第二步
质量守恒约束（压力）	`project()`	第三步
自平流	`advect()`	第四步

✅ 关键说明：
压力项无独立函数：其效果通过project()中间接实现（泊松方程求解）
project()双重调用：首次在扩散后修正粘性导致的质量亏损，二次在平流后修正惯性导致的质量亏损（§4.4代码）
粘滞力与扩散项：严格对应diffuse()的物理意义

六、边界条件与数值稳定性

固体边界（set_bnd()）：强制法向速度为零（§4.4）：

u_wall = (边界类型==1) ? -u_adjacent : u_adjacent; // 反射法向分量

稳定性保障：
- 扩散项用隐式格式（无条件稳定）
- 平流项用半拉格朗日法（逆向追踪避免数值振荡）

此设计确保算法可处理任意时间步长（论文§2强调），突破传统CFL条件限制。

七、参考

1.Navier-Stokes方程式介紹與推導
2.描述流場性質的2種方法 – Eulerian v.s. Lagrange，固定座標系下流動質點的物理性質變化表示