首页 文章

# 基于 LBM (Lattice Boltzmann Method) 流体模拟

相对于基于解 N-S 方程的流体模拟方法，Lattice Boltzmann Method，简称 LBM，LBM 是另外一种截然不同的方法，这种方法和传统解 N-S 方程关系不同，传统方法是纯数值求解，通过逼近域来求解偏微分方程，但大多数求解都局限在二阶精度。LBM 方法避免了传统方法求解流体力学方程的复杂求解及低精度，但也能实现逼真的流体模拟动画，因为不同复杂求解 N-S 方程，所以速度上会快很多，在 3070s 显卡 1024 分辨率仅需 0.06ms，但缺点是显存占用比较多，因为要存 9 个函数，还有对于静止流体计算效率不高，不适用于强压缩性及，不适合在逼真粘度下直接模拟声音的远距离传播，优点就是 LBM 很容易实现，低噪声，出来的图像新闻上的热效应图，合适大型求解。

# LBM 方程

LBM 方法本质上是一种基于拉格朗日视角的模拟方法，形式上却跟欧拉视角下的网格有点相似，它将微观粒子和宏观规律相结合，同时有着较好的精度，在游戏领域上一般都是结算不可压缩流体，所以可以很好的达到需求而且速度也快。

LBM 是先把流体空间分割为一个个平均的网格，但同时网格内又存在通过概率密度函数来描述微观尺度下的粒子分布情况

Untitled

这种看似欧拉网格一样的划分，但还是有点不同，以二维为例，欧拉网格的邻居格子只有上下左右，但 LBM 的邻居是周围一圈 8 个格子，通常用 DdQq 来表示，d 表示维度，q 则表示格子数目，常用的有 D1Q3 、 D2Q9 、 D3Q15 、 D3Q19 和 D3Q27

每个格子都有个概率密度函数来描述粒子分布，这个函数就是 f（x,t），x 就是格子的位置，t 就是时间，每个格子有 9 个流动方向，一个速度方向有一个分布函数，所以每个格子有 9 个分布函数，

有了这些分布函数，就可以这样描述速度场和密度

\begin{align} \rho(x,t)&=\Sigma_i f_i(x,t)\\ u(x,t)&=\frac{1}{\rho(x,t)}\Sigma_i c_i f_i(x,t), \ \ i=0,...,8 \end{align} \tag {1}

其中 p（x，t）是虚拟粒子的密度，就是 9 个函数求和，u（x,t）就是速度，Ci 就是 9 个方向。

LBM 的核心分为碰撞和流动两步骤，这里碰撞就用 BGK 碰撞算子 Ωi (f), 这个算子经过 BGK 模型动力学演化过来的

Ωi(f)=−fi(x,t)−fieq(x,t)τffi(x,t+Δt)=fi(x,t)+Ωi(f)fi(x,t+Δt)=fi(x,t)−fi(x,t)−fieq(x,t)τf(2)Ωi(f) = -\frac{f_i(x, t)-f_i^{eq}(x,t)} {\tau_f} \\f_i(x, t+\Delta t)=f_i(x, t) + Ωi(f)\\f_i(x, t+\Delta t)=f_i(x, t) -\frac{f_i(x, t)-f_i^{eq}(x,t)} {\tau_f} \tag {2}

其中 τf 是松弛时间，和粘度 Mu 有关系，而 Fqe 函数就是流场平衡分布函数

μ=cs2(τf−Δt2)fieq(x,t)=ωiρ(x,t)(1+u(x,t)⋅cics2+(u(x,t)⋅ci)22cs4−u(x,t)⋅u(x,t)2cs2)(3)\mu =c_s^2(\tau_f -\frac{\Delta t}2)\\f_i^{eq}(x,t)=\omega_i \rho(x,t)(1+\frac{u(x,t)\cdot c_i}{c_s^2}+\frac{(u(x,t)\cdot c_i)^2}{2c_s^4}-\frac{u(x,t)\cdot u(x,t)}{2c_s^2}) \tag {3}

其中 C 就是声波在流体的传播速度，wi 为每个方向的权重，Cs = C / 根号 3，为了方便计算我们就设置格子的声速为 1。除 Cs 就是乘他的倒数，那化简后为

fieq(x,t)=ωiρ(x,t)(1+3∗u(x,t)⋅ci+4.5∗(u(x,t)⋅ci)2−1.5∗u(x,t)⋅u(x,t))(4)f_i^{eq}(x,t)=\omega_i \rho(x,t)(1+3*{u(x,t)\cdot c_i}+4.5*{(u(x,t)\cdot c_i)^2}-1.5*{u(x,t)\cdot u(x,t)}) \tag {4}

