CN113158531B - 一种利用形变梯度的单组分与多组分不可压缩流体仿真方法 - Google Patents

Abstract

一种利用形变梯度的单组分与多组分不可压缩流体仿真方法，利用形变梯度，计算单组分或多组分粒子的体积变化量，并基于形变梯度建立流体不可压缩条件的方程并用松弛的雅各比迭代方法求解。本发明方法实现了对多组分流体的不可压缩求解，实现了真实的多组分流体仿真效果，并提供了艺术化控制的方案，对单组分流体的仿真效率和效果与当前业界先进方法的效果可比。

Description

一种利用形变梯度的单组分与多组分不可压缩流体仿真方法

技术领域

本发明涉及计算机图形学物理模拟领域，尤其涉及一种单组分与多组分不可压缩流体的动力学仿真方法与系统。

背景技术

不可压缩的流体仿真是流体动力学的研究重点，它对于流体仿真的视觉真实性，物理可靠性都有重要影响。多组分的流体仿真，是流体仿真中重要的研究领域。

对于不可压缩性的流体仿真，前人有着长时间的研究。借鉴于基于网格的不可压缩流体求解器，早期方法，如S.J.Cummins and M.Rudman在“An sph projection

Method”一文中，通过直接将光滑粒子映射到离散的网格上，并求解不可压缩条件，以维持流体的不可压缩性。但是这种方法对采样方式敏感，而且尺度伸缩性差。近年来，M.Ihmsen等人在“Implicit incompressible sph”提出的隐式不可压缩的sph方法，改进了离散化方式，使得整个方案在并行图形处理单元上可以运行。

在多组分的流体仿真领域，前人也有出色的研究成果。如B.Ren等人在“Multiple-fluid sph simulation using a mixture model”提出的使用混合的SPH模型对多组分流体进行仿真，实现了多组分流体的分层效果。但是由于混合流体的密度不是常量，速度场也不是无散的，该方法难以保证流体的可压缩性。

发明内容

本发明提出了一种新颖的保证流体仿真中流体不可压缩性的算法，借助了形变梯度(deformation gradient)来衡量物理量的体积变化。提出的该算法框架不仅可以实现传统的单组分不可压缩流体仿真任务，还可以实现新的多组分不可压缩流体仿真任务。该算法框架在实现了高真实感的流体仿真的同时，保证了流体运动的高效求解。

本发明的技术方案：

一种利用形变梯度的单组分与多组分不可压缩流体仿真方法，本发明方法可以进行下述流体仿真任务：

1)对单组分流体进行基于光滑粒子动力学方法的不可压缩流体仿真；

2)对多组分流体进行基于光滑粒子动力学方法的不可压缩流体仿真；对于单组分或多组分流体的仿真步骤如下：

第1步：初始化流体场景的粒子参数，包括位置，速度，各组分体积分数，密度；

第2步：在每个时间步开始，根据当前场景中的各粒子参数，预测计算下一时刻的粒子位置、速度，并计算出与预测对应的形变梯度及其雅可比行列式；对于多组分流体仿真任务，额外提前根据混合流模型的质量方程计算粒子的组分与密度变化及每组分的偏移速度；

第3步：利用第2步的计算结果，通过使用松弛雅可比迭代方式求解一个线性方程组的形式，求解获得一个当前时间步内的粒子压强修正值；所述求解线性方程组的步骤如下：

第3.1步：对于每个粒子，利用第2步中的预测计算结果各项，推导出该粒子压强修正值与邻域粒子压强修正值需满足的方程等式，联立所有粒子需满足的方程等式构成一个线性方程组；构建线性方程组的方法如下：

第3.1.1步：借助形变梯度及雅可比行列式计算各流体粒子在仿真过程中的体积变化；

第3.1.2步：对于固体边界处的仿真，在流体粒子邻域存在边界固体粒子时，对所述的线性方程组进行修正以包含边界固体粒子的影响；该修正不改变方程的整体数学形式和求解方法；

第3.2步：对于单组分或多组分流体模拟任务，所求解的线性方程组及其求解方法具有统一的数学形式，二者的区别仅在于少量的具体项的计算方法不同；

