Notes On Stable Fluids

Posted on 2009-12-17 10:11 安泰 阅读(...) 评论(...) 编辑 收藏

 

Stable Fluids and  Real-Time Fluid Dynamics for Games (Jos Stam 1999)

这两篇文章是我最早接触的关于图形学方面的文章,前前后后看了很多遍,每一次都有新的收获。在进入正题之前,首先说下这个阅读理解能力其实是我们很多人的Limitation, 英语和数学作为两门语言来讲是研究的硬功夫,而且you can never be TOO good at these stuffs.

0、Prework

百度一下什么是流体、什么是图形学、Navier, Stokes何许人也:)

1、Introduction

首先Stable Fluids的基础是Navier-Stokes方程组(下称NS方程组)。这个方程组在流体动力系统的计算中举足轻重,可以说是作为流体模拟的首选。而Jos Stam的Stable Fluids一文在近些年的研究中具有里程碑式意义,很多人都是以此为基础进行其研究。

 

我们先看看NS方程组的形式。事实上,NS的表达有多种形式,这里的是出现在Jos文章中的形式。指出一点,这里的流体是不可压缩(incompressible)流体,也就是密度、温度不变。

这里其实一开始我郁闷了好久,因为对数学表达比较生疏。不妨看看展开形式。

回来,

 

(1)式中两项的点乘(dot product)表达的意思是速度场u的散度场为0。具体到场内某一点,各个方向汇聚来的和发散出的量之和为零。(1)式反应了质量守恒(如果你想研究原子弹爆炸,或者黑洞什么的,就不用质量守恒了:)BTW)。

(2)式中各项的物理意义P是压力场,希腊字母p表示密度,希腊字母v表示粘滞系数(viscosity),最后一项f是外力(external force)。 对(2)式的把握是个关键,现在我就对等式右边逐项讨论:第一项叫对流项(advection),速度沿着自身的方向运动。可以看出,速度快的地方就会越来越快,漩涡阿、龙卷风啊的貌似就是这个原理。因为u出现了两次,所以该项非线性。第二项的意思是,速度沿着一个压力梯度场运动,速度在高压力下会溢出去(凭空减少/增多)。这当然是不符合质量守恒的,Jos在自己的ppt中忽略了这个,计算中压根也没算,后文会讲为什么。第三项告诉我们,速度在自身梯度下扩散。扩散越强烈,总的场就越平滑。第三项是扩散项。最后一项,是外力项。大家看哈,(2)式左边是加速度(速度关于时间求导),其实这个式子是符合基本形式f=ma的。(其中推导过程这里就不说了,好事者可以看看这个大哥http://www.cnblogs.com/zxx_1987的推导,我拜服了) 

 

现在看上面这个式子,形式上和(2)式有很多相似的地方。事实上,(2)式可以用来建模很多东西。在Real-Time Fluid Dynamics for games里Jos就用了密度来产生效果。

2、Discretization

仅仅有了公式,对于计算机模拟的操作是不够的。一系列的离散化手段和编程trick需要掌握,另外就是关于simulation有很多注意事项,我个人的感悟而言:

 

i)对导数的离散化理解。二维三维中(对于三维以上的是一个challenging不过三维以上的模型很大程度上已近beyond我们的思考模式),对于各种差分(forward,back ,central)的实施,我们一般会用网格铺分,一般而言,我们是均匀的铺分成小格子,体素,欧拉栅格,what ever it’s called, 也就是说dx=dy=dz.这样在梯度、散度、拉普拉斯算子什么的我们分母就统一了,减少麻烦。

ii)边界条件(Boundary Condition, 下称BC)。Jos 此文提到了两个,Perodic BC和Fixed BC. 简而言之,Perodic BC就是说空间中重复着无限个我们的主模拟空间central box, 从central box 左边出去的粒子,会从central box 右边进入。Fixed BC就是说,边界值取一个固定的值。为什么会有边界问题啊?因为我们处理导数的离散问题时候,比如第(3,8)号点用到用到了(3,9)(3,7)(2,8)(2,9)四个点,但是边界(0,1)就没有四邻域了,(0,0)更没有了,我们需要给他们个交代。(以上描述的仅仅是最最简单的边界情况)

