一种非牛顿流体模拟方法及装置

摘要

本发明提供一种非牛顿流体模拟方法及装置，该方法包括：基于连续性方程和动量方程，分别模拟非牛顿流体的液体相和悬浮颗粒相；通过流体阻力公式和流体压强公式，将液体相和悬浮颗粒相进行耦合，得到介于液体相和悬浮颗粒相之间性质的二相悬浮体模型；根据二相悬浮体模型的塑性剪切应力和体积变化，构建基于悬浮颗粒相的弹塑性本构模型；根据二相悬浮体模型的粒子骨骼结构变化，构建摩擦函数，并将摩擦函数添加到弹塑性本构模型中，以用于模拟非牛顿流体的流变学特性。本发明通过对其底层的真实物理机制的建模，实现非牛顿性质的模拟，提供物理上可解释的参数控制，从而统一地模拟不同种类的非牛顿流体，准确地刻画非牛顿流体的流变学特性。

说明书

技术领域

本发明涉及计算机图形学流体模拟与渲染技术领域，尤其涉及一种非牛顿流体模拟方法及装置。

背景技术

非牛顿流体不同于一般流体的特性在于，其粘滞性是随着剪切流的情况发生变化的，广泛存在于日常生活和工业应用当中，例如，番茄酱、奶油以及各种润滑剂。由于非牛顿流体用途广泛，且展现出丰富的视觉效果，图形学领域对其的关注度也与日俱增。相比一般的牛顿流体，非牛顿流体的复杂流变学特性还未被人们完全理解，其高度的非线性行为也对模拟提出了重大挑战。

目前业界普遍的做法是采用一种指数模型来描述非牛顿流体，例如，Herschel-Bulkley模型，

但是现有这些方法将非牛顿流体视为均匀的单一流体来处理，忽略了其内部的微观结构。但是绝大多数非牛顿流体属于悬浮体，即微小颗粒悬浮于牛顿流体形成的一种复合物系统，举例来说，各种酱汁、奶油、血液等，都属于这类悬浮物，并没有真实的模拟出非牛顿流体的特性。因此，现在亟需一种非牛顿流体模拟方法及装置来解决上述问题。

发明内容

针对现有技术存在的问题，本发明提供一种非牛顿流体模拟方法及装置。

本发明提供一种非牛顿流体模拟方法，包括：

基于连续性方程和动量方程，分别模拟非牛顿流体的液体相和悬浮颗粒相；

通过流体阻力公式和流体压强公式，将所述液体相和所述悬浮颗粒相进行耦合，得到介于液体相和悬浮颗粒相之间性质的二相悬浮体模型；

根据所述二相悬浮体模型的塑性剪切应力和体积变化，构建基于悬浮颗粒相的弹塑性本构模型；

根据所述二相悬浮体模型的粒子骨骼结构变化，构建摩擦函数，并将所述摩擦函数添加到所述弹塑性本构模型中，以用于模拟所述非牛顿流体的流变学特性。

根据本发明提供的一种非牛顿流体模拟方法，所述流体阻力公式为：

其中，k

所述流体压强公式为：

其中，p

根据本发明提供的一种非牛顿流体模拟方法，所述根据所述二相悬浮体模型的塑性剪切应力和体积变化，构建基于悬浮颗粒相的弹塑性本构模型，包括：

根据塑性剪切应力，构建第一屈服条件：

τ≤(μ+β)p；

其中，τ表示剪切应力，μ表示摩擦系数，I

根据塑性拉伸，构建第二屈服条件：

p+∈≥0；

其中，∈表示预设粘黏性阈值；

根据塑性压缩，构建第三屈服条件：

其中，φ表示悬浮颗粒相体积分数，g(φ)表示一个随φ变化的系数；

根据所述第一屈服条件、所述第二屈服条件和所述第三屈服条件各自对应的屈服表面，构建得到基于悬浮颗粒相的弹塑性本构模型。

根据本发明提供的一种非牛顿流体模拟方法，所述摩擦函数的公式为：

φ

其中，

根据本发明提供的一种非牛顿流体模拟方法，所述方法还包括：

通过物质点法，对所述弹塑性本构模型进行离散模拟求解，得到所述非牛顿流体的流变学特性。