第3.3步：从一个初始估计压强修正值开始，根据松弛雅可比方法不断迭代获得更好的压强修正值，直到仿真中所有粒子的方程左右两边的误差的最大值小于给定最小量，或迭代次数超过给定上限；

第4步：利用第3步获得的压强修正值修正预测计算中的粒子位置，速度，获得不可压缩流体的仿真结果；

第5步：重复第2-4步，获得一个时间序列中每时间步对应的不可压缩流体的仿真结果。

本发明的优点和有益效果：

本发明的创新点包括：

1)基于形变梯度的衡量流体形变的技术

物质物理上的压缩性直接与物质的体积相关，这个体积变化可以用形变梯度张量来衡量，它的Jacobian式衡量了体积的变化率。形变梯度的概念起源于固体动力学仿真，用来描述运动过程中的不可压缩固体。本发明引入了形变梯度的概念，并提出了新颖的仿真方案，来保证基于SPH方法的多流体仿真中的不可压缩性。

2)多流体混合物仿真的不可压缩增强技术

本发明提出了一种新颖的保证多流体混合物不可压缩性的技术。借助形变梯度，直接测定流体混合物中，由于物理性和非物理性导致的流体体积变化，从而推导出关于使得流体不可压缩的压强修正值的线性方程组。这项新技术可以被直接应用在现有单流体和多流体混合物的仿真框架中，并实现真实的视觉效果以及高效的计算。与现有的多流体混合物仿真框架相比，本发明能够获得不可压缩的仿真结果；对单流体仿真来说，本发明与现有先进框架的计算性能相当。

综合而言，本发明从理论上，将流体运动中的流体密度，速度散度与非物理性的压缩进行解耦。本发明大大地增强了多流体仿真中的视觉真实性，避免了前人的技术方案中流体压缩的问题。本发明应用在多流体仿真中，在避免了原有求解器的复杂的同时，仿真出更加逼真的流体场景；应用在单流体仿真过程中，本发明可以集成到现有的流体仿真算法框架中，而无需引入过多的计算量。

附图说明

图1为本发明利用相场理论进行二维与三维晶体生长模拟的系统方框图。

图2是二维5角对称情形不同对称破缺的模拟效果。

图3是三维结晶生长结果及其对应的对称函数的图示。

图4是艺术控制模拟效果。

图5是艺术控制模拟效果。

图6是使用本发明方法恢复的二维晶体三角网格输出示意图。

图7是本发明方法在曲面流型上的效果。

图8展示了本发明方法相较对比方法能更好地实现不可压缩性。

图9展示了本发明方法相较对比方法能降低粒子杂乱度，提高分层效果。

具体实施方式

1.初始化

我们的方法可以使用现有标准的粒子方法仿真方案的通用初始化方法。对每个粒子，存储它的位置，速度，各组分体积分数，密度。

2.单/多组分流体仿真方法中的预测计算

根据标准的单组分流体纳维-斯托克斯动量方程(参见下公式2)和多组分流体混合流模型纳维-斯托克斯方程(参见下公式15)的显示求解方法，根据当前时间步开始时所有粒子上存储的参量，进行显式地计算，获得一个预测的流体粒子加速度，并进一步根据该加速度，使用标准流体仿真的平流算法获得预测的下一时间步该粒子的位置、速度，并计算出与预测对应的形变梯度及其雅可比行列式；对于多组分流体仿真任务，额外提前根据混合流模型的质量方程(公式14)计算粒子的组分与密度变化及每组分的偏移速度；

3.单/多组分流体仿真方法中所需求解的线性方程组的推导过程、形式及求解方法

先以单流体为例，推导并解释本发明方案的基本思想，然后再扩展到多流体仿真应用上。

3.1.借助形变梯度衡量体积变化

以单流体成分的仿真为例，流体运动由质量方程和动量方程来描述：

其中，ρ，p，u，M分别表示密度，压力，速度和流体受到的非压力外力。通常情况下，单流体的不可压缩性是通过恒定的密度来实现的，通过公式1可以推出流体无散速度的条件，也就是

