第六讲：流体模拟- 流体力学基础

2022-06-11 本文已影响0人离原春草

虽然绿灯没怎么为我亮过，但我还是对生活充满热情，这就是我理解的年轻 —— 汪曾祺《人间草木》

所谓的流体就是流动的气态、液态甚至固态的物体，比如空气、水、沙子等，这是经典物理里面目前还依然没能成功攻克的难题。

普通物体的形变可以分成三种：

拉伸Stretch
压缩Compression
剪切形变Shear

而流体一个方面很难做到前面两种形变，另一个方面则十分容易发生第三种形变，因为这个性质的存在，导致计算机模拟变得很困难，因为：

前者意味着当受力之后由于抵抗拉伸与压缩的能力很强，这个受力状态将会以扰动的方式很快的传播开来，也就是说任何一个扰动将会导致整个流体中的所有位置都会受到影响（即使应用了前面的Projection Dynamics，也很难计算出结果）
后者意味着扰动会容易引起旋涡（湍流），而旋涡的特性在于大旋涡会容易变成小旋涡，即旋涡会不断分解，从而导致使用离散化方式计算的计算机模拟无法完全模拟出流体动力学的细节。

流体问题可以从空间跟时间两个角度来看待，尺度有：

空间：从大到小，分子->平均自由程（分子相撞之前平均的移动距离）->system scale（空气流过飞机翅膀经过的距离）
时间：t0(分子碰撞时间）-> t1（水分子布朗运动碰撞前的平均经历时间）->tc（空气流过飞机翅膀需要的时间）

根据不同的尺度，有不同的模拟方式：

分子尺度：分子动力学（合成、制药）
自由程尺度：统计意义上的效果，Kinetic theory，这是最近一段时间开始流行的模拟角度，Lattice Boltzmann Model
大尺度（system scale）：NS方程，这是计算机模拟在过去二三十年专注的区域

游戏模拟中常用的方式是后面两种，这两种模拟方式各有千秋，这里会优先介绍大尺度下的模拟方式，后面有时间再介绍Kinetic theory方法。

宏观尺度下流体模拟的一个核心思想是Field View，即场视角，也就是说对于流体中任意一点的物理量（如速度、温度、密度、压强等）并不是固定的，而是与对应点位置相关的一个数值。

NS方程本质上是牛顿第二定律：

$f = m \frac{d u}{d t}$

在流体中，速度是一个向量场 $u(x, t)$ ，数值随时间、位置发生变化。

在流体中，我们对速度u进行求导，得到的不是加速度（回顾下加速度的定义，单位时间内速度的变化量），举个例子，某个小船在水面上行驶，其加速度可以按照下面方式进行推导。

首先原始的速度为 $u_0 = u(t_0, x_0)$ ，经过 $\Delta t$ 之后，其速度就变成了 $u(t_0 + \Delta t, x_0 + u_0 \Delta t)$ ，加速度就变成了：

$a = \lim_{\Delta t \rightarrow 0}\frac{u(t_0 + \Delta t, x_0 + u_0 \Delta t) - u_0}{\Delta t}$

用泰勒展开就变成了：

$a = \frac{D u}{D t} = u \nabla_x u + \frac{\partial u}{\partial t}$

这个加速度叫Material Derivative，这个公式中，前面一项表示的是时间不变，空间上的速度变化，后面一项则是空间不变（即小船待在原地不动）随着时间流逝带来的速度的变化。

加速度有了，那么质量怎么来定义呢，流体没有单个点的只能，只能用单位体积的质量（也就是密度 $\rho$ ）来表示，那么我们有：

$\frac{D u}{D t} * \rho = \frac{f}{V}$

$u \nabla_x u * \rho + \frac{\partial u}{\partial t} * \rho = \frac{f}{V}$

上面的V是流体的体积（ $m = \rho * V$ )，而 $\frac{f}{V}$ 可以看成是力的密度：