根据本发明提供的一种非牛顿流体模拟方法，所述通过物质点法，对所述弹塑性本构模型进行离散模拟求解，所述方法还包括：

基于液体漂移速度模型，对液体和干燥颗粒之间的渗透现象进行模拟，所述液体漂移速度模型为：

其中，c

根据本发明提供的一种非牛顿流体模拟方法，所述基于连续性方程和动量方程，分别模拟非牛顿流体的液体相和悬浮颗粒相，包括：

基于固体连续性方程和固体动量方程，模拟所述非牛顿流体中的悬浮颗粒相，所述固体连续性方程为：

所述固体动量方程为：

其中，

基于液体连续性方程和液体动量方程，模拟所述非牛顿流体中的液体相，所述液体连续性方程为：

所述液体动量方程为：

其中，φ

本发明还提供一种非牛顿流体模拟装置，包括：

第一模拟模块，用于基于连续性方程和动量方程，分别模拟非牛顿流体的液体相和悬浮颗粒相；

耦合模块，用于通过流体阻力公式和流体压强公式，将所述液体相和所述悬浮颗粒相进行耦合，得到介于液体相和悬浮颗粒相之间性质的二相悬浮体模型；

模型构建模块，用于根据所述二相悬浮体模型的塑性剪切应力和体积变化，构建基于悬浮颗粒相的弹塑性本构模型；

第二模拟模块，用于根据所述二相悬浮体模型的粒子骨骼结构变化，构建摩擦函数，并将所述摩擦函数添加到所述弹塑性本构模型中，以用于模拟所述非牛顿流体的流变学特性。

本发明还提供一种电子设备，包括存储器、处理器及存储在存储器上并可在处理器上运行的计算机程序，所述处理器执行所述程序时实现如上述任一种所述非牛顿流体模拟方法的步骤。

本发明还提供一种非暂态计算机可读存储介质，其上存储有计算机程序，该计算机程序被处理器执行时实现如上述任一种所述非牛顿流体模拟方法的步骤。

本发明提供的非牛顿流体模拟方法及装置，提出一个更具扩展性的非牛顿流体模拟方法，通过对其底层的真实物理机制的建模，实现非牛顿性质的模拟，提供物理上可解释的参数控制，从而统一地模拟不同种类的非牛顿流体，准确地刻画非牛顿流体的流变学特性。

附图说明

为了更清楚地说明本发明或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图作一简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图。

图1为本发明提供的非牛顿流体模拟方法的流程示意图；

图2为本发明提供的屈服表面的示意图；

图3为本发明提供的三种剪切变稠流变机制的示意图；

图4为本发明提供的非牛顿流体模拟装置的结构示意图；

图5为本发明提供的电子设备的结构示意图。

具体实施方式

为使本发明的目的、技术方案和优点更加清楚，下面将结合本发明中的附图，对本发明中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有作出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

虽然绝大多数非牛顿流体属于悬浮体，但现有针对非牛顿流体的模拟方法，忽略了其内部的微观结构，将非牛顿流体视为均匀的单一流体来处理，导致模拟得到非牛顿流体并不能完整的描述真实特性。基于此，本发明将非牛顿流体用两相模拟，可以模拟非均匀的非牛顿流体，甚至可以直接从原料(即水和粉末)出发模拟非牛顿流体的形成，因此大大扩展了可以模拟的现象。特别是对于一类非连续的剪切变稠现象，存在于如欧不裂(oobleck)等物质中，这类物质在剪切时可以展现出类似刚体一样的性质，但在静置时又能变回流体，这与其内部的结构变化息息相关。

图1为本发明提供的非牛顿流体模拟方法的流程示意图，如图1所示，本发明提供了一种非牛顿流体模拟方法，包括：

步骤101，基于连续性方程和动量方程，分别模拟非牛顿流体的液体相和悬浮颗粒相。

在本发明中，通过非牛顿流体中液体和悬浮颗粒对应的连续性方程和动量方程，以模拟两个独立的连续介质相为目标，模拟得到非牛顿流体的液体相和悬浮颗粒相，其中，悬浮颗粒相是将非牛顿流体中悬浮的微小颗粒，通过连续性方程和动量方程进行描述。具体地，连续性方程用于更新粒子密度，而动力方程包括重力、粘滞力、压强和剪切应力等，用于更新粒子速度。

