CN115248989A - 基于屈服准则约束的粘性流体仿真方法 - Google Patents

Abstract

本公开的实施例公开了基于屈服准则约束的粘性流体仿真方法。该方法的一具体实施方式包括：初始化粘性流体仿真场景；根据隐式流体粒子模型，确定经过时间步长后的粒子速度；通过模拟热量传导过程，以确定经过时间步长后的粒子温度；根据该粒子温度，对该粒子速度进行修正。从而，实现了与温度相关的基于屈服准则约束的粘性流体流动现象仿真，扩宽了仿真类型的范围。

Description

基于屈服准则约束的粘性流体仿真方法

技术领域

本公开的实施例涉及计算机技术领域，具体涉及基于屈服准则约束的粘性流体仿真方法。

背景技术

近年来，基于计算机图形学的流体仿真领域内越来越多的研究人员把目光投向了粘性流体现象的研究中。一些关于粘性流体仿真的方法相继被提出，包括基于光滑粒子流体动力学(Smoothed ParticleHydrodynamics，SPH)方法粘性项的数值求解、引入物质点模型(Material Point Method，MPM)、引入Shape matching约束条件等等。然而，由于粘性流体的运动过程相对于非粘性流体更加复杂多变，运动过程计算繁琐，因此在计算效率、算法稳定性上仍有很大的研究空间。

常见的粘性流体仿真研究分为两个思路，基于流体动力学方程的粘性项数值求解和基于几何模型修正的粘性行为约束。在数值求解实现流体粘性方面，通过修正SPH核函数的粘度项，可以使SPH可以模拟诸如泡沫、蜂蜜等非牛顿粘性流体的行为；基于霍尔姆霍兹自由能的能量方程，控制流体粘性从高能态向低能态的转化，可以近似模拟诸多粘性流体的交互，如多种颜料的混合、蛋清蛋黄的粘性流体行为等；MPM方法用MAC网格求解粘性流体的压强投影，以实现考虑粘性属性的流体现象，也可以仿真蜂蜜、牙膏、奶油等等众多粘性流体的流动。基于几何模型约束的粘性流体行为仿真在求解流体动力学方程时仍遵循非粘性流体的数值模型，转而在流体速度、位移更新前对流体行为施加额外约束实现粘性流体的仿真，例如，通过shapematching约束条件对流体粒子的位置和速度施加粘弹性约束，可以限制粒子的运动范围进而近似模拟粘性流体的运动过程；利用空间自适应的四面体网格可以实现可变粘性系数和对高粘性表面堆叠现象的支持；利用弹簧-质子模型在牛顿流体动力学求解步骤后对粒子间施加额外约束的方法也是仿真粘性流体的一种简洁手段。而无论是基于数值求解的方法还是现存基于几何约束的方法，大都存在计算模型复杂，仿真效率低下的缺点，同时，对于一些粘性流体如血流、蜂蜜等，其粘性的大小与其温度属性紧密相关，而现有方法很少考虑流体粘性与温度的关系。

屈服准则(Yield Criterion)是复杂应力状态下，用于控制材料是否发生塑性变形的判断条件。使用屈服准则作为图形学中的流动材料粒子间相互约束行为的判断简单可行，而且运算流程清晰易实现。粒子表示的粘性流体是一种连续介质材料，应用流体动力学表达式可以描述其运动学行为，假设粘性流体是塑性流(Plastic Flow)且遵循塑性定律，则可以通过屈服准则作为约束条件描述其粘性属性特征。

发明内容

本公开的内容部分用于以简要的形式介绍构思，这些构思将在后面的具体实施方式部分被详细描述。本公开的内容部分并不旨在标识要求保护的技术方案的关键特征或必要特征，也不旨在用于限制所要求的保护的技术方案的范围。

本公开的一些实施例提出了基于屈服准则约束的粘性流体现象仿真的方法，来解决以上背景技术部分提到的技术问题中的一项或多项。

