SPH(光滑粒子流体动力学)流体模拟实现四:各向异性(Anisotropic)表面光滑(1)

流体模拟四:

各向异性(Anisotropic)表面光滑(1)

 

表面的表示与定义

我们是以隐式表面来表示流体表面。定义一个足够大的参考域\large \overline{D},它完全包含粒子表示的流体所在的空间D,即\large D\subseteq\overline{D}。流体表面\large S=\partial D隐含的定义为空间函数\large \phi(x):\mathbb{R}^{d}\rightarrow \mathbb{R}的一个等值表面,即:

另外,流体的内外区域可表示为:

下图所示为隐式曲面表示表面示意图:

其中,绿色表示曲线内\large (\phi >0),蓝色表示曲线外\large (\phi <0),颜色的深浅代表场强的大小, 红色的曲线是值为0的等值面\large (\phi =0)。上述表示很容易离散化,只需要将\large \overline{D}等分成细小的网格\large x(i,j,k)\subseteq \overline{D},其中(i,j,k)是网格的索引。 离散化后的表示为:

在以前的各项同性方法中,我们采用如下方法定义等值面:

其中,C是一个正的常数,核函数\large W(x-x_{j},h_{j})可以选取任意适合的平滑函数。

 

基于各向异性核函数的表面重构

通常情况下,颜色场定义为:

其中,W 是各向同性核函数,其形式为:

\large \sigma为比例系数,d是维数,r是距离向量,P为有限支持域的衰减样条函数,因为我们的颜色场为密度场,所以该各向同性函数就是光滑核Poly6函数。但是,所得的表面通常具有隆起,这有两个原因。首先,粒子的不规则放置使得很难表现出绝对平坦的表面。尽管不规则采样是任何Lagrangian方案的基本特征,但边界粒子位置的这种不规则可能会使表面显得粗糙。第二,平滑内核的球形形状不适合描述表面附近的密度分布。即,为了正确地对表面几何形状建模,要使近表面粒子的半径沿不同方向以不同的速率降低。

为了处理第一个不规则放置的问题,这里引入了拉普拉斯平滑的3D变体,具有使点云降噪的效果。计算更新的内核中心\large \overline{x_{i}} :

其中w是合适的有限支持权重函数,而λ是0 <λ<1的常数。我们通常使用0.9到1之间的λ以达到最大化平滑效果。请注意,此平滑过程仅用于曲面重建,并且平均位置不会带回到模拟中。通常,拉普拉斯平滑会导致体积缩小,并且该方法还通过将边界粒子的内核向内移动来略微缩小流体体积。

为了处理第二个粒子球形形状问题,我们在颜色场定义函数中引入了矩阵,采用各项异性核函数重构流体表面,其表达式为:

其中,矩阵G将对r进行了旋转和缩放,  \large r=x-x_{i}为网格点到粒子的距离,x视为任意位置。因此通过G矩阵变换可以使我们的流体粒子变换为椭球体,在表面的粒子变形的更加扁平。

该方法构成了新的颜色场:

该方法为每个粒子确定一个各向异性矩阵G,以便更准确地描述粒子周围的密度分布。例如,在流体体积内的粒子附近,密度可能在所有方向上都是恒定的,这使对应的G为单位矩阵的标量倍数,以保持平滑核W各向同性。另一方面,在靠近平坦表面的粒子周围,粒子密度沿法向轴的衰减快于沿切向轴的衰减。然后,G应该沿切线轴拉伸W,并沿法线轴收缩W。在尖锐特征处,密度将在几个方向上急剧衰减,并且G应该缩小W以捕获尖锐特征。关于流体区域附近的各向同性和各向异性内核之间的比较如下图所示:

其中(a)为各项同性方法的粒子,(b)为各项异性的粒子。

为了确定G,我们将在Koren和Carmel中提出的加权分量主成分分析(WPCA)应用于邻近粒子位置。常规PCA的一个缺点是它对异常值的敏感性,当样品数量少且样品位置嘈杂时,通常会产生不准确的信息,这通常发生在基于颗粒的流体中。相反,WPCA通过为数据点分配适当的权重来实现对异常值和嘈杂数据的显着鲁棒性。具体而言,WPCA首先计算数据点的加权平均值。接下来,WPCA构造一个经验平均值为零的加权协方差矩阵C,并对C进行特征分解。所得特征向量给出主轴,特征值表示沿相应特征值的点的方差。然后,我们构造一个各向异性矩阵G,以将平滑核W与WPCA的输出进行匹配。

