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

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

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

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

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

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

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

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

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

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

\(-u_{j,i}n_j+u_{j,k}n_jn_kn_i\)

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

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

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

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

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

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

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

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

随机采样算法

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

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

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

第一步,是基于水平集函数\(\phi\)进行推导。

把\(\phi\)关于P点做二阶泰勒展开:

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

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

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

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

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

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

这里边\(\tilde{\kappa}\)就是和\(\mathbf{t}\)相切这条曲线的曲率了。

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

可以看到,\(\beta\)是直角三角形的直角边,可用三角函数求其长度:

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

然后\(\alpha\)也可以用三角函数推导出来:

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

然后就可以在新的\(\phi\)表达式中消掉\(\beta\)的三阶项:

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

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

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

再代回去,结果就是

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

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

两端点积\(\mathbf{n}\)就可以得到曲率:

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

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

\(\mathbf{x}_P+\beta\mathbf{t}+\alpha\mathbf{n}\).

它的第i维就是:

\(x_i^{new}=x_i+\beta t_i+\alpha n_i\).

我们刚才推出来了\(\alpha\)的公式,略去小量,代入得

\(x_i^{new}=x_i+\beta t_i-\frac{1}{2}n_{k,l}t_k t_l\beta^2 n_i\).

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

\(x_i^{new}=x_i+\beta t_i+\frac{1}{2}\tilde{\kappa}n_i \beta^2 \).

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

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

曲率张量是:

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

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

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

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

发表评论

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