第一方面，本公开的一些实施例提供了一种基于屈服准则约束的粘性流体现象仿真方法，该方法包括：初始化粘性流体仿真场景，其中，上述粘性流体仿真场景包括：粘性流体运动区域、边界和初始条件，上述边界包括半开放边界和封闭边界，上述初始条件包括流体位置、密度、温度和速度；根据隐式流体粒子模型，确定经过时间步长后的粒子速度；通过模拟热量传导过程，以确定经过时间步长后的粒子温度；根据上述粒子温度，对上述粒子速度进行修正。

本公开的上述各个实施例具有如下有益效果：首先，初始化粘性流体仿真场景。然后，根据隐式流体粒子模型，确定经过时间步长后的粒子速度。接着，通过模拟热量传导过程，以确定经过时间步长后的粒子温度。最后，根据上述粒子温度，对上述粒子速度进行修正。由此，实现了对流体温度与流体粘性的相关性的建模。同时，比现有粘性流体仿真模型具有更小的复杂度和计算量，满足了基于物理的粘性流体现象仿真需求。

本发明的原理在于：

本发明提出了一种基于屈服准则约束的粘性流体现象仿真方法，其原理在于：基于粒子表示的粘性流体相较于非粘性流体，在运动过程中表现为粒子之间具有更大的粘滞力，这种粘滞力本质上是一种流体粒子间的相互约束关系，使粘性流体在整体上表现为流动缓慢，凝结现象显著的特征。基于几何模型修正的粘性行为约束方法使用非粘性流体的数值计算模型，可以在进行流体动力学行为求解时具有相对较高的计算效率，在粒子位置和速度更新前，几何约束方法对粒子间的相对位移、应力、速度等表征粘性的属性进行修正，进而可以实现粘性流体的行为建模。塑形屈服准则(Yield Criterion)是复杂应力状态下，用于控制材料是否发生塑性变形的判断条件。粒子表示的粘性流体一种连续介质材料，应用流体动力学表达式可以描述其运动行为，假设粘性流体这种粘性流体是塑性流(Plastic Flow)且遵循Mohr-Coulomb塑性定律，则可以通过塑性屈服准则描述其粘性行为。

同时，粘性流体的粘度属性与其温度具有很高的关联性，粘性流体的温度越高，其构成的分子运动越活跃，宏观上表现为流体的粘性属性特征越不明显；相反，粘性流体的温度越低，其构成的分子运动不活跃，宏观上表现为流体的粘性属性特征明显，即流体温度越低，粘性越大。考虑到温度与粘性属性的对应关系，在屈服准则约束中引入温度属性和控制温度权重的相关系数，可以实现温度敏感的屈服准则约束模型，进而仿真不同温度条件下不同粘性的粘性流体现象。

本发明与现有技术相比的优点在于：

1、本发明提出的基于屈服准则约束的粘性流体现象仿真方法，应用于计算机动画和虚拟现实场景建模领域，创新性的将屈服准则约束引入基于粒子的粘性流体的建模中用以流体粘性的模拟，较于现有的粘性流体仿真方法更加简洁，且易于实现。

2、本发明改进了屈服准则约束条件，通过引入温度特征与权重参数，实现了粘性流体温度属性与粘性属性的关联，可以实现不同温度下不同粘度的流体现象建模。

3本发明提出的基于屈服准则约束的粘性流体现象具有扩展性强，可仿真类型广的优点，诸如血流、蜂蜜、奶油、牙膏等温度敏感的粘性流体都可通过本发明提出的方法进行仿真实现。

附图说明

结合附图并参考以下具体实施方式，本公开各实施例的上述和其他特征、优点及方面将变得更加明显。贯穿附图中，相同或相似的附图标记表示相同或相似的元素。应当理解附图是示意性的，元件和元素不一定按照比例绘制。

图1是根据本公开的一些实施例的基于屈服准则约束的粘性流体仿真方法的流程图；

图2是根据本公开的基于屈服准则约束的粘性流体仿真方法的一些实施例的流程图；

图3是本公开的基于屈服准则约束的粘性流体仿真方法的温度映射过程示意图；

图4是本公开的基于屈服准则约束的粘性流体仿真方法的粘性流体粒子的热传导示意图；

图5是20摄氏度下粘性流体行为效果图；

图6是80摄氏度下粘性流体行为效果图。

具体实施方式

