Jos Stam Stable Fluids中一些问题解读

Method of Characterstics特征方法

Stable Fluids文章中使用特征方法解决Advection(对流)。

维基百科对特征方法的解释:

https://en.wikipedia.org/wiki/Method_of_characteristics

就是找一条让PDE成为ODE的线,沿着这条线解ODE,这样就搞定了原来的PDE问题。

看论文中的解释:

特征方法用来解这样的方程:

这个\(\frac{\partial a(\textbf{x},t)}{\partial t}+\textbf{v}(\textbf{x})\cdot\nabla a(\textbf{x},t)\)其实就是随体导数。

然后论文定义\(\textbf{p}(\textbf{x}_0,t)\)为速度场\textbf{v}的characterstics(特征),满足

这个\(\textbf{p}(\textbf{x}_0,t)\)的意思就是在时刻0时位于\(\textbf{x}_0\)的这个流体点在\(t\)跑到了什么地方。

定义

为场a的“随体场”,也就是采取拉格朗日观点,定义它为0时刻在\(\textbf{x}_0\)的流体点在t时刻的值。

上面说了,我们要解的方程相当于随体导数,那么\(\bar{a}\)就是“随体场”,所以求解的直接就是它的时间导数

具体来说,在做Advection的时候,沿着流线(characterstic line)“上溯”一步即可。

边界条件

Jos Stam代码中的边界条件如下:

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)]);
}

Stable Fluids采用的边界条件是Dirichlet边界条件,也就是:在边界处,速度场的法向分量为零。即,没有任何通量“穿过”边界。

在文章中,速度场的离散点分布如下:

注意这里的坐标是1~N.想要实现“法向为零”的边界条件怎么办呢?我们取边界上的一个点,算该点处的速度值时,会取“外面”的速度点:

就像这样,红点的速度值是格子左上角黑点和格子外面蓝点的平均值。

那么,如果想让红点处速度的y分量为零,只需要把蓝点速度的y分量置为黑点速度的y分量的相反数即可。

蓝点的x分量呢?红点速度的x分量应当和黑点的x分量相同,那么蓝点的x分量也应该等于黑点的x分量。

最后,set_bnd函数的结尾四句是设置四个角的边界条件:

红点的速度是两个蓝点速度的平均值。原因是,例如我们考察y分量,令红点为y00:

y01=-y11

y10=y11

考虑到角上的速度应当为0,也就是说(y00+y01+y10+y11)/4=0,代入得到y00+y11=0,x的情况类似,也就是说红点的速度应该和黑点的速度恰好相反。

但是Jos Stam代码中,这个点的速度应该是0,原因有待实验。

Poisson Solver

在Stable Fluids里面,求解泊松方程的代码长这样:

void lin_solve ( int N, int b, float * x, float * x0, float a, float c )
{
	int i, j, k;

	for ( k=0 ; k<20 ; k++ ) {
		FOR_EACH_CELL
			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)]))/c;
		END_FOR
		set_bnd ( N, b, x );
	}
}

这是啥玩意呢?查了一圈发现,这玩意就是线性方程组的Jacobi Solver.看一眼维基百科:

https://en.wikipedia.org/wiki/Jacobi_method

注意,每次是“基于”对角元去更新解向量,这可以解释为什么更新的是x[i][j].

Jacobi Solver的条件是:对角元最大。这个怎么满足呢?把泊松方程离散化以后,就变成了类似拉普拉斯算子的玩意:

比方说一个大致长这样的算子,中间被加了好多次,可不就是最大了嘛。每个格子一个方程,自然对角线上的元素最大。

发表评论

电子邮件地址不会被公开。 必填项已用*标注