基于网格的解法里，传统的求解无散条件的方法通过求解带条件的动量方程。但是在基于SPH的方法里，由于粒子的位置在不停发生变化，所以一种高效的方式是使用当前时间步已知的周围粒子的数值来求解不可压缩条件。

一般来说，在流体仿真中，视觉上的可压缩效果通常是由两个原因导致的，一个是由于弱可压缩条件或者求解器中的数值误差；另一个是由于流体物理属性发生变化导致的体积变化，例如多流体混合物的反应，溶解等。因此，多流体仿真模型的设计应该允许由于物理现象导致的流体密度变化，速度散度的同时，消除非物理性原因导致的压缩效果。

本发明方案设计了一个新颖的保证不可压缩的SPH仿真技术，通过当前时间步的粒子拓扑计算形变梯度，并用形变梯度来衡量流体体积上的变化。形变梯度是通过形变位置与初始时间步位置的梯度计算的，公式如下：

其中F是形变梯度张量，I表示单位矩阵，x和u是形变物的位置和速度，

表示物质梯度运算符，下标0表示相较于初始帧的梯度。

形变梯度的雅各比式

衡量了体积的变化，其中V和V₀分别表示形变后和原始的元素体积。|F|表示F矩阵的行列式。因此，对于恒定密度的流体来说，要保证不可压缩性就相当于在仿真过程中，设置：J＝1。对于多流体的仿真来说，J指示了流体粒子局部体积的变化量，直接衡量了局部粒子体积变化是否与物理现象导致的体积变化一致。

3.2.单组分流体线性方程组

用x^adv，v^adv和F^adv分别表示预测的粒子位置，速度及形变梯度，这三个物理量是通过已知的当前帧的物理量，向前平流一个时间步得到的预测量。然后我们计算一个使J满足不可压缩条件的压力校正项Δp，来更正预测量。根据公式2和SPH速度更新方案，Δp会带来的速度改变量u^p为：

其中，

NF(i)表示粒子i的邻域，m表示粒子的质量，

表示插值密度，W为核函数。速度上的改变量u^p会连带着改变公式3中定义的F，改变量为：

由于粒子的质量守恒，可知

目标密度可以设置为ρ_t+1＝ρ₀或其它可变量。以下公式7中，我们设置

为当前流体密度并设置目标密度为ρ_t+1＝ρ₀，可得下式：

在后续推导中，任一项都默认是使用t时刻的数值进行计算的，为了简便，下文中的下标t将省略，下标t+1会明确写出。公式7右侧可以通过下式计算：

J＝|F^adv|+d|F^adv|＝|F^adv|+tr(F^adv*dF^adv) 公式8

其中F^adv*是F^adv的伴随矩阵。我们使用SPH的离散化方法来计算

其中N_i是一个3×3矩阵，其中X行Y列的元素为：

其中，对于一个物理量Q，Q_ji＝Q_j-Q_i，

和

分别指对应向量的第X，Y个元素。根据公式6-9，对于每个粒子i，我们都可以计算得到一个关于速度的线性方程：

基于上述式子，我们可以进一步替换u^p中的Δp。首先，为了方便起见，可以重写下式：

将公式11代入到公式10中，可以得到如下关于Δp的线性方程：

使用松弛雅各比方法，如下式迭代地求解Δp：

其中ω是松弛因子。

3.3.多组分流体线性方程组

对于多组分流体，求解混合物的SPH模型需要求解下列公式：

其中，ρ_m是粒子总密度，α_k是流体中第k组分的体积分数，u_m，u_mmk分别是粒子速度和第k组分的漂移速度，M_m表示导致粒子速度改变的力。组分速度为：u_k＝u_m+u_mk。

多组分流体中，每个粒子的组成处于不断变化之中，要为一个流体粒子定义J值并不容易，取而代之的是，我们为每个组分定义J_k和粒子总密度ρ_m＝∑_kα_k，t+1ρ_k，来推导不可压缩求解方法。首先，混合粒子的密度可以表示为：

其中J_k是第k组分的形变梯度张量的雅各比式。公式16的第二个等式是通过假设各个组分都是均匀地分布在SPH粒子上得到的。粒子上的一个元体积中，每个组分的体积分数的升高使得粒子中该组分的比例升高，反之亦然。这就相当于流体局部区域的组分体积的升高或收缩。因此，