下面将参照附图更详细地描述本公开的实施例。虽然附图中显示了本公开的某些实施例，然而应当理解的是，本公开可以通过各种形式来实现，而且不应该被解释为限于这里阐述的实施例。相反，提供这些实施例是为了更加透彻和完整地理解本公开。应当理解的是，本公开的附图及实施例仅用于示例性作用，并非用于限制本公开的保护范围。

另外还需要说明的是，为了便于描述，附图中仅示出了与有关发明相关的部分。在不冲突的情况下，本公开中的实施例及实施例中的特征可以相互组合。

需要注意，本公开中提及的“第一”、“第二”等概念仅用于对不同的装置、模块或单元进行区分，并非用于限定这些装置、模块或单元所执行的功能的顺序或者相互依存关系。

需要注意，本公开中提及的“一个”、“多个”的修饰是示意性而非限制性的，本领域技术人员应当理解，除非在上下文另有明确指出，否则应该理解为“一个或多个”。

本公开实施方式中的多个装置之间所交互的消息或者信息的名称仅用于说明性的目的，而并不是用于对这些消息或信息的范围进行限制。

下面将参考附图并结合实施例来详细说明本公开。

图1示出了根据本公开的一些实施例的基于屈服准则约束的粘性流体仿真方法的流程图。图1给出了基于屈服准则约束的粘性流体现象仿真方法的总体处理流程，下面结合其他附图及具体实施方式进一步说明本发明。

继续参考图2，示出了根据本公开的基于屈服准则约束的粘性流体仿真方法的一些实施例的流程200。该基于屈服准则约束的粘性流体仿真方法，包括以下步骤：

步骤201，初始化粘性流体仿真场景。

在一些实施例中，基于屈服准则约束的粘性流体仿真方法的执行主体可以通过接收预先设定的初始参数来初始化粘性流体仿真场景。其中，上述粘性流体仿真场景可以用于对粘性流体进行流体模拟。上述粘性流体仿真场景可以包括：粘性流体运动区域、边界和初始条件，上述边界可以包括半开放边界和封闭边界，上述初始条件可以包括流体位置、密度、温度和速度。上述粘性流体运动区域可以是粘性流体活动的空间范围。上述粘性流体可以包括至少一个流体粒子。上述半开放边界可以是开放区域的边界。上述封闭边界可以是封闭区域的边界。上述流体位置可以是流体粒子在上述粘性流体仿真场景中的位置。

作为示例，上述粘性流体可以是血液。上述粘性流体运动区域可以是杯子内的空间范围。

步骤202，根据隐式流体粒子模型，确定经过时间步长后的粒子速度。

在一些实施例中，上述执行主体可以通过求解隐式流体粒子模型，确定经过时间步长后的粒子速度。其中，上述隐式流体粒子模型可以是FLIP(Fluid implicitparticles，流体隐式粒子)模型。上述粒子速度可以是流体粒子的速度。

可选地，上述执行主体根据隐式流体粒子模型，确定经过时间步长后的粒子速度，可以包括以下步骤：

第一步，将上述粘性流体仿真场景中的初始条件包括的速度插值到三维网络上。其中，上述三维网络可以包括至少一个三维网格。上述三维网络中的三维网格与流体粒子在上述粘性流体仿真场景中的流体位置相对应。上述三维网络可以是三维空间。上述三维网格可以是用网格组成的一个三维物体的形状。上述三维网络中的各个三维网格的大小相等。上述执行主体可以将流体粒子在仿真场景内的速度投影到上述三维网络中的对应位置。

第二步，通过求解以下方程，确定经过时间步长后的流体粒子在上述三维网格上的速度：

其中，t表示时间。ρ表示在t时刻流体的密度。u表示流体粒子在t时刻的速度。p表示在t时刻流体的预设压强。f表示在t时刻流体所受的外力。

第三步，将上述流体粒子在上述三维网格上的速度与上述粘性流体仿真场景中的初始条件包括的速度的差确定为速度变化量。

进而，根据上述隐式流体粒子模型的插值方法将上述速度变化量插值回流体粒子。

第四步，通过以下公式，确定粒子速度：

v＝αv_FLIP+(1-α)v_PIC。

