基于 WebGPU 的欧拉网格流体模拟 (Eulerian Grid Fluid)
基于 WebGPU 的欧拉网格流体模拟(AI编写)
在计算机图形学中,流体模拟主要分为两大流派:拉格朗日视角 (Lagrangian) 和 欧拉视角 (Eulerian)。
- 拉格朗日视角:将流体看作无数个离散的粒子 (Particle)。每个粒子携带位置、速度等信息。SPH (Smoothed Particle Hydrodynamics) 就是典型的代表。
- 欧拉视角:将空间划分为固定的网格 (Grid)。我们关注的是每个网格点上的物理量(如密度、速度、压力)随时间的变化。
本文介绍的是基于 WebGPU Compute Shader 实现的 欧拉网格流体 (Eulerian Grid Fluid)。这种方法非常适合模拟烟雾、火焰等气体效果。
核心原理:Navier-Stokes 方程
流体运动由 Navier-Stokes 方程描述。对于不可压缩流体(如水、低速烟雾),我们主要关注两个方程:
动量方程 (Momentum Equation):
简单来说,速度的变化由以下几部分组成:
- : 平流 (Advection),流体顺着自己的速度移动。
- : 压力 (Pressure),流体从高压区流向低压区。
- : 扩散 (Diffusion),粘性导致的动量扩散(在烟雾模拟中通常忽略或由数值误差代替)。
- : 外力 (External Forces),如重力、浮力、人为施加的力。
不可压缩条件 (Incompressibility Condition):
即速度场的散度 (Divergence) 为零。流体不会凭空产生或消失,流入一个格子的量必须等于流出的量。
代码实现步骤
我们的实现使用了 WebGPU 的计算着色器 (Compute Shader),利用 GPU 的并行能力处理 512x512 甚至更高分辨率的网格。
渲染管线详解 (Render Pipeline Breakdown)
我们的 frame() 函数每一帧都会按顺序执行以下 Compute Pass。所有的物理量都存储在 rgba32float 纹理中,并使用 Ping-Pong (双缓冲) 技术进行读写交换。
1. 速度平流 (Advect Velocity)
- Pipeline:
advectPipeline - 资源 (RTs):
- Input:
velocityA(采样速度场),velocityA(被搬运的量) - Output:
velocityB
- Input:
- Swap:
velocityAvelocityB - 原理: 速度场自己搬运自己。我们根据当前速度回溯采样,计算下一时刻的速度分布。
let backUV = uv - vel * dt; let newVel = bilerp(velocityTex, backUV);
2. 密度平流 (Advect Density)
- Pipeline:
advectPipeline - 资源 (RTs):
- Input:
velocityA(采样速度场),densityA(被搬运的颜色) - Output:
densityB
- Input:
- Swap:
densityAdensityB - 原理: 烟雾(颜色)随着速度场移动。
3. 注入速度 (Splat Velocity)
- Pipeline:
splatPipeline - 资源 (RTs):
- Input:
velocityA - Output:
velocityB
- Input:
- Swap:
velocityAvelocityB - 原理: 当鼠标移动时,向速度场中叠加额外的速度向量。
4. 涡度约束 (Vorticity Confinement)
为了保持流体的旋转细节,防止数值耗散导致流体变“粘”,我们需要两个 Pass:
Pass 4.1: 计算旋度 (Compute Curl)
- Pipeline:
curlPipeline - 资源 (RTs):
- Input:
velocityA - Output:
curlTex(临时纹理,无需 Swap)
- Input:
- 原理: 计算每个点的旋转强度 。
Pass 4.2: 施加涡度力 (Apply Force)
- Pipeline:
vorticityPipeline - 资源 (RTs):
- Input:
velocityA,curlTex - Output:
velocityB
- Input:
- Swap:
velocityAvelocityB - 原理: 根据旋度梯度的方向计算增强力,叠加到速度场上。
5. 压力投影 (Pressure Projection)
这是确保流体“不可压缩”的关键步骤,包含三个阶段:
Pass 5.1: 计算散度 (Divergence)
- Pipeline:
divergencePipeline - 资源 (RTs):
- Input:
velocityA - Output:
divergence(临时纹理)
- Input:
- 原理: 计算速度场哪里“多了”或“少了”流体。
Pass 5.2: 求解压力 (Jacobi Iteration)
- Pipeline:
pressurePipeline(循环 20 次) - 资源 (RTs):
- Input:
pressureA,divergence - Output:
pressureB
- Input:
- Swap: 每次迭代
pressureApressureB - 原理: 求解泊松方程 ,计算出能抵消散度的压力场。
Pass 5.3: 减去梯度 (Gradient Subtract)
- Pipeline:
gradientSubtractPipeline - 资源 (RTs):
- Input:
pressureA,velocityA - Output:
velocityB
- Input:
- Swap:
velocityAvelocityB - 原理: ,修正速度场使其无散度。
6. 注入密度 (Splat Density)
- Pipeline:
splatPipeline - 资源 (RTs):
- Input:
densityA - Output:
densityB
- Input:
- Swap:
densityAdensityB - 原理: 在鼠标位置注入颜色(烟雾源)。
7. 渲染 (Render)
- Pipeline:
renderPipeline - 资源 (RTs):
- Input:
densityA - Output: Screen (Canvas)
- Input:
- 原理: 将最终的密度场绘制到屏幕上。
进阶:关于 Density 与水面渲染
densityA 代表什么?
在我们的模拟中,densityA(以及它的双缓冲伙伴 densityB)代表着流体中可见物质的密度或颜色。
- “被搬运的货物”:
Velocity(速度场) 是“运输工具”,Density(密度场) 是“货物”(如烟雾、墨水)。 - 屏幕像素:在渲染阶段,我们直接把
densityA纹理画到了屏幕上。开启“彩色流体”时存的是 RGB,关闭时存的是灰度。
如何改造成水面效果?
如果我们想模拟平静的水面被扰动的效果,不需要改变核心流体解算代码,只需要修改 渲染 (Render) 逻辑:
- 数据含义改变:把
densityA当作 波高 (Height),而不是颜色。 - 渲染 Shader 修改:
- 法线计算: 根据
densityA的梯度(邻域高度差)计算法线。 - 光照计算: 使用 Blinn-Phong 模型计算高光(Specular),模拟水面的波光粼粼。
- 泡沫 (Foam): 利用
Velocity场。如果length(velocity)很大,说明水流湍急,叠加白色泡沫。
- 法线计算: 根据
这样,原本飘动的“烟雾”就会变成水面上起伏的“波浪”。我们在下方的演示中增加了一个同步的水面渲染视口来展示这一效果。
总结
整个模拟循环如下:
- Advect Velocity: 让速度场自己搬运自己。
- Advect Density: 让速度场搬运颜色。
- Splat: 注入鼠标干扰。
- Vorticity: 增强旋转细节。
- Project:
- Calc Divergence
- Solve Pressure (Loop)
- Subtract Gradient
- Render: 显示结果。
通过 WebGPU,我们可以在浏览器中轻松实现 512x512 分辨率下 60FPS 的流畅模拟,这在 CPU 上是难以想象的。