文章略读:A self-adaptive oriented particles Level-Set method for tracking interfaces

在初始时间[latex]t_0[/latex],在界面上放[latex]N_p[/latex]个粒子,用[latex]\mathbf{x}_p[/latex]表示粒子的位置。

显然,粒子跟着流场走,也就是:

由于是粒子,所以就是直接微分,没有什么物质导数之类的。

我们用有符号距离函数[latex]\phi[/latex] 做水平集的函数(level function)。那么根据这个[latex]\phi[/latex]就能求法向:

上面求出了粒子的时间演化方程,利用[latex]\phi[/latex]的性质可以把法向场的物质导数展开:

然后文章又对这个法向场的物质导数做了一下变换:

这个记号比较难懂,单独注明一下。首先是用了偏导的角标简写。在角标中逗号后面的做偏导。比如说[latex]\phi_{,i}[/latex]就是[latex]\phi[/latex]对第i维坐标(x,y,z之一)做偏导,如果 [latex]\phi_{,it}[/latex]就是对第i维和时间求二阶导, [latex]\phi_{j,i}[/latex]就是[latex]\frac{\partial\phi_j}{\partial\phi_i}[/latex]的意思。另外一个就是爱因斯坦求和,所有重复的角标都表示求和。比如说[latex]\phi_{,k}\phi_{,kt}[/latex]中这个重复的k就暗示一个求和,也就是说它实际是[latex]\Sigma_k\phi_{,k}\phi_{,kt}[/latex].

这个推导左端是i,因为法向场的物质导数是一个矢量,这是其中的第i维。前几个等号都可以通过微分运算法则推出来。倒数第二个等号,是专门凑了一下[latex]\phi_{,t}+\mathbf{u}_j\phi_{,j}[/latex]这个形式,把它用全微分展开即得正确性。为什么要凑这个形式,是因为文章中的水平集函数[latex]\phi[/latex]满足守恒方程:

这样凑出来的项都是零了,就只剩下两项。然后根据前面说的, [latex]\phi[/latex] 和法向的关系,很容易验证最后的结果,也就是:

[latex]-u_{j,i}n_j+u_{j,k}n_jn_kn_i[/latex]

注意这只是一维,三维加起来的结果是:

值得高兴的是,在这个式子中 [latex]\phi[/latex] 消失了。也就是说,法向的演化仅仅和速度场有关。文章采用一个三阶Runge-Kutta算法去演化法向场。

文章中采用的是有向粒子法,也就是说粒子上面携带位置信息和当地法向信息。

每一步法向都会重新normalize,放缩到长度为1.

虽然粒子提供了表面信息,但并不知道粒子之间的相邻属性,所以不能直接用于生成表面。为了生成表面,还是需要搞一个水平集函数[latex]\phi[/latex].

方法是:遍历所有粒子,确定interface cells(至少包含一个粒子的格子)。在遍历过程中,计算离每个格点最近的粒子距离。然后用这个距离和该粒子的法向做点积,即得这些点处水平集函数(有向距离)。

然后水平集函数满足Eikonan方程:

就可以把它延拓到整个计算域上。

随机采样算法

文章在此处先讨论了“表面上的两个点”之间的信息,或者说从一个点去找“下一个点”。

如图。P是一个点,然后分析表面上的另一个点的信息。这个点就是右半边图中,P右下角处的那个点。它的位置就是[latex]\mathbf{x}_P+\beta\mathbf{t}+\alpha\mathbf{n}[/latex].其中[latex]\mathbf{n}[/latex]很好理解,就是P点处的法向。至于[latex]\mathbf{t}[/latex],则是沿表面的一个切向量。当然这样的切向量有许多,文章应该是先假设取了一个平面(图中淡黄色平面,这个平面包含[latex]\mathbf{n}[/latex]),然后全部在这个平面上讨论。

然后文章基于这个关系推导了一下。

第一步,是基于水平集函数[latex]\phi[/latex]进行推导。

把[latex]\phi[/latex]关于P点做二阶泰勒展开:

文章在这里推了一下中间那个二阶导项的具体表示:

注意这里用到了[latex]\phi[/latex]和法向的关系,尽量使用[latex]\mathbf{n}[/latex]表示。

