Chango的数学Shader世界（十一）流体模拟-Helmholtz分解，边界，场性质

最新推荐文章于 2021-02-01 15:39:25 发布

目竞

最新推荐文章于 2021-02-01 15:39:25 发布

阅读量2.1k

点赞数 5

分类专栏： Shader 计算机图形学数学物理

本文链接：https://blog.csdn.net/qq_41524721/article/details/97624904

版权

计算机图形学同时被 3 个专栏收录

42 篇文章 8 订阅

订阅专栏

Shader

23 篇文章 10 订阅

订阅专栏

数学

17 篇文章 2 订阅

订阅专栏

目的：

参考《GPU Gems》，尝试重现2D流体模拟。

本节解释文中对上上节Navier–Stokes equations的变形，及需要变形的原因，导出最终算法的结构。

参考：

《GPU Gems》

观察：

Why so serious

分析：

1.边界

归于实际，我们的这些场是要应用在一个有限的2D区域上的，不妨称此区域为D。一般来说，只有在其边界 $\partial D$ 上处处可微，才能应用这些定理和公式。不过我们由于是离散形式，又可以不考虑4个角，也能在边界上进行一般操作，问题不大，后面会谈到。

不过，我发现了一个从大学时期就没懂的问题，为啥区域D的边界总是写成 $\partial D$ ？这跟偏导有啥关系？

经过科学搜索，竟然发现了陶哲轩?的回答：

其实本质道理十分显然，下图既可以适用于2D，也可以想象成3D:

（注意：已更换符号）

当 $|\Delta \boldsymbol{r}|$ 趋向0，D区域增加的体积趋向于D区域的表面积（D区域边界的模），去除掉绝对值符号，把 $\partial D$ 当作D的边界也合情合理。

对于答主下面拿圆举例的公式，我觉得有失细致，乍看会误导人，以为圆面积对r导数一下就是周长，其他图像也一样。其误导之处在于这里的半径r与上面的r-neighbourhoods不是同一个东西。在圆这个例子中，如果要和上面的思想保持一致，那么就该这么写，这里我根据第一段的r（r-neighbourhoods）的思想，引出圆的位矢变化长度p。重设圆的半径为R。