步骤102，通过流体阻力公式和流体压强公式，将所述液体相和所述悬浮颗粒相进行耦合，得到介于液体相和悬浮颗粒相之间性质的二相悬浮体模型。

在本发明中，将流体阻力和流体压强施加到两相上，产生固液强耦合的效应，展现出介于固体和液体之间的性质，从而得到介于液体相和悬浮颗粒相之间性质的二相悬浮体模型。

步骤103，根据所述二相悬浮体模型的塑性剪切应力和体积变化，构建基于悬浮颗粒相的弹塑性本构模型。

在本发明中，根据二相悬浮体模型的特性，通过设计3个塑性屈服条件，从而构建得到基于悬浮颗粒相的弹塑性本构模型。

步骤104，根据所述二相悬浮体模型的粒子骨骼结构变化，构建摩擦函数，并将所述摩擦函数添加到所述弹塑性本构模型中，以用于模拟所述非牛顿流体的流变学特性。

在本发明中，在分别模拟非牛顿流体的液体相和悬浮颗粒相时，首先预设适合目标非牛顿流体的模型参数，包括屈服条件参数，摩擦函数的更新参数，悬浮颗粒直径以及液体粘滞系数等，其中，摩擦函数用于描述微观的粒子骨骼结构在剪切流中发生的变化，从而将非连续剪切变稠现象模拟出来；进一步地，本发明采用物质点法，将待模拟的物质(即液体相和悬浮颗粒相)转化为离散的物质点粒子，将流体阻力和流体压强施加到两相上，产生固液强耦合的效应，每个粒子具有合适的密度、质量和体积分数，并为模拟区域生成适当尺寸和分辨率的网格，用于物质点法的计算；然后，通过弹塑性本构模型，计算空间中任意位置的两相各自的体积分数；在更新了悬浮颗粒的应力张量之后，将更新后的悬浮颗粒相粒子的应力张量，投射到弹塑性本构模型中3个屈服条件对应的屈服表面上，从而产生非牛顿流体的特性，模拟得到非牛顿流体的运动；最后，本发明将模拟时得到的每个时间步信息在显示终端上绘制，得到非牛顿流体的模拟图像，并导出为数据文件，用于后续处理。

本发明提供的非牛顿流体模拟方法，提出一个更具扩展性的非牛顿流体模拟方法，通过对其底层的真实物理机制的建模，实现非牛顿性质的模拟，提供物理上可解释的参数控制，从而统一地模拟不同种类的非牛顿流体，准确地刻画非牛顿流体的流变学特性。

在上述实施例的基础上，所述流体阻力公式为：

其中，k

所述流体压强公式为：

其中，p

在本发明中，将上述实施例模拟得到的两相，通过流体阻力和流体压强耦合，在流体阻力公式中，阻力与两相速度差成正比，与悬浮颗粒直径的平方成反比。在悬浮颗粒微小的情况下，阻力变得非常大，本发明采用反向欧拉方法计算。本发明将两个相通过阻力和压强耦合，从而模拟非牛顿流体处于流体和固体之间的复杂流变学特性。

在上述实施例的基础上，所述根据所述二相悬浮体模型的塑性剪切应力和体积变化，构建基于悬浮颗粒相的弹塑性本构模型，包括：

根据塑性剪切应力，构建第一屈服条件：

τ≤(μ+β)p；

其中，τ表示剪切应力；μ表示摩擦系数，I

根据塑性拉伸，构建第二屈服条件：

p+∈≥0；

其中，∈表示预设粘黏性阈值；

根据塑性压缩，构建第三屈服条件：

其中，φ表示悬浮颗粒相体积分数，g(φ)表示一个随φ变化的系数；