在我们的方法中,粒子i的加权均值\large x^{w}_{i}和协方差矩阵\large C_{i}表示为:

其中是的粒子i到粒子j的距离ri的各向同性(注意)加权函数:

的有限支持域下,计算仅限于半径\large r_{i}内的邻域粒子。 在我们的示例中,我们选择\large r_{i}\large 2h_{i}以便包含足够的邻域粒子并获得合理的各向异性信息。这里选择\large 2h_{i}是因为,我们的WPCA需要样本,我们如果选择\large r_{i}\large h_{i}作为样本会因为包含太少的周围粒子信息造成较差的效果。

对于每个粒子,相关联的C的奇异值分解(SVD)给出了拉伸或压缩方向,以根据特征向量和特征值使平滑内核W变形。 SVD生成:

论文中采用了奇异值分解来获取特征向量和特征值,由于我们的矩阵都是3*3的所以其实是和特征分解达到同样的效果的。由于我们的R是正交阵,所以这里的\large R\Sigma R^{T}等于\large R\Sigma R^{-1}

我们假设:

\large R=\begin{bmatrix} \frac{\sqrt{2}}{2} &-\frac{\sqrt{2}}{2} \\ \frac{\sqrt{2}}{2}& \frac{\sqrt{2}}{2} \end{bmatrix}\large \sum =\begin{bmatrix} \frac{1}{2} & 0\\ 0& 1 \end{bmatrix}\large R^{-1}=\begin{bmatrix} \frac{\sqrt{2}}{2} &\frac{\sqrt{2}}{2} \\ -\frac{\sqrt{2}}{2}& \frac{\sqrt{2}}{2} \end{bmatrix}\large \overrightarrow{v}=(1,0)

 

开始时,\large \overrightarrow{v}=(1,0)在图中是:

\large R^{-1}*\overrightarrow{v}在图中:

\large \sum *R^{-1}*\overrightarrow{v}为:

\large R*\sum *R^{-1}*\overrightarrow{v}为:

根据上面的图,根据我另一篇文章“https://blog.csdn.net/qq_39300235/article/details/103181337”在基变换中讲的,该变换恰好为:

1.将所有向量变换为R基中的向量。

2.将所有的向量按照特征值进行缩放。

3.将向量转换为原来基中的向量。

以此完成我们流体粒子的变形。然后为了处理奇异矩阵并防止极端变形,我们检查是否\large \sigma _{1}\geq k_{r}\sigma _{d},其中\large k_{r}>1。当一个主轴上的最大方差大于另一轴线上的最小方差\large k_{r}倍时,我们认为此为极端变形情况。 在这种情况下,我们修改矩阵C,以使任意两个特征值之间的比率小于\large k_{r}。 另外,当附近的粒子数量较少时,我们通过设置\large G=k_{n}I(I为单位阵)将其重置为球形,以防止几乎孤立的粒子出现较差的粒子变形。 另外,我们将C与比例因子\large k_{s}相乘,以使相关流体内部粒子的∥ksC∥≈1(即ksC的矩阵范数约等于1),“矩阵范数反映了线性映射把一个向量映射为另一个向量,向量的“长度”缩放的比例。”以便使具有完整邻域的粒子的体积保持恒定。 前述过程公式如下,以获得修正的协方差矩阵\large \widetilde{C}

其中􏰆\large \sigma _{k}= max(\large \sigma _{k},\large \sigma _{1}/ \large k_{r}),N是相邻粒子的数量,Nε是阈值常数。 在我们的示例中,我们使用kr = 4,ks = 1400,kn = 0.5,Nε= 25。

为了使粒子i的核W根据􏰆特征值变形,\large G_{i }必须是􏰆的倒数,并按1 / \large h_{i}缩放以反映粒子i的原始半径。 然后我们的方法将Gi生成为以下形式的对称矩阵:

到这里,我们各项异性流体平滑中,最重要的一步就完成了,我们接下来用获取到的G矩阵代入颜色场函数(流体模拟中为密度场):

即在每次求粒子和网格的距离向量时,乘上各项异性矩阵G即可。

到这里我们的理论知识就完成了,下一篇则分析代码如何实现。

  • 1
    点赞
  • 5
    收藏
    觉得还不错? 一键收藏
  • 3
    评论
评论 3
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值