Real-Time Fluid Dynamics for Games 笔记
Stam J. Real-Time Fluid Dynamics for Games[J]. Journal of Graphics Tools, 2003, 6.
感谢朱老师的指导..本文对论文的实现代码进行分析,并试图补全论文中说得不够详细的部分,侧重于对solver的分析,opengl的部分比较简单,附件为已添加注释的代码。
我使用canvas api实现了一个基于web的版本,限于浏览器性能,在帧率上不如原文。可以在此链接中浏览效果:http://pigcage.tech/MyProjects/stablefluid/
一、系统的建立
1、网格系统
这里使用的是欧氏方法,代码中创建了一个N*N的网格(其中N=64,越大的N值意味着越高的性能开销)。在网格的四周有一层边界(0与N+1),亦即实际网格的大小为size=(N+2)*(N+2)。对于网格中的每一个单元格(cell),假定N*N网格的边长为1,那么每个单元格的边长为h=1/N。
基于此网格系统,程序使用一些全局变量来储存流体的基本属性,见附件中demo.c的注释如图:
时间步长为dt,一般来说,应满足dt < h * |u|,以保证模拟的精度。
每个单元格对应有密度(density)与速度(velocity)两个属性。其中速度为矢量,此处为二维模拟,故分为u、v两个方向表示,此处u指上下方,v指左右方。对于这些属性,分别使用形如dens、dens_prev形式的两个数组, 表示网格该属性的当前状态和前一状态。使用SWAP(x0,x)来交替这两个数组以达到状态的更新,其中x0为上一状态。这些数组均为一维数组,程序中使用形如dens[IX(i , j)]的形式,表示第j行i列的单元格属性。
2、准备知识
1.为了便于求解,使用到了Splitting方法,大体来说是这样的,感觉比较显然就不具体描述了吧..
2.N-S方程
3.顺便贴一下物质导数D
二、密度步骤
由于这里是处理烟雾,密度的高低等价于烟雾的浓度。使用密度数据即可绘制对应的烟雾视图。对密度的处理分为上图几个步骤,在每一帧中,依次添加力、处理扩散和移动。
1、力和源的添加
这一步比较简单,使用add_source函数。
void add_source ( int N, float * x, float * s, float dt ) { int i, size=(N+2)*(N+2); for ( i=0 ; i<size ; i++ ) x[i] += dt*s[i]; }
用户输入的源(source)或力(force)被储存在 s[ ] 中,那么把这些状态直接加到新一帧的对应属性 x[ ] 即可。至于乘以dt也就是因为使用了splitting的方法,并不是冲量..下同。
2、扩散(diffusion)
上面已经处理了外部输入源的情况。对于接下来的扩散步骤(烟雾的自发扩散)。根据N-S方程,我们知道扩散项为 k · ∇^2 u ,其中k为扩散系数,u为密度场。于是我们对某个局部的区域 x 而言,有 x - x0 = k · ∇^2 u 。这里是对x、y两个方向求偏导,也就相当于只考虑每个cell与其四领域的交互了。
设每一个帧的时间dt,扩散系数diff,对于上式我们有:
x[IX(i,j)] - x0[IX(i,j)] = a*(x[IX(i-1,j)]+x[IX(i+1,j)]+x[IX(i,j-1)]+x[IX(i,j+1)]-4*x[IX(i,j)]);
其中a=dt*diff*N*N,由于这里和后面projection用到了相同的方法,我想把对这条式子右值的解释留到最后。对上式进行化简整理,可以得到:
x[IX(i,j)] = (x0[IX(i,j)] + a*(x[IX(i-1,j)]+x[IX(i+1,j)]+x[IX(i,j-1)]+x[IX(i,j+1)]))/(1+4*a);
令t=a/(1+4*a),数组x0为已知的上一状态。我们可以把这个等式写成线性方程组,如下图:
(勘误:上图-1应为-t/a)
对于最左矩阵的第IX(i,j)行而言,对应的四个t的位置分别为第IX(i-1,j),IX(i,j-1),IX(i,j+1),IX(i+1,j)列,而-t/a的位置为第IX(i,j)列。这是一个稀疏矩阵,求解这个线性方程组用到了高斯-赛德尔迭代。关于这个方法的本身可参考https://www.cfd-online.com/Wiki/Gauss-Seidel_method。代码实现如下:
for ( k=0 ; k<20 ; k++ ) { for ( i=1 ; i<=N ; i++ ) { for ( j=1 ; j<=N ; j++ ) { x[IX(i,j)] = (x0[IX(i,j)] + a*(x[IX(i-1,j)]+x[IX(i+1,j)]+x[IX(i,j-1)]+x[IX(i,j+1)]))/(1+4*a); } } set_bnd ( N, b, x ); }
其中k为迭代的次数,此处设为20,较大的迭代次数可以保证好的精度,但性能的开销也更大。set_bnd()函数的作用为设置边界条件,后面再具体分析。
3、对流(advection)
由于速度场的存在,流体不仅有自发的扩散运动,也有因速度带来的流动。按照通常的思维角度,也就是N-S方程的对流项 - (u · ∇) ρ,我们可以根据每个单元格当前的u、v分量(下图a)计算出步长dt后的对应位置(下图b)。但这个位置并不是总在单元格的中心点上(下图b),这将导致计算出的结果难以使用。
因此,作者在此处采取了前向追溯的方式,使用半拉格朗日法(semi-lagrangian)通过当前单元格的速度来计算该处状态在前一dt时的位置(下图c)。这里讨论的是密度的对流,对于速度的对流也是一样的道理,只是速度是个矢量场,拆分成多个维度来计算。
实现的代码如下:
void advect ( int N, int b, float * d, float * d0, float * u, float * v, float dt ) { int i, j, i0, j0, i1, j1; float x, y, s0, t0, s1, t1, dt0; dt0 = dt*N; FOR_EACH_CELL //线性回溯上一时刻的位置(x,y) x = i-dt0*u[IX(i,j)]; y = j-dt0*v[IX(i,j)]; if (x<0.5f) x=0.5f; if (x>N+0.5f) x=N+0.5f; i0=(int)x; i1=i0+1; if (y<0.5f) y=0.5f; if (y>N+0.5f) y=N+0.5f; j0=(int)y; j1=j0+1; //四个邻居到点(x,y)线性插值,求得点(x,y)的值 s1 = x-i0; s0 = 1-s1; t1 = y-j0; t0 = 1-t1; d[IX(i,j)] = s0*(t0*d0[IX(i0,j0)]+t1*d0[IX(i0,j1)])+s1*(t0*d0[IX(i1,j0)]+t1*d0[IX(i1,j1)]); END_FOR set_bnd ( N, b, d ); }
这里主要解释一下d、d0的等量关系。
前一位置(x,y)可能位于网格的任意一点上,由于其不一定出于中心点,我们无法直接确定该点的密度。作者提出通过线性插值的方式,使用点(x,y)相邻四个中心点的密度值,来算出点(x,y)处的密度。具体如下图:
其中,红点处为点(x,y),相邻的四个网格内的黑点为中心点,左上角网格为假定点(x,y)经过时间步长dt后到达的当前位置(i,j),s0、s1、t0、t1为点(x,y)到四个相邻中心点的横纵距离。根据线性插值的规则不难得出(x,y)处的实际密度,而此处的密度,经过dt时长后到达点(i,j)处,故实际上就是当前中心点(i,j)的密度,如上述程序片段所示。对于后面速度步骤中,处理u、v分量,亦是相同的原理。
至此,对当前一帧密度状态的处理完成,如下:
void dens_step ( int N, float * x, float * x0, float * u, float * v, float diff, float dt ) { add_source ( N, x, x0, dt ); SWAP ( x0, x ); diffuse ( N, 0, x, x0, diff, dt ); SWAP ( x0, x ); advect ( N, 0, x, x0, u, v, dt ); }
三、速度步骤
对速度的处理与对密度的处理大致类似,由于速度有u、v两个方向的分量,因此对每一分量要经过同样的处理。不同的是,对速度的处理需要经历投影(project)步骤,得到一个不可压缩场(质量守恒场),实际上就是处理N-S方程中的压力项。
首先有以下概念:
至于为什么是这样来呢,以及project部分的详细数学过程,我想留到最后解释。这里Δx = Δy =1/N = h。根据霍奇分解定理(Helmholtz-Hodge Decomposition),速度场w可以被分解为不可压场u与梯度场∇p的和:
w = u + ∇p
u是不可压场,有∇ · u = 0,对于等式(1)两边分别进行 ∇ · 操作,得到:
∇ · w = Δ p
离散化:
等式左侧实际上就是散度场div[],我们可以通过已知的u、v分量求得,代码如下。
for ( i=1 ; i<=N ; i++ ) { for ( j=1 ; j<=N ; j++ ) { div[IX(i,j)] = -0.5*h*(u[IX(i+1,j)]-u[IX(i-1,j)]+v[IX(i,j+1)]-v[IX(i,j-1)]); p[IX(i,j)] = 0; } }
于是有:
类似前面在advection中所述,此处的式(1)同样可以写成线性方程组的形式,不再列出。因而这里同样可以用高斯-赛德尔迭代求解压力场 p[] ,如下:
for ( k=0 ; k<20 ; k++ ) { for ( i=1 ; i<=N ; i++ ) { for ( j=1 ; j<=N ; j++ ) { p[IX(i,j)] = (div[IX(i,j)]+p[IX(i-1,j)]+p[IX(i+1,j)]+p[IX(i,j-1)]+p[IX(i,j+1)])/4; } } set_bnd ( N, 0, p ); }
解出p以后,就可以算出梯度场 ∇p ,最后用速度场 w - ∇p = u,得到不可压场。
for ( i=1 ; i<=N ; i++ ) { for ( j=1 ; j<=N ; j++ ) { u[IX(i,j)] -= 0.5*(p[IX(i+1,j)]-p[IX(i-1,j)])/h; v[IX(i,j)] -= 0.5*(p[IX(i,j+1)]-p[IX(i,j-1)])/h; } }
至此,project步骤完成。
四、边界条件
所谓边界,就是网格最外层的一圈(0和N+1)。在计算网格的某一cell的具体状态时,我们通常要使用到与之相邻的几个cell的状态,那么对位于网格边缘的cell而言,它们在一些方向上自然是没有邻居的。为了解决这个问题,我们在网格的外面额外增加一层边界,并根据实际情况为其赋值。
先看实现:
void set_bnd ( int N, int b, float * x ) { int i; for ( i=1 ; i<=N ; i++ ) { x[IX(0 ,i)] = b==1 ? -x[IX(1,i)] : x[IX(1,i)]; x[IX(N+1,i)] = b==1 ? -x[IX(N,i)] : x[IX(N,i)]; x[IX(i,0 )] = b==2 ? -x[IX(i,1)] : x[IX(i,1)]; x[IX(i,N+1)] = b==2 ? -x[IX(i,N)] : x[IX(i,N)]; } x[IX(0 ,0 )] = 0.5f*(x[IX(1,0 )]+x[IX(0 ,1)]); x[IX(0 ,N+1)] = 0.5f*(x[IX(1,N+1)]+x[IX(0 ,N)]); x[IX(N+1,0 )] = 0.5f*(x[IX(N,0 )]+x[IX(N+1,1)]); x[IX(N+1,N+1)] = 0.5f*(x[IX(N,N+1)]+x[IX(N+1,N)]); }
参数b等于1、2时分别处理u、v两个方向的分量,也就是该速度方向到达对应的两个边界时发生碰撞,速度方向取反。当b为0时,处理的是密度值,密度不存在碰撞问题,故四个边界均保留原值。
对于网格的四个顶点而言,这里的状态为其相邻两个单元格的平均值。
五、完整过程
最后放一下实际执行模拟的代码。这里实际是先处理速度再处理密度,因为需要保证其不可压缩性。
while ( simulating ) { get_from_UI ( dens_prev, u_prev, v_prev ); vel_step ( N, u, v, u_prev, v_prev, visc, dt ); dens_step ( N, dens, dens_prev, u, v, diff, dt ); draw_dens ( N, dens ); }
六、散度、梯度及其它
现在我们来把散度、梯度以及上面提及的使用高斯-赛德尔迭代求解的等式具体解释一下。
首先说散度(divergence),某位置上的散度形容的是这个位置上向量场的发散程度。举个可能不太恰当的例子,你周围静止地站着一些人,这时候你这个位置的散度div=0;由于你今天喷了香水,周围的人纷纷向你凑过来闻,这时候你的散度就是负的了;然后你又放了一个屁,大家不堪其臭纷纷散开,于是你又成了一个正源,即div > 0。
回到流体中,对于一个可压缩的流体来说,某点周围的物质由于速度的存在,可能总体会朝这个点聚拢,这时候div<0,表现为该点附近的物质被压缩了;类似地,当div>0时,物质纷纷离开这个点,看起来就是这里的物质膨胀或者爆炸了。一般而言,我们研究的是不可压流体,或者换句话说,我们希望这个流体看起来是质量守恒的(mass conserving),物质不会凭空压缩或爆炸。也就是,我们需要保准散度为0,物质不能压缩也不能爆炸。更准确地说,流体的体积应该是不变的,而控制这一过程的正是压力p,也就是N-S方程的压力项。
现在来看散度算子 。div = ∇ · w ,其中∇ · 是散度算子,w是一个梯度场,运算的结果是纯量。我们讨论二维的情况。令w=( u , v ), ∇ · w = ∂u/∂x +∂v/∂y 。中心差分法,对二维平面上的速度场u( i , j )、v( i , j ) 有:
∂u/∂x = ( u(i+1,j) - u(i-1,j) )/2Δx
同理有:
∂v/∂y = ( v(i,j+1) - v(i,j-1) )/2Δy
于是我们得到散度计算公式:
∇ · w = ∂u/∂x +∂v/∂y = ( u(i+1,j) - u(i-1,j) )/2Δx + ( v(i,j+1) - v(i,j-1) )/2Δy
已知Δx=Δy=h=1/N,得到:
散度div = ( u(i+1,j) - u(i-1,j) + v(i,j+1) - v(i,j-1) ) / 2h
2. 接下来如何去除速度场的散度,使得 ∇ · w = 0,即实现质量守恒、不可压性呢?
根据亥姆赫兹 - 霍奇分解定理,任一速度场 w(向量场),可以唯一分解为一个无散度场u与一个梯度场的和。即:
w = u + ∇ p
所谓projection,也就是去掉∇ p,即取 w 在 u 上的投射。
那么对于上式,u是无散场,∇ · u = 0,对等式两边∇ · 有:
∇ · w = Δ p
对于此式,∇ · w = ∂u/∂x +∂v/∂y ,此处已知,也就是我们前面求得散度场div。
对于等式右边,Δ p = ∇ · ∇ p = ∂(∂p/∂x) / ∂x +∂(∂p/∂y) / ∂y ,不妨以x方向来举例,设x方向上有三个离散的点 x(i-1) , x(i) , x(i+1) ,此外x(i-1) , x(i)中点处有一点a, x(i) , x(i+1) 中点处有一点b。
中心差分法,对于a处我们有 f' (a) = [ f(i) - f(i-1) ]/h,同理对b有f' (b) = [ f(i+1) - f(i) ]/h
那么,对于i处有f'' (i) = [ f' (b) - f' (a) ] / h = [ f(i+1) - 2f(i) + f(i-1) ] / h^2。令 ∇ · w = Δ p = div,于是:
于是得到p的迭代公式:
这里注意,因为原程序中lin_solve()函数的写法关系,上式div对应的变量系数只能为1,因此程序在计算散度div时实际计算的时d = - h^2 · div,再代入上式计算。于是有:
至此可以用高斯迭代计算出压力场p。得到p后,求压力梯度场 ∇ p = ( ∂p/∂x , ∂p/∂y ),最后速度场减去此梯度场,得到不可压场。
七、附件
附件为添加了一定注释的完整实现代码demo.c与solver.c,由于为边分析边添加注释,难免错漏。本文亦出于整理思路的需要写成,错漏之处还望指正。
附件: