基于动态压力校正的高粘性流体流变模拟方法

技术领域

本发明涉及高粘性流体流变技术领域，特别涉及一种基于动态压力校正的高粘性流体流变模拟方法。

背景技术

SPH(Smoothed Particle Hydrodynamics)方法一种流体仿真方法，这种仿真的方法主要在游戏交互、电影特效、CFD数值模拟、天文学研究、海洋学研究等多种领域有着广泛的应用。经过很多的年的发展历程，该方法已经能够实现不同流体的特征模拟，以及多相流耦合交互的模拟，并且能够很好的满足模拟场景所需的实时性和稳定性。

然而从1977年发明至今，SPH方法一直都存在着，随时间的推移出现密度震荡和体积衰减的问题，原因有以下四个方面：(1)SPH方法中粒子的坐标更新是基于拉格朗日法实现的，而粒子间的作用力项存在欧拉积分离散化的表述，由于核函数的引入，这将会导致带有逆压力梯度的粒子会随着相互距离的减小而吸引力更强，当前流场的密度则会增大；(2)SPH方法中引入的密度误差来校正压强，静态的校正系数将会导致每次校正的量都会超调，从而引起较大的稳态误差；(3)SPH方法的拉格朗日连续方程中唯独只对密度的误差做了限制，而无散度约束并没有相关的公式加以限制，这将导致模拟随时间的推移，体积衰减的问题。(4)传统SPH方法中无法克服复杂模型和形变模型的流固耦合中，边界条件的角动量损失：当边界粒子生成的不足时引起的压力和粘性力的不对称，从而造成角动量损失，视觉上将反映出部分粒子率先接触到边界时，靠近边界的粒子会因为压力的不对称，受到上层粒子的挤压，而永远停留在边界处无法回弹，粘性力不对称对导致边界处的速度损失，该现场通常被称为虚假的视觉效果。

发明内容

基于此，本发明的目的是提出一种基于动态压力校正的高粘性流体流变模拟方法，通过设置密度恒定约束与无散度约束共同作用的动态校正系数对传统的SPH方法进行改进，以解决现有SPH方法的不足，从而显著提高了流体仿真的稳定性。

根据本发明提出的一种基于动态压力校正的高粘性流体流变模拟方法，所述方法包括：

初始化高粘性流体的属性，所述属性包括速度、坐标、邻域矩阵、密度、校正因子,并根据速度场的拉普拉斯近似与隐式线性矩阵方程组的组合算法以及雅可比迭代进行隐式粘性力的解算，得到非压力下的速度模拟结果；

在压力校正项引入之前，根据CFL条件计算得到高粘性流体的全局时间步长；

将预测密度与静态密度的残差作为第一源项，并根据所述第一源项来回馈控制密度恒定约束条件的校正系数，以建立密度恒定约束条件下的密度恒定约束模型，以根据密度恒定约束模型预测得到经过密度恒定约束的压力校正项后的速度场；

根据经过密度恒定约束的压力校正项后的速度场更新一次粒子的位置，并基于更新后的位置更新速度场属性，速度场属性包括邻域矩阵、密度和校正因子；

将更新后的密度的随体导数作为第二源项，并利用第二源项来反馈控制无散度约束条件的校正系数，从而建立无散度约束的无散度约束模型，以根据无散度约束模型预测得到经过无散度约束的压力校正项后的速度场；

当模拟时间大于预设时间阈值时，则输出速度场，并根据速度场和全局时间步长来更新下一次的流场粒子坐标。

优选地，所述根据速度场的拉普拉斯近似与隐式线性矩阵方程组的组合算法以及雅可比迭代进行隐式粘性力的解算，得到非压力下的速度模拟结果的步骤包括：

根据以下公式计算得到粘性力项对应的速度预测：

其中，

根据以下公式计算得到n+1时刻的流体粒子i的速度拉普拉斯：

其中，d为空间维度，m

优选地，所述所述根据速度场的拉普拉斯近似与隐式线性矩阵方程组的组合算法以及雅可比迭代进行隐式粘性力的解算，得到非压力下的速度模拟结果的步骤还包括：

根据以下公式进行雅可比迭代迭代求解，以收敛出下一时刻粘性力项引起的速度：

其中

其中，Δt表示全局时间步长，A表示速度场拉普拉斯的邻域矩阵，A

优选地，所述根据CFL条件计算得到高粘性流体的全局时间步长的步骤包括：

根据以下公式计算全局时间步长：

其中，v

优选地，所述将预测密度与静态密度的残差作为第一源项，并根据所述第一源项来回馈控制密度恒定约束条件的校正系数，以建立密度恒定约束条件下的密度恒定约束模型的步骤包括：

根据以下公式计算得到预测密度：

其中，

根据以下公式构建密度恒定约束模型下的校正参数的求解方程：

其中，α

根据以下公式构建密度的随体导数的求解方程：

其中，

根据以下公式确定经过密度恒定约束的压力校正项：

其中，

将压力项带入密度恒定约束模型下的校正参数的求解方程中，以得到密度恒定约束模型下的校正参数：

表示核函数范围内所有流体粒子j的质量权重和，/>

优选地，所述根据密度恒定约束模型预测得到经过密度恒定约束的压力校正项后的速度场的步骤包括：

其中，

优选地，所述将更新后的密度的随体导数作为第二源项，并利用第二源项来反馈控制无散度约束条件的校正系数，从而建立无散度约束的无散度约束模型：

对密度的随体导数的求解方程进行离散化得到：

根据以下公式计算得到无散度约束条件的校正系数：

其中，

优选地，所述根据无散度约束模型预测得到经过无散度约束的压力校正项后的速度场的步骤包括：

根据以下公式计算得到速度场：

其中，

与现有技术相比，本发明具有以下优点：

(1)本发明中粘性力的解算采用了隐式雅可比迭代，并且粘性力项中不携带应变率张量属性，避免了边界粒子的粘性力的速度梯度不为零的问题，从而消除边界的角动量损失与阻尼运动，保证了粒子在高粘性运动下的边界稳定性；

(2)本发明通过密度和散度的动态校正来保证模拟的稳定性，根据CFL条件，在保证动量守恒的前提下，时间步长增大也不会发生散度误差随着时间的积累，可以实现更大的时间步长的模拟，从而降低了模拟所消耗的计算机性能；

(3)本发明选用了预测密度和静态密度的残差作为源项，通过动态密度校正系数来调节流场密度的超调，从而保证最后的稳态误差范围限制在1.2％；

(4)本发明选用了密度的随体导数作为源项，通过动态散度校正系数来限制流场体积的收缩，使流场的体积不会随时间的推移而减少，从而保证了流场的精确性；

(5)本发明选用的动态校正系数中包含的校正因子可以在密度约束和散度约束中同时使用，降低了模拟的消耗，并且作为单粒子属性，便于局部控制；

(6)本发明中存储流场属性的粒子在模拟中分布均匀，可以有效的提高模拟时的收敛速度，提高模拟的稳定性；

(7)本发明的密度和散度约束可以通过GPU硬件设备并行处理，进一步的减少高粘度流体流变模拟过程中的时间开销。

本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实施例了解到。

附图说明

图1为本发明一实施例中基于动态压力校正的高粘性流体流变模拟方法的流程图；

图2为不同粘性系数下的血液流变模拟对比；

图3为流场平均密度结果对比示意图；

图4为流场靠近边界时的角动量对比示意图；

图5为在牙科手术系统中的应用实例图。

如下具体实施方式将结合上述附图进一步说明本发明。

具体实施方式

为了便于理解本发明，下面将参照相关附图对本发明进行更全面的描述。附图中给出了本发明的若干实施例。但是，本发明可以以许多不同的形式来实现，并不限于本文所描述的实施例。相反地，提供这些实施例的目的是使对本发明的公开内容更加透彻全面。

除非另有定义，本文所使用的所有的技术和科学术语与属于本发明的技术领域的技术人员通常理解的含义相同。本文中在本发明的说明书中所使用的术语只是为了描述具体的实施例的目的，不是旨在于限制本发明。本文所使用的术语“及/或”包括一个或多个相关的所列项目的任意的和所有的组合。

请参阅图1，所示为本发明一实施例中基于动态压力校正的高粘性流体流变模拟方法的流程图，该方法包括步骤S01至步骤S06，其中：

步骤S01：初始化高粘性流体的属性，所述属性包括速度、坐标、邻域矩阵、密度、校正因子,并根据速度场的拉普拉斯近似与隐式线性矩阵方程组的组合算法以及雅可比迭代进行隐式粘性力的解算，得到非压力下的速度模拟结果；

需要说明的是，本实施例中的高粘度流体可以是蜂蜜、血液、面糊、糖浆、胶水等物质，首先在正常的大气压强条件下，血液的流变模拟是不可压缩的，且流动具有高粘度的特性，因此需要使用到Navier-Stokes方程作为不可压缩牛顿流体模拟的基础，该流体运动物理模型的构建基于三个假设：第一是流体在空间中的运动是连续的，即在单位时间内的流进和流出是相同的。第二是流体速度场中的各项是可微分的，如密度、压强等。第三是满足牛顿三大守恒定律。N-S方程也可以称作是流体的控制方程，其中包括了连续方程与动量方程：

其中第一项为连续方程，也叫做质量守恒方程，描述的是单位时间内有限控制单元体的流体质量增量等于流出有限控制单元面的流体质量。而式的第二项为动量方程，描述的是流体单元的加速度与受力关系。

在本发明一些可选的实施例中，对N-S方程中的非压力项，即重力项和粘性力项进行解算，重力项为重力加速度g，而粘性力项，在高粘性力血液系统的模拟中，一定需要使用到隐式粘性力的解算，给出的速度场的拉普拉斯近似与隐式线性矩阵方程组的组合算法，来实现隐式粘性力的解算，对高粘度血液的模拟，首先给出粘性力项对应的速度预测：

其中，

再将速度场的拉普拉斯近似带入粘性力项对应的速度预测：

其中，d为空间维度，m

将上述方程写成线性系统的形式，其中下一刻的预测速度设为未知量，通过迭代求解该线性系统，就能收敛出下一时刻粘性力项引起的速度：

其中

其中，Δt表示全局时间步长，A表示速度场拉普拉斯的邻域矩阵，A

最后通过雅可比迭代，在高粘性力血液系统的模拟中，通过隐式积分求解线性方程组，不断迭代直到下一刻的速度与上一刻邻域粒子速度的加权残差小于一定阀值，则跳出循环，预测出下一时刻由粘性力导致的速度，例如下式为坐标迭代的示例：

x

其中，x

参见图2，为改进型隐式粘性力的模拟，当随着μ的增大，粒子会因为粘性系数的增大而缓慢的减速，形成了第三幅和第四副图中粒子堆叠的情况，且当μ增大到0.75时，速度场也不会出现瞬间崩溃的现象。由此看出，在牙科手术中建议使用μ＝0.5的隐式粘性力模拟血液喷射。

图3所示μ＝0.5为时，应用于牙科手术牙龈切割中的血液喷溅模拟，在隐式粘性力的作用下，先喷先出的血液先接触到空气，而后喷射出的血液因为出速度较大，所以经过拉普拉斯平滑后，会与先喷射出来的血液形成堆叠效应，达到了预期的牙科手术血液模拟中高粘性力的模拟效果

步骤S02：在压力校正项引入之前，根据CFL条件计算得到高粘性流体的全局时间步长；

在每个模拟循环过程中更新一次，由两个压力约束共同使用，这一点很大程度上的降低了迭代的计算量。在压力校正项引入之前，计算隐式粘度力、重力等所有非压力提供的加速度。再根据CFL条件实时调整模拟中的时间步长：

其中，v

步骤S03：将预测密度与静态密度的残差作为第一源项，并根据所述第一源项来回馈控制密度恒定约束条件的校正系数，以建立密度恒定约束条件下的密度恒定约束模型，以根据密度恒定约束模型预测得到经过密度恒定约束的压力校正项后的速度场；

为了实现密度恒定压力校正，目标是实现实际密度与静态密度的残差

其中，

其中在初始化之前对校正参数

其中，α

再根据以下公式构建密度的随体导数的求解方程：

其中，

根据以下公式确定经过密度恒定约束的压力校正项：

其中，

将压力项带入拉格朗日连续方程中，并将校正参数

表示核函数范围内所有流体粒子j的质量权重和，/>

得到密度恒定约束的压力校正项后，预测出下一个时间步长的速度场：

其中，

步骤S04：根据经过密度恒定约束的压力校正项后的速度场更新一次粒子的位置，并基于更新后的位置更新速度场属性，速度场属性包括邻域矩阵、密度和校正因子；

其中预先计算的因子α

步骤S05：将更新后的密度的随体导数作为第二源项，并利用第二源项来反馈控制无散度约束条件的校正系数，从而建立无散度约束的无散度约束模型，以根据无散度约束模型预测得到经过无散度约束的压力校正项后的速度场；

需要说明的是，在本发明一些可选的实施例中，在完成了密度偏差校正后，速度场通常不是无散度化的。最后，需要由无散度约束的介入，为了实现无散度压力校正，目标是实现密度的随体导数的最小化，密度的随体导数的计算是通过对拉格朗日连续方程进行离散化来确定的：

而无散度约束的压力校正是通过无散度校正参数来实现的，而无散度校正参数根据以下公式计算得到：

其中，

将密度的随体导数带入无散度校正参数中，即可得到无散度约束的压力校正项，从而预测出下一个时间步长的速度场：

即根据以下公式计算得到速度场：

其中，

参见图4，图4(a)中蓝色线为静态密度曲线，红色线为传统SPH方法的平均密度曲线，随着时间的增加，100帧之前平均密度处于持续增加的状态，100帧之后流体粒子的平均密度呈现不断校正超调的震荡情况，100帧左右超调量为62％，当运行到300帧左右时稳态误差维持在了22％的水平，这会直接导致流体在不同时间段内的体积不等。图4(b)红色线为密度恒定无散度SPH的平均密度曲线，随着时间的推移，100帧之前，增长率为0.013kg/m

参见图5，虚拟牙科手术的血液流变模拟离不开与牙龈组织的耦合作用，传统SPH方法中，因为压力和粘性力的无法克服流固耦合边界条件处的角动量损失，如图5(a)，颜色越深，代表应变率张量的F范数越靠近零。由于边界处粒子缺失，导致应变率张量的F范数不为零。一种基于密度恒定无散度SPH的血液流变模拟方法很大程度上克服边界处角动量损失问题，如图5(b)所示，该实验证明了本发明应用在虚拟牙科手术领域的精确性、鲁棒性和可扩展性。

步骤S06：当模拟时间大于预设时间阈值时，则输出速度场，并根据速度场和全局时间步长来更新下一次的流场粒子坐标。

当模拟时间(程序运行时间)大于认为设置的一个时间阀值时，由无散度约束模型输出的速度场即为最终速度场，通过最终速度场乘以时间步长等于最终坐标，再根据该最终坐标更新完坐标后，再将程序运行时间清零，重复执行上述步骤S01至步骤S06，若模拟时间小于或等于时间阈值时，此时重复步骤S01至步骤S05，直至模拟时间大于时间阈值，以输出最终速度场。

与现有技术相比，本发明具有以下优点：

(1)本发明中粘性力的解算采用了隐式雅可比迭代，并且粘性力项中不携带应变率张量属性，避免了边界粒子的粘性力的速度梯度不为零的问题，从而消除边界的角动量损失与阻尼运动，保证了粒子在高粘性运动下的边界稳定性；

(2)本发明通过密度和散度的动态校正来保证模拟的稳定性，根据CFL条件，在保证动量守恒的前提下，时间步长增大也不会发生散度误差随着时间的积累，可以实现更大的时间步长的模拟，从而降低了模拟所消耗的计算机性能；

(3)本发明选用了预测密度和静态密度的残差作为源项，通过动态密度校正系数来调节流场密度的超调，从而保证最后的稳态误差范围限制在1.2％；

(4)本发明选用了密度的随体导数作为源项，通过动态散度校正系数来限制流场体积的收缩，使流场的体积不会随时间的推移而减少，从而保证了流场的精确性；

(5)本发明选用的动态校正系数中包含的校正因子可以在密度约束和散度约束中同时使用，降低了模拟的消耗，并且作为单粒子属性，便于局部控制；

(6)本发明中存储流场属性的粒子在模拟中分布均匀，可以有效的提高模拟时的收敛速度，提高模拟的稳定性；

(7)本发明的密度和散度约束可以通过GPU硬件设备并行处理，进一步的减少高粘度流体流变模拟过程中的时间开销。

在本说明书的描述中，参考术语“一个实施例”、“一些实施例”、“示例”、“具体示例”、或“一些示例”等的描述意指结合该实施例或示例描述的具体特征、结构、材料或者特点包含于本发明的至少一个实施例或示例中。在本说明书中，对上述术语的示意性表述不一定指的是相同的实施例或示例。而且，描述的具体特征、结构、材料或者特点可以在任何的一个或多个实施例或示例中以合适的方式结合。

以上所述实施例仅表达了本发明的几种实施方式，其描述较为具体和详细，但并不能因此而理解为对本发明专利范围的限制。应当指出的是，对于本领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干变形和改进，这些都属于本发明的保护范围。因此，本发明专利的保护范围应以所附权利要求为准。