其中，v表示上述粒子速度。α表示第一权重。且α取值范围为[0，1]。v_FLIP表示根据上述隐式流体粒子模型求得的速度。v_PIC表示根据PIC(particle in cell，质点网格法)方法求得的速度。

其中，上述第一权重表征根据上述隐式流体粒子模型求得的速度所占的比重。

步骤203，通过模拟热量传导过程，以确定经过时间步长后的粒子温度。

在一些实施例中，上述执行主体可以通过模拟热量传导过程，以确定经过时间步长后的粒子温度。其中，上述热量传导过程可以是热量传递过程。

可选地，上述执行主体通过模拟热量传导过程，以确定经过时间步长后的粒子温度，可以包括以下步骤：

第一步，将上述粘性流体仿真场景中的初始条件包括的温度插值到上述三维网格上。

第二步，通过求解以下方程，模拟热量传导过程，以确定经过时间步长后粒子在上述三维网络中的三维网格的温度：

其中，T表示温度。b表示热传导模型的热扩散系数。t表示时间。Δt表示时间步长。x表示在上述三维网格中的网格点的坐标的横坐标，y表示在三维网格中的网格点的坐标的纵坐标。Z表示在三维网格中的网格点的坐标的第三维坐标。

第三步，将经过时间步长后流体粒子在三维网络中的三维网格的温度与上述粘性流体仿真场景中的初始条件包括的温度的差确定温度改变量。

第四步，根据上述隐式流体粒子模型，将上述温度改变量插值回粒子。

第五步，通过以下公式，确定流体粒子的粒子温度：

TN＝αF+(1-α)P。

其中，TN表示粒子温度。α表示上述第一权重。且α取值范围为[0，1]。F表示根据上述FLIP模型方法求得的温度。P表示根据上述PIC方法求得的温度。

步骤204，根据粒子温度，对上述粒子速度进行修正。

在一些实施例中，上述执行主体可以将上述粒子温度作为输入参数，根据屈服准则约束控制方程计算流体粒子所受的摩擦应力，从而通过应力对切向速度施加的增量近似实现流体粒子之间粘性属性的仿真，以实现对上述粒子速度进行修正。

可选地，上述执行主体根据粒子温度，对上述粒子速度进行修正，可以包括以下步骤：

第一步，通过以下公式，确定第一特征，其中，上述第一特征可以是上述三维网络中每个三维网格的应变率张量：

其中，D表示上述第一特征。u表示流体粒子在t时刻的速度。

表示梯度。

表示上述梯度

的转置。

第二步，根据上述第一特征，通过以下公式，确定摩擦应力：

其中，σ_f表示上述摩擦应力。p表示压力。D表示第一特征。|D|_F表示第一特征D的Frobenius范数。

第三步，响应于流体粒子的位置在场景内部，根据上述摩擦应力和上述粒子温度，通过以下公式，对流体粒子的粒子速度进行修正：

其中，u表示流体粒子的粒子速度。σ_f表示上述摩擦应力。

表示用中心差分法计算得到的滑动摩擦力σ_f的散度。β表示预设权重系数。

第四步，响应于流体粒子的位置在场景边界上，根据流体粒子的粒子速度，通过以下公式，修正流体粒子的切向速度：

其中，UT表示上述切向速度。μ表示摩擦系数。n表示法线。|u·n|表示法向速度的模。

实践中，本发明提出的基于屈服准则约束的粘性流体现象仿真方法，具体实施为基于FLIP模型的流体动力学仿真和基于屈服准则约束的粘性行为修正，主要步骤介绍如下：

1、流体动力学仿真建模

为了计算流体粒子的属性和运动过程，采用基于粒子的方法对离散的N-S方程进行求解，N-S方程包含如上述公式两个重要的方程形式。其中，第一个公式

称之为连续性方程，主要作用是保持流体的质量守恒，第二个公式