・是 dot，不是乘。接下来看一下流动步骤的核心公式

fi(x+ciΔt,t+Δt)=fi(x,t)(5)f_i(x+c_i \Delta t, t+\Delta t) = f_i(x,t) \tag {5}

这公式理解就非常简单了，就是相邻格子之间传播 F 函数

整理一下，那么每个格子只要解算 9 次这公式就可以了，然后把他们加起来就可以得到密度和速度了，是不是非常简单

fi(x,t+Δt)=fi(x−ciΔt,t)−Δtτf(fi(x−ciΔt,t)−fieq(x−ciΔt,t))(6)f_i(x, t+\Delta t)=f_i(x-c_i \Delta t,t)-\frac{\Delta t}{\tau_f}(f_i(x-c_i \Delta t,t)-f_i^{eq}(x-c_i \Delta t,t)) \tag {6}

然后整理下，提取下公因式就是

fi(x,t+Δt)=（1−Δtτf）∗fi(x−ciΔt）+Δtτf∗fieq(x−ciΔt,t))(6)f_i(x, t+\Delta t)=（1-\frac{\Delta t}{\tau_f}）*f_i(x-c_i \Delta t）+\frac{\Delta t}{\tau_f}*f_i^{eq}(x-c_i \Delta t,t)) \tag {6}

这样看是不是更简单了，就是说平衡函数和分布函数做个线性插值

# Niagara 实现

说完公式了，就下来说说迭代步骤：

  1. 解算密度和速度场，第一帧的话密度等于 1，速度定于定值。
  2. 解算平衡函数 Feq
  3. 每个方向的 Feq 和当前 F 函数做 lerp，松弛系数为 alpha，然后这个值就是新的 F 存起来

直接上 Niagara，因为解算过于简单，直接放图和 Hlsl

首先是初始化

https://pic2.zhimg.com/80/v2-faf1abff3655ee971b10491a0cc520d9_1440w.webp

https://pic3.zhimg.com/80/v2-36825ec7bcc17d53ecf743ee357b4572_1440w.webp

注意这里 Execute Behavior 要选 On Simulation Reset，才是第一执行

""
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
**int** XIndex;
**int** YIndex;
Grid2D.ExecutionIndexToGridIndex(XIndex, YIndex);

float2 Velocity **=** 0;
**float** Dentisy **=** 1.f;
**float** Solid **=** 0.f;
**if**(distance(float2((**float**)XIndex,(**float**)YIndex),CylinderPosition)**<**CylinderRadius)
{
    Solid **=** 1.f;
    Velocity **=** IniVelocity;
}
*/*if(distance(float2((float)XIndex,(float)YIndex),float2(300.f,300.f))<50.f)
{
    Solid = 1.f;
}*/*Grid2D.SetValueAtIndex(XIndex, YIndex, DensityIndex, Dentisy);
Grid2D.SetValueAtIndex(XIndex, YIndex, VelocityIndex, Velocity.x);
Grid2D.SetValueAtIndex(XIndex, YIndex, VelocityIndex**+**1, Velocity.y);*//SolidIndex*Grid2D.SetValueAtIndex(XIndex, YIndex, SolidIndex,Solid);

float2 Dir[9] **=** {float2(0, 0), float2(1, 0), float2(0, 1), 
float2(**-**1, 0), float2(0, **-**1), 
float2(1, 1),float2(**-**1, 1), 
float2(**-**1, **-**1), float2(1, **-**1)};