接下来，通过求解混合流模型中的质量方程(公式14)，可以得到α_k，t+1并且将ρ_m，t+1在后续计算中作为已知值。通过类似的步骤，可以计算得到

计算得到的数值用来估算机械运动导致的形变量J_k的变化。

通过一阶泰勒展开，

公式16改写为：

通过公式6和公式9中的u^p替换为

可以得到：

上式与单组分仿真的对应公式10相比具有完全一致的形式。进一步地推导得到公式12变为：

公式13变为：

比较公式12与公式19，公式13与公式20可见，单组份流体与多组分流体的线性方程组与其求解方法具有统一的数学形式。

3.4.边界处理

本发明采用了与Versatile rigid-fluid coupling for incompressible sph类似的流固耦合方案。边界粒子对流体粒子i施加的压力为：

与单组分流体模拟相比，通过用ψ中的ρ_mi＝∑_kα_kiρ_ki替换ρ_0i，得到：

对于单流体求解器，将公式22代入到公式10可以得到：

对于多组分流体仿真，经过完全类似地推导可得，可将公式23式左边替换为公式19左边，并用G_i替换

即得到适用于多流体仿真的公式。

4.修正预测计算，保证流体不可压缩性

在第3步中已经计算出Δp后，重新进行第2步的显式计算过程，但在此时用到粒子压强p的时候以p+Δp作为取值，即利用Δp修正原先的p值。这样获得的下一时间步的位置、速度即可保证流体的不可压缩性。

图1第一行展示了本发明方案的效果，第二行展示的是Multiple-fluid sphsimulation using a mixture model中方法的结果。展示的场景是三种不可互溶的液体的溃堤场景。可以从黄线的指示看出，本发明的方法更好地保证了流体体积上的不可压缩性。

图2展示仿真过程中测量到的平均密度误差。本发明方案显著地减少了密度误差，并实现了更好的视觉效果。本发明方案的仿真效果，使得流体体积得到更好的维持，流体运动过程中的外观更有序。

图3设置了一个立方体形状的流体掉入一个由两种流体构成的水池中，四种组分的液体可互溶。第一行展示了本发明方案的效果，本发明方案实现了具有真实感的不可压缩流体仿真。第二行展示了对比方法的结果，可以观察到，流体压缩性导致了粒子边界的剧烈动荡，此外，可压缩性和剧烈的粒子动荡导致了混合后不自然的颜色。第三行展示了使用对比方法但使用了1/100的时间步大小和100倍的stiffness。可以观察到，通过调小时间步，流体整体体积得到了维持，但是更大的stiffness导致了边界动荡问题更加严重，甚至会使仿真过程崩溃，粒子边界动荡同样地导致了不同流体组分的不自然的混合。在实验统计中，第二行的仿真过程中，有48.7％-41.2％的体积损失，第三行的仿真过程中，通过使用更小的时间步，体积损失减少到21.9％-9.1％，但是第三行的结果仍然存在着明显的不自然的痕迹：剧烈的粒子间动荡导致了不同组分间流体的不自然的混合。本发明方案提供了自然的不可压缩仿真方案。

图4展示了一个离心力场景。图中的圆柱形容器是完全密封并且完全注满流体的，也就是说，如果流体的不可压缩性能被完美维持，那么容器中不会有任何额外空间存在。图中第一行展示的是本发明方案效果，第二行展示的是对比方案的效果。可见，本发明方案的仿真实现了高度的不可压缩性。而对比方案的仿真结果中，由于流体中出现了体积压缩，导致容器中心出现了明显空缺。

图5的例子中，通过改变形变梯度来进行艺术化效果的控制。也就是说，在不改变组分密度的前提下，通过改变J，粒子体积会相应发生改变。并将下一个时间步的粒子体积设置为期望值，

通过

重写公式16。替换J_m，t为J_contJ_m，t，并将公式18中的ρ_mi，t+1替换为

在仿真过程中，将J_cont设置为小于1的数值会使相应位置的压力值更大，反之亦然。所以通过在场景中S型区域，在仿真过程中途，人为设置J_cont＝2，密度更大的深色组分向上运动，形成一个S型的喷泉，而其它流体区域仍维持平稳。

