流体模拟《Fluid Simulation for Computer Graphics 》——从入门到放弃(四)烟雾模拟/涡度限制/表面追踪/水平集方法

从入门到放弃系列

学的都是个啥,,开始怀疑人生ಥ_ಥ

有些提到的东西可以在前置文章里找到,哪里理解错误啥的请指出3Q


开头再放上之前提到的Navier-Stokes方程和变量:

  • 向量u表示流体速度
  • ρ表示流体密度
  • p表示压力,1\ρ ∇p体现了压力对流体的影响
  • 向量g表示重力加速度
  • ν表示运动粘度
  • q表示通过流体传输的模拟量之一(比如温度)
  • t表示时间

7 烟雾

在流体是烟雾的情况下,模拟烟雾需要两个额外的变量:

  • 空气的温度T(摄氏度)
  • 烟雾颗粒的浓度s(0-1之间)

二者材料导数都为0。然后研究Ts对流体速度有什么影响,我们知道热空气会上升,冷空气会下沉;好像载有较重烟雾颗粒的空气在重力作用下会被往下拉。我们可以用之前那个动量方程,然后用浮力加速度代替原来的重力加速度g

其中αβ是非负系数,T_amb是环境温度,比如s=0T=T_amb的时候,这个就相当于环境里的空气,自然有b=0

7.1 烟雾来源

通常要定义一个烟雾源头区域,对这个区域每个时间间隔都要增加一点烟雾和热量。对烟雾源头的每个网格点的更新公式如下:

其中r_tr_s用于控制添加烟雾和热量的速率(在烟雾源头外这两个量都为0),T_target表示源头的温度。

7.2 热量和烟雾的扩散

热量和烟雾也可以扩散,在这种情况下,小规模的现象(比如传导和布朗运动)或是大规模的现象(比如紊流混合),用于平滑梯度。这可以用拉普拉斯算子建模:

其中k_Tk_s是非负扩散常数。

比如温度扩散就是下面这个式子:

7.3 vorticity confinement(涡度限制)

涡度限制是对Navier-Stokes方程的一种试图保留涡度的修改。引入Δx,这样在极限中涡度就没了。基本思想是检测旋涡的位置,并且增加一个力来满足每个旋涡周围的旋转运动,这样旋涡位置速度快于周围所有流体,可以简单的通过标准化𝛚的梯度来构造指向这些漩涡中心的单位向量N

现在N指向旋涡的旋转中心,𝛚方向沿旋转轴,所以要获得增加的那个旋转力,只要一个叉积:

其中ε是一个可以调整涡度限制效果的参数。

7.4 Divergence Control(发散控制)

通过使用这种浮力模型,假设流体密度是温度和烟雾浓度的函数(烟雾溶解在空气中)。在热膨胀问题中,流体应该在加热时膨胀在冷却时收缩,即在DT/Dt≠0时,我们不需要无散度的速度场。如果温度升高导致密度降低,那么我们需要正散度来实施这个热膨胀。可以通过在网格单元中心定义一个控制场d(x),它等于在整个流体体积中的体积分数变化率(ΔV/V)·Δt,然后通过压力来求解这种差异。

文中是在使用MAC网格的数值离散的压力求解部分,将d加到线性系统的流体单元的右手侧。

8 表面追踪

对于液体(比如水),需要找到一种方法来跟踪液体表面。这是通过使用隐式表面表示或者不同隐式表面组合来实现的。

8.1 Tracker particles(跟踪器粒子)

首先用粒子填充水体积并将其视为水的几何形状(如果存在水源添加水,那么也需要从水源处添加粒子)。在每个对流步骤中,根据网格速度场移动粒子,最后包含粒子的单元格都是水,其余单元格留空。在渲染过程中,在水和空气之间形成光滑的表面,也就是构建一层包裹在粒子周围的光滑的隐式表面:

其中是周围粒子位置的加权平均值,杠r是周围粒子半径的加权平均值,曲面定义为φ(x) = 0的所有向量x代表的点,也就是φ的零等高线或者水平集。


(利用粒子跟踪和渲染流体的二维模拟图)

8.2 水平集方法

在数学中,n个变量的实值函数f的水平集是形式为Lc(f)={(x1,…,xn)∣f(x1,…, xn)=c}的集合,也就是取定一个常数c然后取满足f(X)=c的集合就好了。

所以不用像8.1里一样围绕粒子构建隐式函数,而是直接使用网格。在网格单元的中心定义隐式曲面函数φ_(i,j,k),可以用插值来估计φ(x),表面就是满足φ(x)=0的点集,表面内部有φ(x)<0,外部有φ(x)>0

Signed distance(距离)

给定一个点集S,定义一个距离函数,可以获得任意一点到S中一点的最短距离:

如果S将空间划分为内部和外部的话,那么对应的符号情况如下:(外部为distance(x),内部为-distance(x)

综上,S就是φ(x)=0distance(x)=0的点集。

Signed distance properties(一些性质)

  1. 在表面内的某个点x处,设n为朝表面点集距离最近点的方向的单位向量,对于一定范围的正数ε,有φ(x+εn)=φ(x)+ε。因为如果沿n方向移动,在表面上的那个最近点是不变的。所以有下面这个式子:

  1. 增大或减小到表面的距离的最快方法就是沿着n,所以φ的梯度必须是n的方向(到最近点)。

  1. 将1、2中的两个式子(25)(26)结合得:

  1. 对于任意点,离表面的最近点为:(后面那坨其实就是求出来一个那个方向上的距离向量了)

Calculating the signed distance(计算)

给出一种基于几何的计算距离的算法,这个算法可以将有关表面的信息有效的传播到远处的网格点,而不用对表面上的每个点进行暴搜:

  1. 找到临近网格点的最近表面点,设定他们的距离,其他网格点的距离设为未知
  2. 以选定的顺序遍历周围的未知的那些网格点(i,j,k),如果新的点到表面的距离更近,那么将这个点周围的点标记为未知(这里的所有距离都按前文说的是带符号的距离)

此外,这个循环顺序是按照所有未知点从近到远的顺序,靠优先队列来实现即可。

(英语太差了,,这部分感觉没翻译明白,,,总之我觉得意思就是不用暴力,每次可以贪心的向那个距离下降最快的方向找到一个更优的参考点,然后把它周围的点加入未知集合中,在下一轮循环中继续从近到远循环未知点,,应该是这个意思吧)

Advecting Level Set

水-空气之间的表面已经解决了,还要考虑实心墙的边界条件,一种可行的方法是将φ从水-空气区域外推到固体中。

(下面咋外推的就没看懂了,,略过)

Current Fluid Simulation Loop

(没看懂0.0,,)

8.3 粒子水平集(PLS)

PLS用辅助标记粒子优化了上面8.2的欧拉水平集方法,在表面的正反面都添加辅助粒子,每个粒子除了自身位置外,还有一个φ值估计。在每个时间间隔粒子都会移动。此外还需要对那些溢出的粒子进行一点纠错操作:比如空气颗粒落入水中或者水颗粒落入空气中等等。最后,删除那些离零水平集太远的粒子,并在靠近边界的地方添加新粒子。


(利用粒子水平集方法跟踪流体的二维模拟图,红色粒子是空气边界,绿色粒子是水边界)


to be continued…

打赏
  • 版权声明: 本博客所有文章除特别声明外,著作权归作者所有。转载请注明出处!
  • Copyrights © 2018-2020 LeFlacon

奶茶一杯 快乐起飞

支付宝
微信