那么，所谓的 $|S_{r}|$ 就是函数 $G(p)=\pi ((R+p)^{2}-R^{2})$ ，对p求导并令p=0，得 $G^{'}(p)|_{p=0}=2\pi R$ ，这样，这个形式和p就和上面的r-neighbourhoods对应了。

注意这里的函数G也不是唯一的。因为把一个图形扩出去，不一定要沿着法线方向，你想怎么来就怎么来，函数也千变万化，但是对微变化量求导，取0就为原图形周长是一定的。所以没把这个p规定死为位矢变化量或法线方向大小，就说是neighbourhoods。

另一种解释有关拓扑，咱也不懂，就不说了。

总的来说，用 $\partial$ 表示其边界没有什么硬道理，但也不是没道理。

2.Helmholtz-Hodge分解

既然1个矢量可以分成任意个矢量的和，那么矢量场也可以。

那从中就会有一些规律，比如Helmholtz分解，指出D上矢量场w能唯一分解为一个无散场和一个无旋场

无散场==无通量源场。

无旋场==保守场==某数值场的梯度场

所以在《GPU Gems》此文中，采用Helmholtz分解公式：

$\boldsymbol{w}=\boldsymbol{u}+\nabla p$

其中w为原矢量场，u为无散场（注意这里符号u只是占位符，不一定是之前的速度场u），p为数值场。

而且在 $\partial D$ 上有 $\boldsymbol{u}\cdot \boldsymbol{n}=0$ ，n为边界法线场。

这里我们可以粗略看出为何要用Helmholtz分解对我们的原方程一进行变形了。还记得上上节我们说，在这里需要的方程有2：

$\frac{\partial \boldsymbol{u}}{\partial t}=-(\boldsymbol{u}\cdot \nabla)\boldsymbol{u}-\frac{1}{\rho }\nabla P+k\nabla^{2}\boldsymbol{u}+\boldsymbol{F}$

$\nabla \cdot \boldsymbol{u}=0$

如果对方程一使用Helmholtz分解，不仅将方程二成立的条件隐含在其中，而且考虑到了边界。

并且，这个分解公式是具有唯一性和线性的。可惜文中并未提及或证明。

这两个性质是下面推论的前提。

现在，对分解公式进行进一步研究，得出几个结论，方便之后变形。

3.推论一

要将某个有散场w转为无散场u，可以减去某梯度场：

$\boldsymbol{u}=\boldsymbol{w}-\nabla p$

4.推论二

推论一中我们无法得知 $\nabla p$ 怎么求。在分解公式两边应用散度算子：

$\nabla \cdot \boldsymbol{w}=\nabla \cdot (\boldsymbol{u}+\nabla p)=\nabla \cdot \boldsymbol{u}+\nabla^{2}p$

又因为 $\nabla \cdot \boldsymbol{u}=0$ ，所以得到Poisson压力方程：

$\nabla^{2}p=\nabla\cdot \boldsymbol{w}$

这样，通过解此方程，可得到p，再代入推论一。

5.推论三

将分解公式定义转化成算子 $\mathbb{P}$ ，即 $\mathbb{P}\boldsymbol{w}=\boldsymbol{u}$ ,将一个有散场转为了无散场。

由于其线性和唯一性，可对分解公式本身两边应用 $\mathbb{P}$

$\mathbb{P}\boldsymbol{w}=\mathbb{P}\boldsymbol{u}+\mathbb{P}(\nabla p)$

显然又有 $\mathbb{P}\boldsymbol{u}=\boldsymbol{u}$

则 $\mathbb{P}(\nabla p)=0$

即对于梯度场，其无散部分只能是处处为0矢量的矢量场

6.推论四

对方程一两边都应用 $\mathbb{P}$

$\mathbb{P}\frac{\partial \boldsymbol{u}}{\partial t}=\mathbb{P}(-(\boldsymbol{u}\cdot \nabla)\boldsymbol{u}-\frac{1}{\rho }\nabla P+k\nabla^{2}\boldsymbol{u}+\boldsymbol{F})$

因为 $\boldsymbol{u}$ 是无散场，所以 $\frac{\partial \boldsymbol{u}}{\partial t}$ 也是无散场，其证明只需一步：

因为 $\frac{\partial u1}{\partial x}+\frac{\partial u2}{\partial y}=0$

所以 $\frac{\partial u1/\partial t}{\partial x}+\frac{\partial u2/\partial t}{\partial y}=0$ ,即 $\nabla \cdot \frac{\partial \boldsymbol{u}}{\partial t}=0$

所以左边去除符号，再结合推论三，右边去除压力场梯度项：

$\frac{\partial \boldsymbol{u}}{\partial t}=\mathbb{P}(-(\boldsymbol{u}\cdot \nabla)\boldsymbol{u}+k\nabla^{2}\boldsymbol{u}+\boldsymbol{F})$

7.整合

现在，我们结合边界条件，方程一，方程二，得出新的方程：

$\frac{\partial \boldsymbol{u}}{\partial t}=\mathbb{P}(-(\boldsymbol{u}\cdot \nabla)\boldsymbol{u}+k\nabla^{2}\boldsymbol{u}+\boldsymbol{F})$

回忆起上上节的解释平流时，实现起来其实操作的是染料场d。若再次将其看成x,y分量的染料浓度矢量场d=(d1,d2)。

则有：

$\frac{\partial \boldsymbol{d}}{\partial t}=\mathbb{P}(-(\boldsymbol{u}\cdot \nabla)\boldsymbol{d}+k\nabla^{2}\boldsymbol{d}+\boldsymbol{F})$

注意其中压力项被隐含到了 $\mathbb{P}$ 中，也就是推论一和推论二导出的将有散场转化为无散场的过程中。

这就是我们算法的理论基础。

我们现在假设速度场不变。有了初始的染料浓度场（单张图RenderTarget,每点的r,g颜色分量用于存储d1,d2）后，每帧更新：加上染料变化场 $\frac{\partial \boldsymbol{d}}{\partial t}$ 的离散形式，也就是另一张RT，就能得到新的染料场。

对于这张染料变化场RT，希望利用4个函数操作，每个函数的输入和输出都是一个场（一张RT图）。

平流A,扩散D，外力F,投射P：

RT functionA(RT d);

RT functionD(RT Ad);

RT functionF(RT DAd);

RT functionP(RT FDAd);

其中functionP涉及到解Poission压力方程，比较复杂，下次再谈。这些函数的作用,换成算子就是：

$\mathbb{A}\boldsymbol{d}=-(\boldsymbol{u}\cdot \nabla)\boldsymbol{d}$

$\mathbb{D}\boldsymbol{d}=\boldsymbol{d}+k\nabla^{2}\boldsymbol{d}$

$\mathbb{F}\boldsymbol{d}=\boldsymbol{d}+\boldsymbol{F}$

$\mathbb{P}\boldsymbol{d}=\boldsymbol{d}-\nabla (SolvePoission(\boldsymbol{d}))$

你问我为啥要换成算子写一遍?可能这样论文看起来帅吧。

然后，整个一帧的染料变化场，可写成对旧染料场 $\boldsymbol{d}$ 施加一个算子 $\mathbb{S}$ 。

$\frac{\partial \boldsymbol{d}}{\partial t}=\mathbb{S}\boldsymbol{d}=\mathbb{P}\mathbb{F}\mathbb{D}\mathbb{A}\boldsymbol{d}$

算子从右到左进行运算。这里省略了时间t，但t参与运算。

结语：

本节解释《GPU Gems》文中对上上节Navier–Stokes equations的变形，及需要变形的原因。将其算法最终导出4个函数/算子的结构，并引出了解Poission压力方程的问题。下节结合实际情况对各部分算子离散化，代码化。

目竞

关注

5
点赞
踩
11

收藏

觉得还不错? 一键收藏
1
评论
Chango的数学Shader世界（十一）流体模拟-Helmholtz分解，边界，场性质

目的：参考《GPU Gems》，尝试重现2D流体模拟。本节解释文中对上上节Navier–Stokes equations的变形，及需要变形的原因，导出最终算法的结构。参考：《GPU Gems》观察：Why so serious分析：1.边界归于实际，我们的这些场是要应用在一个有限的2D区域上的，不妨称此区域为D。一般来说，只有在其边界上处处可微...
复制链接

扫一扫