称之为流体的动量方程，表示流体速度在压强、粘性力和外力的共同作用下随时间的变化规律。FLIP模型作为一种粒子-网格混合模型，本质上也是通过求解N-S方程实现流体动力学行为的仿真。与直接求解N-S方程的欧拉网格方法不同的是，FLIP方法中流体的表现基于离散的粒子模型，流体粒子的属性首先投影到网格中进行求解，即FLIP通过网格求解N-S方程，然后将速度的变化量从网格插值回流体粒子，进而驱动流体粒子的运动。FLIP方法由PIC方法发展而来，不同之处在于PIC方法直接将求得的速度值插值回流体粒子，相较而言，FLIP方法只传递速度变化量的方式，避免了误差的累积，比PIC方法具有更高的精度。通常，使用PIC和FLIP的速度加权平均值作为流体粒子的新速度既可以保证流体仿真的稳定，也可以尽量减少插值误差的累积。

2、热传导过程仿真

基于热传导模型的粘性流体现象仿真可以表现不同温度条件对流体粘性的影响，得益于FLIP模型中粒子属性的计算全部基于网格求解，非常适合与基于网格求解的简化热传导模型进行结合。

为了模拟热传导过程，需要为所有粒子额外增加一项温度属性，在场景初始化时给定一个温度值，温度的改变基于热传导模型的计算结果。在每个时间步长内，网格格点的温度由粒子插值而来。

在网格格点处更新过温度值之后，将温度的变化映射回粒子(映射过程如图3所示)。然后对粒子的当前温度进行更新，温度的更新规则参考FLIP模型方法与PIC方法结合的模式。热量在不同粒子间的热量传导过程如图4所示，展示了流体粒子的热传导过程，高温流体冲击固体模型，由于固体模型温度较低导致流体粒子温度逐渐降低，同时流体粒子在相互碰撞中也有热量传导，为了清晰地表现出热传导的过程，图中的粒子模型根据温度高低进行着色，颜色信息表示粒子的不同温度，其中浅色粒子表示高温粒子，深色粒子表示低温粒子。

3、粘性属性仿真

为了模拟粉末材料的沙子状态，本发明在FLIP流体求解器的压强投影步骤之后增加了一个用于对粒子摩擦力和塑性进行控制进而模拟粘性属性的步骤。

首先，使用标准的中心差分法(Standard Central Differences)来估算网格中每个单元格的3×3应变率张量D。

通常，通过以下公式计算流体粒子流动时的摩擦应力：

其中，σ_f表示上述摩擦应力。

表示摩擦角度，代表粒子材料静止堆积时的最大坡度，

越小意味着堆积的状态越扁平。p表示压力，是压强投影步骤中计算得到的压力，其梯度表示对速度的影响。D表示上述应变率张量。|D|_F表示上述应变率张量D的Frobenius范数。

对粘性流体而言，只需考虑粘性流体的流动运动，即认为流体粒子无法形成静态堆积效应，粒子间切向应力权重最大，故摩擦角度

取90°。

接着，对于所有粘性流体单元格，根据已经求得的温度TN，对速度进行更新。上述预设权重系数是可以控制温度权重的可变系数，上述预设权重系数越大，温度对粒子间粘性影响越大，温度越高，流体粒子间的粘性也越大。

粘性流体内部的速度处理完以后，还需将边界条件换成摩擦边界条件，通过这种方式来施加粘性流体粒子与外部物体(如墙壁、障碍物)的摩擦。本发明只对法向速度指向物体内部(也就是粘性流体内部)的交界处单元格进行摩擦处理，修正其切向速度。

以上步骤进行完成后，进入下一个仿真的时间步长内，重复以上步骤，实现与温度相关的基于屈服准则约束的粘性流体流动现象仿真。

为证明本发明在计算机动画领域的正确性和有效性，设计了一个封闭边界的三维流场场景，粘性流体模型在重力作用下自由落体掉入场景的地板上。作为对比，图5与图6中均采用本发明提出的基于屈服准则约束的流体仿真方法，图5中流体粒子初始温度为20℃，图6中流体粒子初始温度为80℃。通过这两个图可以明显看出，本发明的方法能够有效的对粘性流体现象进行仿真，同时，不同温度下流体的粘性表现具有明显差异，说明本发明能够基于屈服准则约束对粘性流体的行为进行可视化仿真，并且能够模拟不同温度下不同流体的粘性表现。

本发明未详细阐述的技术内容属于本领域技术人员的公知技术。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发明创造均在保护之列。