图6，展示了本发明方法以及IISPH方法对单组分流体仿真的效果。为了比较两种方法的效果，我们设置了一样的仿真参数，并设置了固定的迭代步数。在图6场景中，使用了52000个流体粒子，在一个Intel Xeon 2.6GHz CPU的单核上，本发明方案每个时间步使用1.73s，IISPH方法耗时1.53s。仿真过程中的平均密度误差和最大密度误差展示在图7中。本发明方案与IISPH的收敛速度相近。

图8展示了与对比方案的效果。第一行，设置了更高的stiffness系数，使用了更小的时间步，对比方案的流体仿真仍然开始有出现崩溃的倾向。第二行，使用对比方法并调小时间步为原来1/4，并调大stiffness系数为原来14倍，但是容器中间由于流体压缩仍出现明显空缺。

图9中左侧展示了本发明方案的效果，左侧展示了对比方法的效果。从粒子排布上可以观察到，本发明方案的粒子排布更加规整有序。而右侧中，相邻的粒子之间的颜色差异很大，流体静止后的分层间也存在着无序粒子。本发明方法在仿真过程中，粒子处于有序状态，流体静止后，分层效果更好。

Claims (3)

1.一种利用形变梯度的单组分与多组分不可压缩流体仿真方法，其特征在于，该仿真方法能够进行下述流体仿真任务：1)：对单组分流体进行基于光滑粒子动力学方法的不可压缩流体仿真；2)：对多组分流体进行基于光滑粒子动力学方法的不可压缩流体仿真；对于单组分或多组分流体的仿真步骤如下：

第1步：初始化流体场景的粒子参数，包括位置，速度，各组分体积分数，密度；

第2步：在每个时间步开始，根据当前场景中的各粒子参数，预测计算下一时刻的粒子位置、速度，并计算出与预测对应的形变梯度及其雅可比行列式；对于多组分流体仿真任务，额外提前根据混合流模型的质量方程计算粒子的组分与密度变化及每组分的偏移速度；

第3步：利用第2步的计算结果，通过使用松弛雅可比迭代方式求解一个线性方程组的形式，求解获得一个当前时间步内的粒子压强修正值；具体方法是，使用松弛雅各比方法，如下式迭代地求解Δp：

其中上标s,s+1代表第s,s+1次迭代；下标i代表粒子编号；ω是松弛因子，

为当前流体密度，ρ₀为目标密度，F^adv表示形变梯度，a_ii是自修正压强系数，dterm(Δp)是邻域线性修正压强；

对于多组分流体，按如下式迭代地求解Δp：

其中ρ_mi，t+1代表下一时间的预测流体密度，

为本时间的组分体积分数，

为形变梯度的雅可比行列式；

第4步：利用第3步获得的压强修正值修正预测计算中的粒子位置，速度，获得不可压缩流体的仿真结果；

第5步：重复第2-4步，获得一个时间序列中每时间步对应的不可压缩流体的仿真结果。

2.根据权利要求1所述利用形变梯度的单组分与多组分不可压缩流体仿真方法，其特征在于，第3步所述求解线性方程组的步骤如下：

第3.1：对于每个粒子，利用第2步中的预测计算结果各项，推导出该粒子压强修正值与邻域粒子压强修正值需满足的方程等式，联立所有粒子需满足的方程等式构成一个线性方程组；

第3.2：从一个初始估计压强修正值开始，根据松弛雅可比方法不断迭代获得更好的压强修正值，直到仿真中所有粒子的方程左右两边的误差的最大值小于给定最小量，或迭代次数超过给定上限。

3.根据权利2所述利用形变梯度的单组分与多组分不可压缩流体仿真方法中，其特征在于：第3.1步的构建线性方程组的方法如下：

第3.1.1步：借助形变梯度及雅可比行列式计算各流体粒子在仿真过程中的体积变化；

第3.1.2步：对于固体边界处的仿真，在流体粒子邻域存在边界固体粒子时，对权利要求2中的线性方程组进行修正以包含边界固体粒子的影响；该修正不改变方程的整体数学形式和求解方法。

Images