iii)合适的dt。很多前人的算法比不上Jos Stam就是因为随着dt取值出现了很多不希望出现的局面,比如不stable,或者数据跳跃。这里指出一个经常用到的条件,Courant-Friedrichs-Lewy Condition,下称CFL。CFL的说白了就是,dt不能太大,dt*u<dx,只有这样才能保证信息能够随着栅格传播。实施中dt的取值有多种方法。

iv)数值方法。因为离散计算不得不考虑精度,不同阶的精度在时间累积和错误累积下的效果是迥异的。Upwind Differencing, ENO, WENO, TVD Runge-Kutta。又是龙格库塔、切比雪夫、雅克比、欧拉这伙牛人:)

3、Implementation

 

根据Helmholtz-Hodge Decomposition,任意矢量场可以被唯一的分解为一个散度为0的场和一个标量场(scalar field)。也就是说任意的一个场,是一个质量守恒场和一个梯度场的叠加。我们要求解的正是质量守恒的场,所以我们要从总场中减掉那个梯度场。我们要进行投射操作P,把总场w投射到散度为零的场u. 上面是我的一页幻灯片。Jos的求解速度场就是按照(5)式的思路,一步步进行(i)加力、(ii)扩散、(iii)对流、(iv)投射。

这里谈下我的理解,为什么Jos忽略了(2)式中压力场项。因为压力梯度场在投影操作中会被消除,所以与其先前算出来再被消除,还不如直接不用算他。

i) 加力/源:对每个单元 x+=f*dt

ii)扩散。也就是 拉普拉斯算子(u), 每个单元和周围的邻域交换。

x0: previous state, x: current state, a: diffuse rate

num : number of its neighbors

sum_p: the sum of its neighbor (previous moment)

sum_c: the sum of its neighbor (current moment)

一个简单的实施

For Each Cell:

x = x0 +a*(sum_p – num* x0 )*dt

(第二项可能为负,产生摇摆和发散,所以不这么做)

Jos’ implement: trace BACKWARD in time

For every cell:

x0 = x - a*(sum_c – num*x)*dt

x0= x – a*sum_c*dt+a*num*x*dt

(1+a*num*dt)x = x0+a*sum_c*dt

x = (x0+a*sum_c*dt)/(1+a*num*dt)

上式的分母衡大于1,所以呢比较稳定

Jos 认为对这个的求解没必要进行求逆(matrix inversion routine)(He thinks this is overkill for this problem because the matrix is very SPARSE). 他用了一个比较实用的方法就是Gauss-Seidel relaxation, 迭带若干次,松散化一下。我的体验就是,其实没必要带很多次,稍微几次效果可以很不错了:)

(iii)对流

这里要求速度沿自身方向运动了,前面说了该项已经不是线性的了,所以不能用Gauss-Seidel relaxation来搞。这会用上了半拉个朗日方法,貌似现在纯粹的欧拉方法已经木有了。我们把每个网格的中心建模为一个粒子,这里还是用了trace BACKWARD in time的思想,沿着速度场回溯,看dt时刻之前位于哪个位置的粒子在会运动到当前的位置(s=so+v*dt, so=s-v*dt),用线性插值的方法求出dt时间前那个粒子,再把该粒子的属性值赋给当前粒子。线性插值么,就用最简单的类似求质心贡献的方法搞搞。

(iv)投射

首先说明白了,标量场比如密度,是没有必要投射的,投直接投木有了。

对于矢量场,比如速度场,

Step1) call div for every cell

div = (ui+1 – ui-1+ vi+1 – vi-1 + wi+1 – wi-1)/2dx

Step2) Gauss-Seidel relaxation

Find p[ ], which is the gradient filed

Step3) subtract the gradient filed

u -=(p[i+1,,]-p[i-1,,])/2dx

v -=(p[,j+1,]-p[,j-1,])/2dx

w -=(p[,,k+1]-p[,,k-1])/2dx

 

4、summary

本来以为一会就写完了,没想写了这么久。还是这么觉得,每一次看大牛的文章都有新的收获和感悟,不断地查询和搜索相关的知识又学到了新的东西,好比滚雪球的过程。要从根本的把握,还是数学层面的东西。实践时,离散化的思想及其重要。至于在什么环境、用什么编程语言不重要了。这个问题说白了就是个计算物理的过程,数学就是个工具,而所谓图形学就是个壳。