对于重力而言，其密度就是单位质量上的力，也就是 $\rho g = \frac{m g}{V}$
压力，以2D为例，在流体中任意框选出一个box中，四边都有压强，且左右上下压强不一定相同，从而形成一个力，这个力就是压强的梯度： $- \nabla P$
viscosity（抵抗剪切形变的力），这一项就直接给出公式为 $\mu \nabla \cdot \nabla u$ ，被称为viscosity term，其中的 $\mu$ 称为dynamic viscosity。

在上面这几项力的作用下，我们就有最终的NS方程：

$\frac{\partial u}{\partial t} * \rho+ \rho(u \nabla_x) u = \rho g - \nabla P + \mu \nabla \cdot \nabla u$

在很多流体模拟中， $\mu$ 很小趋近于0，所以通常直接省略最后面一项。

这个方程中，u是不知道的，压强P也是不知道的，密度也是不知道的，但是只有一个方程，因此是不足以求解的，要想求解，还需要补上两个方程。

根据质量守恒方程，我们有（难受，这里没听到，看了半天也没看出来是如何推导的，后面有机会再补上对应的说明）：

$\frac{d}{dt}\int_{V_0}\rho(x) dx = - \oint \rho u dA$

这个公式可以变化成（也错过了这个变化的推导，-_-||）：

$\int_{V_0} \frac{d}{dt}\rho(x) dx = - \int_{V_0} \nabla(\rho u) dx$

最后我们将两者合并，得到：

$\int_{V_0} [\frac{d}{dt}\rho(x) + \nabla(\rho u)] dx = 0$

由于质量不能为负，所以我们有：

$\frac{d}{dt}\rho(x) + \nabla(\rho u) = 0$

也就是：

$\frac{d}{dt}\rho(x) + u \nabla \rho + \rho \nabla u = 0$

而等式左边前面两部分（为什么？参考前面提到的Material Derivative）又可以表示成： $\frac{D \rho}{D t}$ ，所以我们有：

$\frac{D \rho}{D t} + \rho \nabla u = 0$

而基于流体的不可压缩性原理，我们可以得到速度的梯度为0：

$\nabla u = 0$

好了，有了上面这些条件，我们就可以着手来求解NS方程了，这里介绍一种经典的解法（也就是Houndini中的解法），即通过time integrator的方式来求解，对NS方程进行分析，我们可以知道，整个方程中只有第一项 $\frac{\partial u}{\partial t}$ 是跟时间有关系的，也就是说，我们可以将这一项按照时间写成如下形式：

$\frac{\partial u}{\partial t} \approx \frac{u^{n+1} - u^n}{\Delta t}$

接下来只需要将NS方程中剩下的几项估计出来，就能求得 $u^{n+1}$ ，但是由于这里有多项需要估计，且压强是不可知的，我们要如何进行估计与求解呢？

这里常用的一种方法叫做Operator Splitting，简单来说，如果我们有一个简单的微分方程：

$\frac{dx}{dt} = f(x) + g(x)$

根据显示欧拉方程，我们转换一下就是：

$x^{n+1} = x^n + f(x^n) \Delta t+ g(x^n) \Delta t$

这个公式可以将之拆分成两步来求得：

$x^* = x^n + f(x^n) \Delta t$

$x^{n+1} = x^* + g(x^n) \Delta t$

我们知道，显式欧拉跟隐式欧拉之间的区别在于右边的函数使用的x是 $x^n$ 还是 $x^{n+1}$ ，其实都是一种近似，而在这里我们也可以用其他的x来近似，所以这里我们可以将上面的第二个式子改写成：

$x^{n+1} \approx x^* + g(x^*) \Delta t$

所以问题就变成了，我们需要求解两个简单的微分方程：

$\frac{dx}{dt} = f(x)$

$\frac{dx}{dt} = g(x)$

回到NS方程来说，我们就可以采用同样的思想将之拆解成四个步骤的简单微分方程的求解（实际上，最后一项粘稠性非常小，可以先暂时忽略）。

首先，我们先看左边的一项：

$\frac{\partial u}{\partial t} \rho = -u \nabla u$

右边的 $\nabla u$ 是u的梯度，也就是空间中的微分，一种最简单的方式（经典软件使用的方法会更复杂）就是将一块区域想象成由格子来定义的场，每个采样点都落在格子的顶点上，那么通过 $\frac{u_{n+1} - u_n}{\Delta s}$ 就能求解出 $\nabla u$