根据所述第一屈服条件、所述第二屈服条件和所述第三屈服条件各自对应的屈服表面，构建得到基于悬浮颗粒相的弹塑性本构模型。图2为本发明提供的屈服表面的示意图，三个屈服条件各自所对应的屈服表面(即f

本发明根据非牛顿流体应有的塑性剪切应力和体积变化，设计了悬浮颗粒相应服从的弹塑性本构模型，在这一模型中，悬浮颗粒相的内部应力受限于压强和摩擦系数的乘积，同时还受剪切率、悬浮颗粒相的体积分数等影响，从而构成了对悬浮体物理性质比较完整的描述。

在上述实施例的基础上，所述摩擦函数的公式为：

φ

其中，

在本发明中，采用Ψ

本发明采用体积分数描述空间中悬浮颗粒相和液体相的占比，体积分数对悬浮体的特性具有关键的作用，其中，悬浮颗粒相的体积分数记为φ，它描述了悬浮颗粒相占据空间的密集程度，根据定义，在同一位置液体相的体积分数为1-φ。在本发明中，悬浮体的等效粘滞度随φ的增加而增加，并在φ接近φ

在上述实施例的基础上，所述方法还包括：

通过物质点法，对所述弹塑性本构模型进行离散模拟求解，得到所述非牛顿流体的流变学特性。

在本发明中，采用物质点法进行模拟，具体步骤如下：

步骤201，将粒子的质量和速度插值到网格上，计算压强梯度和内应力；

步骤202，更新网格速度，应用边界条件；

步骤203，将更新后的网格速度插值回粒子，更新粒子的位置，并计算粒子所在位置的速度梯度；

步骤204，更新悬浮相粒子携带的应力张量，计算液体相粒子的压强；

步骤205，根据悬浮相的弹塑性屈服条件，获得剪切率

步骤206，重复步骤201至步骤205，直至算法结束。

在上述实施例的基础上，所述通过物质点法，对所述弹塑性本构模型进行离散模拟求解，所述方法还包括：

基于液体漂移速度模型，对液体和干燥颗粒之间的渗透现象进行模拟，所述液体漂移速度模型为：

其中，c

在本发明中，为实现模拟液体与干燥颗粒之间的渗透，本申请在将更新后的网格速度插值回粒子时(即上述实施例中的步骤203)，引入基于漂移速度模型的扩散，从而模拟两相之间的渗透扩散现象。

在上述实施例的基础上，所述基于连续性方程和动量方程，分别模拟非牛顿流体的液体相和悬浮颗粒相，包括：

基于固体连续性方程和固体动量方程，模拟所述非牛顿流体中的悬浮颗粒相，所述固体连续性方程为：

所述固体动量方程为：

其中，

基于液体连续性方程和液体动量方程，模拟所述非牛顿流体中的液体相，所述液体连续性方程为：

所述液体动量方程为：

其中，φ

本发明构建模拟非牛顿流体的统一框架，通过物理可解释的参数模拟多种非牛顿流体。在模拟任意体积分数的悬浮体时，可以容易地模拟它们之间的耦合。本发明所采用的屈服准则和摩擦函数提供了准确的物理工具，可模拟非连续的剪切变稠现象。最后，本发明包含的渗透扩散机制，可以模拟初始状态为分离的相之间的混合过程，从而模拟非牛顿流体的形成。

图4为本发明提供的非牛顿流体模拟装置的结构示意图，如图4所示，本发明提供了一种非牛顿流体模拟装置，包括第一模拟模块401、耦合模块402、模型构建模块403和第二模拟模块404，其中，第一模拟模块401用于基于连续性方程和动量方程，分别模拟非牛顿流体的液体相和悬浮颗粒相；耦合模块402用于通过流体阻力公式和流体压强公式，将所述液体相和所述悬浮颗粒相进行耦合，得到介于液体相和悬浮颗粒相之间性质的二相悬浮体模型；模型构建模块403用于根据所述二相悬浮体模型的塑性剪切应力和体积变化，构建基于悬浮颗粒相的弹塑性本构模型；第二模拟模块404用于根据所述二相悬浮体模型的粒子骨骼结构变化，构建摩擦函数，并将所述摩擦函数添加到所述弹塑性本构模型中，以用于模拟所述非牛顿流体的流变学特性。

本发明提供的非牛顿流体模拟装置，提出一个更具扩展性的非牛顿流体模拟方法，通过对其底层的真实物理机制的建模，实现非牛顿性质的模拟，提供物理上可解释的参数控制，从而统一地模拟不同种类的非牛顿流体，准确地刻画非牛顿流体的流变学特性。

本发明提供的装置是用于执行上述各方法实施例的，具体流程和详细内容请参照上述实施例，此处不再赘述。

图5为本发明提供的电子设备的结构示意图，如图5所示，该电子设备可以包括：处理器(processor)501、通信接口(CommunicationsInterface)502、存储器(memory)503和通信总线504，其中，处理器501，通信接口502，存储器503通过通信总线504完成相互间的通信。处理器501可以调用存储器503中的逻辑指令，以执行非牛顿流体模拟方法，该方法包括：基于连续性方程和动量方程，分别模拟非牛顿流体的液体相和悬浮颗粒相；通过流体阻力公式和流体压强公式，将所述液体相和所述悬浮颗粒相进行耦合，得到介于液体相和悬浮颗粒相之间性质的二相悬浮体模型；根据所述二相悬浮体模型的塑性剪切应力和体积变化，构建基于悬浮颗粒相的弹塑性本构模型；根据所述二相悬浮体模型的粒子骨骼结构变化，构建摩擦函数，并将所述摩擦函数添加到所述弹塑性本构模型中，以用于模拟所述非牛顿流体的流变学特性。

此外，上述的存储器503中的逻辑指令可以通过软件功能单元的形式实现并作为独立的产品销售或使用时，可以存储在一个计算机可读取存储介质中。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分或者该技术方案的部分可以以软件产品的形式体现出来，该计算机软件产品存储在一个存储介质中，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例所述方法的全部或部分步骤。而前述的存储介质包括：U盘、移动硬盘、只读存储器(ROM，Read-OnlyMemory)、随机存取存储器(RAM，RandomAccessMemory)、磁碟或者光盘等各种可以存储程序代码的介质。

另一方面，本发明还提供一种计算机程序产品，所述计算机程序产品包括存储在非暂态计算机可读存储介质上的计算机程序，所述计算机程序包括程序指令，当所述程序指令被计算机执行时，计算机能够执行上述各方法所提供的非牛顿流体模拟方法，该方法包括：基于连续性方程和动量方程，分别模拟非牛顿流体的液体相和悬浮颗粒相；通过流体阻力公式和流体压强公式，将所述液体相和所述悬浮颗粒相进行耦合，得到介于液体相和悬浮颗粒相之间性质的二相悬浮体模型；根据所述二相悬浮体模型的塑性剪切应力和体积变化，构建基于悬浮颗粒相的弹塑性本构模型；根据所述二相悬浮体模型的粒子骨骼结构变化，构建摩擦函数，并将所述摩擦函数添加到所述弹塑性本构模型中，以用于模拟所述非牛顿流体的流变学特性。

又一方面，本发明还提供一种非暂态计算机可读存储介质，其上存储有计算机程序，该计算机程序被处理器执行时实现以执行上述各实施例提供的非牛顿流体模拟方法，该方法包括：基于连续性方程和动量方程，分别模拟非牛顿流体的液体相和悬浮颗粒相；通过流体阻力公式和流体压强公式，将所述液体相和所述悬浮颗粒相进行耦合，得到介于液体相和悬浮颗粒相之间性质的二相悬浮体模型；根据所述二相悬浮体模型的塑性剪切应力和体积变化，构建基于悬浮颗粒相的弹塑性本构模型；根据所述二相悬浮体模型的粒子骨骼结构变化，构建摩擦函数，并将所述摩擦函数添加到所述弹塑性本构模型中，以用于模拟所述非牛顿流体的流变学特性。

以上所描述的装置实施例仅仅是示意性的，其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本实施例方案的目的。本领域普通技术人员在不付出创造性的劳动的情况下，即可以理解并实施。

通过以上的实施方式的描述，本领域的技术人员可以清楚地了解到各实施方式可借助软件加必需的通用硬件平台的方式来实现，当然也可以通过硬件。基于这样的理解，上述技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在计算机可读存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行各个实施例或者实施例的某些部分所述的方法。

最后应说明的是：以上实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。