Claims (5)

1.一种基于屈服准则约束的粘性流体仿真方法，包括：

初始化粘性流体仿真场景，其中，所述粘性流体仿真场景包括：粘性流体运动区域、边界和初始条件，所述边界包括半开放边界和封闭边界，所述初始条件包括流体位置、密度、温度和速度；

根据隐式流体粒子模型，确定经过时间步长后的粒子速度；

通过模拟热量传导过程，以确定经过时间步长后的粒子温度；

根据所述粒子温度，对所述粒子速度进行修正。

2.根据权利要求1述的方法，其中，所述根据隐式流体粒子模型，确定经过时间步长后的粒子速度，包括：

将所述粘性流体仿真场景中的初始条件包括的速度插值到三维网络上，其中，所述三维网络包括至少一个三维网格，所述三维网络中的三维网格与流体粒子在所述粘性流体仿真场景中的流体位置相对应；

通过求解以下方程，确定经过时间步长后的流体粒子在所述三维网络上的速度：

其中，t表示时间，ρ表示在t时刻流体的密度，u表示流体粒子在t时刻的速度，p表示在t时刻流体的预设压强，f表示在t时刻流体粒子所受的外力；

将所述流体粒子在所述三维网络上的速度与所述粘性流体仿真场景中的初始条件包括的速度的差确定为速度变化量；

根据所述隐式流体粒子模型的插值方法，将所述速度变化量插值回流体粒子；

通过以下公式，确定粒子速度：

v＝αv_FLIP+(1-α)v_PIC，

其中，v表示所述粒子速度，α表示第一权重，且α取值范围为[0，1]，v_FLIP表示根据所述隐式流体粒子模型求得的速度，v_PIC表示根据PIC方法求得的速度。

3.根据权利要求2所述的方法，其中，所述通过模拟热量传导过程，以确定经过时间步长后的粒子温度，包括：

将所述粘性流体仿真场景中的初始条件包括的温度插值到所述三维网格上；

通过求解以下方程，模拟热量传导过程，以确定经过时间步长后流体粒子在所述三维网络中的三维网格的温度：

其中，T表示温度，b表示热传导模型的热扩散系数，t表示时间，Δt表示时间步长，x表示在所述三维网格中的网格点的坐标的横坐标，y表示在三维网格中的网格点的坐标的纵坐标，Z表示在三维网格中的网格点的坐标的第三维坐标；

将所述经过时间步长后流体粒子在三维网络中的三维网格的温度与所述粘性流体仿真场景中的初始条件包括的温度的差确定温度改变量。

4.根据权利要求3所述的方法，其中，所述通过模拟热量传导过程，以确定经过时间步长后的粒子温度，还包括：

根据所述隐式流体粒子模型，将所述温度改变量插值回粒子；

通过以下公式，确定流体粒子的粒子温度：

TN＝αF十(1-α)P，

其中，TN表示粒子温度，α表示所述第一权重，且α取值范围为[0，1]，F表示根据所述FLIP方法求得的温度，P表示根据所述PIC方法求得的温度。

5.根据权利要求4所述的方法，其中，所述根据所述粒子温度，对所述粒子速度进行修正，包括：

通过以下公式，确定第一特征：

其中，D表示所述第一特征，u表示流体粒子在t时刻的速度，

表示梯度，

表示上述梯度

的转置；

根据所述第一特征，通过以下公式，确定摩擦应力：

其中，σ_f表示所述摩擦应力，p表示压力，D表示第一特征，|D|_F表示第一特征D的Frobenius范数；

响应于流体粒子的位置在场景内部，根据所述摩擦应力和所述粒子温度，通过以下公式，对流体粒子的粒子速度进行修正：

其中，u表示流体粒子的粒子速度，σ_f表示所述摩擦应力，

表示用中心差分法计算得到的滑动摩擦力σ_f的散度，β表示预设权重系数；

响应于流体粒子的位置在场景边界上，根据流体粒子的粒子速度，通过以下公式，修正流体粒子的切向速度：

其中，UT表示所述切向速度，μ表示摩擦系数，n表示法线，|u·n|表示法向速度的模。

Images