之后就可以用显式欧拉的方法来求解出 $u^{n+1}$ 了，不过由于显式欧拉本身不稳定，会不断震荡跑飞，而如果用隐式欧拉，则需要求解非线性方程，加大了复杂度。

而上面这个方程其实可以表示成：

$\frac{Du}{Dt} = 0$

可以想象有个particle在流体中流动，如上图中所示，从格子上顶点流动到右上角的红点上，也就是说，上面公式的物理意义在于，在经过这个流动之后，两者的速度并不会发生变化，也就是说我们可以采用一种叫做Semi-Lagrangian的方法来求解。

以上图进行说明，假设我们想要知道上图中绿色方框中的一点的速度，我们假设这一点是来自于图上五角星位置，那么我们就可以知道上一个timestep上五角星位置跟当前timestep上绿色box位置的速度是相同的，这里有两个问题：

如何从绿色方框位置反推出五角星位置，这就是semi的由来，也就是根据当前绿色方框的位置朝着时间相反的方向以当前的速度移动一段距离，就知道五角星位置
找到五角星位置后，我们怎么知道五角星位置的速度；由于我们知道顶点上的速度，我们就可以通过插值得到五角星位置的速度了。

这种算法的好处在于：

整个系统是稳定的，因为我们只是将某一点的速度拷贝到另一点，而不会出现速度的增加，所以不会爆掉
求解方便，不需要求解非线性方程
Energy Disappearance会小一些，隐式欧拉会随着时间的推动，能量会不断磨损，而如果这个过程过于迅速的话，流体的效果就一闪而逝了，看起来就不真实，而这里的Semi-Lagrangian方法就可以稍微降低消失的速度，使得效果更好一点。

经过这一步计算之后，我们就根据 $u^n$ 得到了 $u^*$ ，接下来我们考虑右边的第一个式子：

$\frac{\partial u}{\partial t} = g$

由于右边是一个常量，因此这个公式可以直接使用显式欧拉求解得到，从而得到了 $\tilde u$ ，接着再看下一个式子：

$\frac{\partial u}{\partial t} = -\frac{1}{\rho} \nabla P$

由于我们不知道压强P，所以需要将流体的不可压缩特性用起来，也就是说，用近似公式，我们有：

$u^{n+1} = \tilde u -\frac{1}{\rho} \nabla P \Delta t$

而由于 $\nabla u = 0$ （注意，这里的梯度是空间上的差分或者微分），也就是：

$\nabla u^{n+1} = \nabla \tilde u -\frac{\Delta t}{\rho} \nabla \cdot \nabla P = 0$

根据这个式子我们有：
$\frac{\rho \cdot \nabla \tilde u}{\Delta t} = \nabla \cdot \nabla P$

这个公式叫做Pressure Projection，也就是说，我们有一个压强P，这个参数存在于上图中的格子的顶点上，而这个数据我们就可以通过有限差分公式估计出其二阶导数，而由于这个方程中我们不知道压强P，而密度 $\rho$ 是常量（不可压缩）， $\Delta t$ 也是已知常量，而 $\tilde u$ 是上一步计算得到的结果，所以我们就可以求解得到P（线性方程组求解），将这个代入前面的公式，我们就能够得到 $u^{n+1}$ 。

$u^{n+1} = \tilde u -\frac{1}{\rho} \nabla P \Delta t$

这了需要注意的是，我们用P来估计 $u^{n+1}$ ，采用的是 $\nabla P$ ，而这个有限差分估计出来的结果实际上是上图中顶点中心处的数值，用这个数值估算出来的 $u^{n+1}$ 也就是中心位置的速度，这就是为什么所有的流体算法将上图中的格子叫做Staggered Grid，也就是速度是存放在格子的每条边的中点位置上的。

求解就是这个过程，这里需要说几点：

$u^*$ 跟 $\tilde u$ 的求解是解耦的，因此可以并行完成，且可以放在GPU上来实现（每个格点是独立的）
最后P的求解是需要求解一个线性方程组，需要将整个Grid上的每个格点都耦合在一起，因为流体的不可压缩性，因此在任意位置给流体一个压强或者扰动，就会在无限短的时间内导致全盘的影响，因此需要求解线性方程组，这个计算会很麻烦；目前经典解算方法都是采用这种方式，比如Houndini，虽然可以通过GPU来甲酸，但是因为不太好并行，是迭代的过程，在这一步就会很慢。
弹性体的时候，我们说过，一个约束条件下的求解，可以看成是一个优化问题；这里我们可以采用同样的方式来考虑这个问题，假设我们没有流体的不可压缩条件，那么我们就可以忽视NS方程中的P那一项，也就是说我们只需要经过前面两步，不考虑最后P的求解即可，但实际上不可压缩性是实际存在的，因此我们需要求解如下优化方程：

$min \frac{1}{2}(u^{n+1} - \tilde u)\rho(u^{n+1} - \tilde u)$

也就是说 $(u^{n+1} - \tilde u)$ 最小的时候我们得到最优解，但是由于不可压缩性，我们需要满足：

$\nabla u^{n+1} = 0$

而这种条件我们就需要使用Lagrangian Multiplier方法来求解，其求解得到的结果跟前面Operator Splitting得到的结论是一样的，从这个角度来看，压强P不是一个真正的物理量，而是一个力，其作用是使得上面的流体具有不可压缩性，也就是说用来反抗流体中的每一点的速度变化的（惯性力）。

另外需要介绍的是流体在一个容器内流动的边界条件：

边界上每一点的法线方向上的速度为0（从而避免流体穿过容器）
如果流体中有一个扰动物，比如一个叶片，那么叶片上的速度与叶片上流体的速度是一样的

下面来介绍一下游戏中常用的SPH（Smoothed Particle Hydrodynamics，最开始来自于一些核试验模拟领域）算法，因为NS中的线性方程组的求解很难做到实时，而很多电影特效中常用的就是NS Solver方法（为什么，因为精度更高吗？）

SPH（Smoothed Particle Hydrodynamics）

SPH的思想是，我们将流体看成是一堆particle，想象particle可以在其中移动，我们可以为每个particle赋予一套物理量，假设我们称之为 $f_i$

假设我们知道了每个粒子上的物理量，那么对于任意一个位置x，我们就需要能计算出其物理量 $f_x$ ，从而计算得到其gradient，甚至二阶gradient（Laplacian），我们就能将这些数据代入NS方程中，完成粒子上速度的更新，从而用这个速度对粒子进行移动，就能得到下一个timestep的粒子。

也就是说，其中心思想是根据粒子的物理量，求解出整个场域的物理量。

这个算法中最核心的东西就是一个 $\delta(x)$ ，这不是一个函数，而是一个分布：

在x=0处有一个非零值，其他位置结果为0，定义可以表示为：

$f(x) = \int f(y) \delta(y-x) dx$

实际上这个分布是一个理想的分布，现实中是没有的，也不好给出表达式，而我们可以近似一下，给出如下的近似分布：

$f(x) \approx \int f(y) W(|y-x|,h) dx$

这里的W称之为kernel函数，h用来控制上面近似分布的宽度，h越大，分布越宽。

有了这个近似的分布函数，我们就可以表达f(x)，也就是对于任何的x，我们可以将f(x)表达为一个求和项：

$f(x) = \sum_i f_i W(|x-x_i|, h) w_i$

最后的 $w_i$ 是权重，而W是kernel函数，我们是知道其表达式的，而这个函数在|y-x|大于某个范围的时候，这个函数返回为0。

这里我们先不管f表示的是什么，就先看成是任意的物理量，我们推导出任意位置的f(x)， $\nabla f$ 跟 $\nabla \cdot \nabla f$ ，就可以将实际的物理量（比如速度、压强）代入进去，并放在NS方程中，就能得到 $u^{n+1}$ ，并根据这个速度对粒子进行移动，就能得到下一个timestep的初始数据，完成了一轮迭代。

令r = y - x，我们可以给出W的一个示例函数：