**float** W[9] **=** {4.0 **/** 9.0, 1.0 **/** 9.0, 1.0 **/** 9.0, 1.0 **/** 9.0, 1.0 **/** 9.0, 1.0 **/** 36.0,
        1.0 **/** 36.0, 1.0 **/** 36.0, 1.0 **/** 36.0};

**for**(**int** index **=** 0;index**<**9;index**++**)
{
   
    **float** eu **=** dot(Dir[index],Velocity);
    **float** Dotvel **=** dot(Velocity,Velocity);
    
    **float** balance **=** W[index]*****Dentisy*****(1.f**+**3.f*****eu**+**4.5f*****eu*****eu**-**1.5f*****Dotvel);
    FunctionsGrid.SetValueAtIndex(XIndex, YIndex, 0**+**index, balance);
}

这里我们把周围虚拟格子 function 用 index 写入 grid2d，还有速度和密度

然后就是核心主要解算部分

https://pic3.zhimg.com/80/v2-657a88ab25701ccbf03631cfa1b15312_1440w.webp

""
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
int XIndex;
int YIndex;
int MaxXCell;
int MaxYCell;
Grid2D.ExecutionIndexToGridIndex(XIndex, YIndex);
Grid2D.GetNumCells(MaxXCell,MaxYCell);
float2 Velocity = 0;
float Dentisy = 0.f;

//Grid2D.GetPreviousValueAtIndex(XIndex, YIndex, DensityIndex, Dentisy);

float2 Dir[9] = {float2(0, 0), float2(1, 0), float2(0, 1), 
float2(-1, 0), float2(0, -1), 
float2(1, 1),float2(-1, 1), 
float2(-1, -1), float2(1, -1)};

float W[9] = {4.0 / 9.0, 1.0 / 9.0, 1.0 / 9.0, 1.0 / 9.0, 1.0 / 9.0, 1.0 / 36.0,
        1.0 / 36.0, 1.0 / 36.0, 1.0 / 36.0};
float Solid;

float Alpha = Inv_tau;
Grid2D.GetPreviousValueAtIndex(XIndex, YIndex, SolidIndex, Solid);
Grid2D.SetValueAtIndex(XIndex, YIndex, SolidIndex, Solid);
for(int index = 0;index<9;index++)
{
    int NewIndexX = XIndex+(Dir[index].x);
    int NewIndexY = YIndex+(Dir[index].y);
    if(NewIndexX>MaxXCell-1)
    {
        NewIndexX = 0;
    }
    if(NewIndexY>MaxYCell-1)
    {
        NewIndexY =0;
    }
    if(NewIndexX<0)
    {
        NewIndexX = NewIndexX - 1;
    }
    if(NewIndexY<0)
    {
        NewIndexY = NewIndexX - 1;
    }    
    float OtherBalance;
    FunctionsGrid.GetPreviousValueAtIndex(NewIndexX, NewIndexY, index, OtherBalance);
    Dentisy += OtherBalance;
    Velocity += Dir[index]*OtherBalance;       
}
Velocity /=Dentisy;

if(length(Velocity)>0.5)
{
    Velocity*= 0.5/length(Velocity);
}

if(XIndex==0||XIndex==MaxXCell-1)
{
    Dentisy = 1;
    Velocity =CollsionVelocity;
    Alpha = 1;
}

if(Solid>0.5)
{
    Dentisy = 1;
    Velocity = float2(0.0,0.0);
    Alpha = 1;
}
    
for(int index = 0;index<9;index++)
{
    
    int NewIndexX = XIndex+(Dir[index].x);
    int NewIndexY = YIndex+(Dir[index].y);
    if(NewIndexX>MaxXCell-1)
    {
        NewIndexX = 0;
    }
    if(NewIndexY>MaxYCell-1)
    {
        NewIndexY = 0;
    }
    if(NewIndexX<0)
    {
        NewIndexX = NewIndexX - 1;
    }
    if(NewIndexY<0)
    {
        NewIndexY = MaxYCell-1;
    }
    float CurrentBalance;
    FunctionsGrid.GetPreviousValueAtIndex(NewIndexX, NewIndexY, index, CurrentBalance);
    
    float eu = dot(Dir[index],Velocity);
    float Dotvel = dot(Velocity,Velocity);
    
    //float CurrentBalance;
    float balance = W[index]*Dentisy*(1.f+3.f*eu+4.5f*eu*eu-1.5f*Dotvel);

    CurrentBalance = lerp(CurrentBalance,balance,Alpha);
    FunctionsGrid.SetValueAtIndex(XIndex, YIndex, index, CurrentBalance);
}

Grid2D.SetValueAtIndex(XIndex, YIndex, DensityIndex, Dentisy);
Grid2D.SetValueAtIndex(XIndex, YIndex, VelocityIndex, Velocity.x);
Grid2D.SetValueAtIndex(XIndex, YIndex, VelocityIndex+1, Velocity.y);

最后在输出到 RT 上，绑定到材质上，粒子这时候就可以渲染了

https://pic2.zhimg.com/80/v2-0cec0c0663e7ff416123c09f382b77b5_1440w.webp

最后测试下来，只需要 0.22ms，比直接解 NS 公式要快 10 倍，但缺点是只能模拟简单的效果并没有 NS 那么精细，可能模拟河流这些还是非常好的

https://pic1.zhimg.com/80/v2-c228ff86c7330d0f81bdce12b95f258c_1440w.webp