然后把 [latex]\mathbf{x}_P+\beta\mathbf{t}+\alpha\mathbf{n}[/latex] 这个点的位置信息代入,一番推导后可得:

此处注意,新的点也应该在表面上,其[latex]\phi=0[/latex],所以这个式子等于零。

第二步,基于几何进行推导。

还是看上面那张图的右半部分。两点之间的弧长为

这里边[latex]\tilde{\kappa}[/latex]就是和[latex]\mathbf{t}[/latex]相切这条曲线的曲率了。

进而,使用这个式子,可以把两点之间的弦长表示为

可以看到,[latex]\beta[/latex]是直角三角形的直角边,可用三角函数求其长度:

这里经过了一步三角恒等式代换和一步泰勒近似。我们可以发现这个式子给出了[latex]\beta[/latex]和s的一个关系,也就是说沿着曲面走s,能让[latex]\beta[/latex]变化多大:

然后[latex]\alpha[/latex]也可以用三角函数推导出来:

这里可以发现,[latex]\alpha[/latex]是 [latex]\alpha^2[/latex]这个量级的。把这个关系代入第一步推出来的[latex]\phi[/latex]。取一个比较小的[latex]\beta[/latex]:

然后就可以在新的[latex]\phi[/latex]表达式中消掉[latex]\beta[/latex]的三阶项:

中间这个[latex]\phi_{,ij}\mathbf{t}_i\mathbf{t}_j[/latex]也可以用法向表示。文章里说是根据第一步里面[latex]\phi_{,ij}[/latex]得来,不过我感觉好像实际的推导应该是根据[latex]\phi[/latex]和法向的关系:

首先[latex]\mathbf{n}=\nabla\phi/|\nabla\phi|[/latex],取一维也就是说[latex]n_i |\nabla\phi|=\phi_{,i}[/latex],两边再对第j维求导得[latex]\phi_{,ij}=n_{i,j}|\nabla\phi|+n_i\partial|\nabla\phi|/\partial x_j[/latex].而水平集函数的性质,要求[latex] |\nabla\phi| =1[/latex],偏导总是零,所以第二项就没了,只剩 [latex]\phi_{,ij}=n_{i,j}|\nabla\phi| [/latex].

这只是一个式子,乘上[latex]t_i t_j[/latex],再求和(别忘了上面那些都是自带爱因斯坦求和的),就得到了:

再代回去,结果就是

这样就没有[latex]\phi[/latex]了,直接用P点处的几何信息,就能用[latex]\beta[/latex]得到[latex]\alpha[/latex].

接下来,使用Frenet第一公式可以得到切向和法向的关系(公式就长这样):

两端点积[latex]\mathbf{n}[/latex]就可以得到曲率:

第二个等号用到了Frenet第二公式。第三个等号,把ds写开推一推可以得到。

再回忆一下之前说的,新的点的位置:

[latex]\mathbf{x}_P+\beta\mathbf{t}+\alpha\mathbf{n}[/latex].

它的第i维就是:

[latex]x_i^{new}=x_i+\beta t_i+\alpha n_i[/latex].

我们刚才推出来了[latex]\alpha[/latex]的公式,略去小量,代入得

[latex]x_i^{new}=x_i+\beta t_i-\frac{1}{2}n_{k,l}t_k t_l\beta^2 n_i[/latex].

最后一项可以用曲率表示:

[latex]x_i^{new}=x_i+\beta t_i+\frac{1}{2}\tilde{\kappa}n_i \beta^2 [/latex].

注意此处文章里的公式最后一项没有[latex]\beta^2[/latex],我认为是作者搞错了,否则说不通。

相应的,新的点处的法向是:

曲率张量是:

还记得我们之前说的,文章假定了切向[latex]\mathbf{t}[/latex],但之前没说明怎么选取这个切向。选取的方法包括:

1、计算shape operator的特征向量,用来确定主曲率(principal curvatures)方向,并把切线设成沿这些方向。

2、简单地选定两个垂直的方向。

文章里用的是第二种。对于每一个初始种子,在其四周用上面讨论的方法,取四个采样点,也就是两个垂直的方向,及其相反方向。

发表回复

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