$W(r,h)=\sigma \begin{cases} 6(q^3-q^2)+1 & {0 \leq q \leq \frac{1}{2}}\\ 2(1-q)^3 & { \frac{1}{2} < q \leq 1}\\ 0 & {otherwise} \end{cases}$

其中 $q=\frac{|r|}{h}$ 。

上面已经给出了f(x)的求和项，那么其gradient就可以在前面公式两边加一个gradient就行，那么由于 $f_i$ 是常量， $w_i$ 也是不变的，那么gradient就施加在W上，但是这样的做法虽然形式上没问题，实际上我们在执行中就会发现运行出来的结果就会爆掉（表现跟显式欧拉一样），这是因为粒子在移动过程中，有些地方粒子会密集一点，有些会稀疏。

在密集的位置，这个积分变成求和的近似是相对准确的，但是稀疏区域的近似就不准确，就会导致数值上的震荡，即误差会不断累计，导致爆掉。

所以，这个给我们的启发就是在密集区域，可以按照这种方式来求取，但是在稀疏的区域就需要另外处理，也就是首先我们得知道各个位置的粒子密度，我们需要先计算 $\nabla (\frac{f}{\rho})$ ，也就是用粒子的密度来对物理量进行归一化。

假设我们定义粒子的密度为：

$\rho_i = \sum_j m_j W_{ij}$

其中 $m_j$ 表示的是j点的质量， $W_{ij} = W(|x_i - x_j|,h)$

那么用链式法则（复合函数的微分）展开就可以得到：

$\nabla (\frac{f}{\rho}) = \frac{\nabla f}{\rho} - \frac{f}{\rho^2}\nabla \rho$

变换一下就得到了：

$\nabla f = \rho \nabla (\frac{f}{\rho}) + \frac{f}{\rho}\nabla \rho = \rho \cdot (\nabla (\frac{f}{\rho}) + \frac{f}{\rho^2}\nabla \rho)$

右边有两个gradient，我们就可以用前面的求和（可能导致不稳定的）公式分别来近似这两个式子

$\nabla (\frac{f}{\rho}) = \sum_j \frac{f_j}{\rho_j^2} \nabla W_{j} w_j \rho_j = \sum_j \frac{f_j}{\rho_j^2} \nabla W_{j} m_j$

$w_j \rho_j = m_j$ 表示的是某个位置的质量，可以用这个位置的密度乘上其权重算得。

$\frac{f}{\rho^2}\nabla \rho = \sum_j m_j \frac{f}{\rho^2} \nabla W_{j}$

最后得到一个表达式(定义在任意一个粒子上的位置，而非任意位置）：

$\nabla f(x_i) = \nabla f_i = \rho_i \sum_j m_j(\frac{f_j}{\rho_j^2} + \frac{f_i}{\rho_i^2})\nabla_i W_{ij}$

而这个式子则能够避免简单gradient导致的不稳定问题，因为用这个公式计算，可以保证动量是守恒的。对于上面的式子，每个粒子的动量之和有：

$\sum_i m_i u_i = \sum_i m_i(u_i -\frac{1}{\rho_i} \nabla P_i)$

上面式子中的压强来自于NS方程（除以了等号左边的密度）， $\frac{1}{\rho_i} \nabla P_i$ 表示的是由于压强（内力）导致的速度变化。因为如果流体的速度的变化是由于内力（即压强，准确来说是压强的梯度）导致的，那么就可以保证动量是守恒的。

移除同类项：

$\sum_i m_i \frac{1}{\rho_i} \nabla P_i = 0$

将前面 $\nabla f_i$ 的求和式子代入 $\nabla P_i$ 就有：

$\sum_i m_i \sum_j m_j (\frac{f_j}{\rho_j^2} + \frac{f_i}{\rho_i^2})\nabla_i W_{ij} = 0$

那么这个等式为什么成立呢？

同样，我们就可以求得Laplacian式子，这个有很多个版本：

$\nabla^2 f(x_i) = \sum_j \frac{m_j}{\rho_j} f_j \nabla^2 W_{ij}$

另一个常用的版本为：

$\nabla^2 f(x_i) \approx 2(d+2) \sum_j \frac{m_j}{\rho_j} \frac{A_{ij} \cdot r_{ij}}{|r_{ij}|^2} \nabla_i W_{ij}$

其中d表示的是维度，2维则取值为2，另：

$A_{ij} = \frac{f_j}{\rho_j^2} + \frac{f_i}{\rho_i^2}$

有了一阶二阶梯度公式之后，我们再来求解NS方程：

$\frac{\partial u}{\partial t} * \rho+ \rho(u \nabla_x) u = \rho g - \nabla P + \mu \nabla \cdot \nabla u$

根据Operator Splitting的思想，我们分步进行处理，首先：

$\frac{Du}{Dt} = \frac{\partial u}{\partial t}+ (u \nabla_x) u = \frac{\mu}{\rho} \nabla \cdot \nabla u + g$

根据前面的推导，右边的式子中我们都是知道的（ $\mu$ 是dynamic viscosity），我们就能够通过显式欧拉的方法完成速度的求解。此外，这个过程对于每个粒子而言是完全独立的，因此可以放在GPU上来计算。

之后根据 $\frac{D \rho}{Dt}=0$ 计算出 $\nabla P$ ，压强跟密度有很多种映射表达方式，这里选取最简单的一种：

$P = k (\rho - \rho_0)$

其中 $\rho_0$ 表示的是期望的密度，选取这种方式的好处在于各个粒子之间的计算是相互独立的方便GPU并行。

另外，这个公式也就说明了虽然我们一直在说流体的不可压缩性，但实际上如气体在空间中的流动，其密度是会发生变化的，而这个公式实际上相当于添加了一个惩罚项，目的是维持流体的不可压缩性。

当然，学术界也有其他的表达式是能确保流体的不可压缩性的，但是这类的方法会导致需要进行多轮迭代才能输出结果，从而打破了GPU的并行性，这也是SPH方法尴尬的地方，虽然是想要设计成GPU并行的，但是却没有办法真正并行。

再之后就求解：

$\frac{D u}{Dt} = -\frac{\nabla P}{\rho}$

并根据这个式子完成速度u的更新。

最后根据：

$\frac{D x}{Dt} = u$

求得更新后的位置。

下面来介绍一下Lattice Boltzman Method。

Lattice Boltzman Method

这个方法是上个世纪七十年代就已经提出来的方法，这个方法来自于Kinetic Theory（统计物理学），之所以这么久以来一直不温不火是因为这个方法不够精确，但是2006年左右出现了一系列的方法提升了这个方法的精确性。

而且这个方法本身没有迭代的过程，十分适合GPU并行计算。

统计学中的一个概念叫做Distribution Function（分布函数） $f(x, \xi, t)$ ，这个函数的三个参数分别对应位置、速度与时间，表示的是随着三个参数的不同，某种现象发生的概率，这里我们考虑的是流体中微观层面（如布朗运动）现象的概率，比如在任何时间t在位置x处观察到的以速度 $\xi$ 移动的空气分子的概率（或者密度）。

对这个函数沿着速度（三维）进行积分，就能得到当前时刻对应位置处的空气分子的密度（这是一个标量）：

$\int f(x, \xi, t) d^3 \xi = \rho(x,t)$

而这个公式可以称之为0阶动量（0-th momentum)，而一阶动量则可以表示为：

$\int \xi f(x, \xi, t) d^3 \xi = \rho u$

这里的u是宏观层面的速度（即NS方程中的速度），其实左边项可以看成是微观层面的速度乘以微观的速度（密度），积分之后就是宏观层面的动量了，这个结果是一个向量。

从上面两个公式也可以窥见，LBM方法的要点就是抛开宏观层面的计算，我们关心微观层面的概率分布函数是什么样的。

我们先来看下，概率分布函数的微分形式应该是什么样的：

$\frac{df}{dt} = \frac{\partial f}{\partial t} + \nabla_x f \cdot \frac{\partial x}{\partial t} + \nabla_{\xi} f \cdot \frac{\partial \xi}{\partial t} = \frac{\partial f}{\partial t} + \nabla_x f \cdot \xi + \nabla_{\xi} f \cdot F = \Omega(f)$

上面式子中的F是外力（其实是加速度，这里假设质量为1）， $\Omega$ 叫做Collision Operator，如果这个函数是知道的，那么后面这个等式被叫做Boltzman Equation，不过历史上大家是不知道这个公式是什么样的，因此大家会通过一些条件来推测一些经验公式。

质量守恒
$\int \Omega(f) d^3 \xi = 0$
动量守恒
$\int \xi \Omega(f) d^3 \xi = 0$
能量守恒
$\int |\xi|^2 \Omega(f) d^3 \xi = 0$

根据这三个条件，大家凑出一个模型，叫做BGK Model：
$\Omega(f) = -\frac{1}{\tau}(f - f^{eq}) = \frac{df}{dt}$

这个公式中 $\tau$ 是时间的常量，而 $f^{eq}$ 则是将气体在静止很长一段时间之后，温度达到平衡之后，对应的distribution function。

右边的等式就是一个常微分方程，手算之后其结果输出为：
$f(x, \xi, t) = f^{eq}+[f(x, \xi, 0)-f^{eq}]e^{-\frac{t}{\tau}}$

这个公式的意思是，随着时间的推移，概率将会以指数的形式递减，最终达到 $f^{eq}$ ，流体密度越大， $\tau$ 越小，越快达到平衡，也就是说， $\tau$ 控制的是流体的粘性。

BGK方法虽然简单，但是有一个比较显著的问题在于准确性与真实情况相去甚远，毕竟是拼凑出来的公式，后面大家观察到，这个不准确性来自于，我们这里将f看成是一个随着时间均匀趋向平衡的函数，但是实际上f代表的各个特性，如质量、动量以及能量（也就是不同阶数的梯度）其趋向平衡的速度是不同的。

基于上述发现，2006年之后出现了一系列新的方法大大提升了整个方案的精确性（当然也会变得更为复杂），时间关系，这里就不做扩展。

不管怎么说，我们都有了一个Collision Operator，接下来要做的事情就是根据timestep对其进行积分以求得f。

先来看一个叫做Hermite Interpolation的方法，以2D Domain为例，我们将流体场域分成Grid，每个顶点上对应的就是位置x，对于每个顶点而言，我们存储的不是一个数值，而是多个数值，每个数值对应的f在同一时刻下不同方向下的数值 $f(x, \xi_i, t)$ （即不同速度的概率？），这些方向分别对应于当前顶点指向相邻顶点（加上自身），总共9个方向，因此总结为D2Q9，在3D情况下，就是D3Q27等等。

下面我们看下，如何根据Collision Operator求得f。

$\frac{\partial f}{\partial t} + \nabla_x f \cdot \xi + \nabla_{\xi} f \cdot F = \Omega(f)$

同样，这里采用Splitting Operator。

先考虑前两项：

$\frac{\partial f}{\partial t} + \nabla_x f \cdot \xi = 0$

这个就是advection equation，在NS方程求解的时候，advection是速度场，当时是使用Semi-Lagrangian方法求解的，而这里更为简单，考虑刚才的格子D2Q9，如果我们只考虑一个方向的分量，按照Semi-Lagrangian方法，由于顶点上存储的指向各个相邻顶点的速度其长度正好在一个timestep上走到了相邻顶点上，因此对于相邻顶点在当前timestep的物理量就根本不用求解，只需要将此顶点上的物理量拷贝过去即可，从而就能完成这个方程的求解（文献中称之为Streaming）。

最后考虑剩下的部分，在刚才计算的结果上加上Collision Operator（这个是已知的，前面说过的经验模型）以及外力的作用，就得到了下面的式子

$\frac{\partial f}{\partial t} = \Omega(f) - \nabla_{\xi} f \cdot F$

这个方法由于在每个顶点上进行的计算都是完全独立的，因此适合在GPU上计算，适合用在实时场景。

第六讲：流体模拟- 流体力学基础

SPH（Smoothed Particle Hydrodynamics）

Lattice Boltzman Method

猜你喜欢

热点阅读