CN110909472A - 一种基于混合模型的粉末材料仿真方法 - Google Patents

一种基于混合模型的粉末材料仿真方法 Download PDF

Info

Publication number
CN110909472A
CN110909472A CN201911181144.0A CN201911181144A CN110909472A CN 110909472 A CN110909472 A CN 110909472A CN 201911181144 A CN201911181144 A CN 201911181144A CN 110909472 A CN110909472 A CN 110909472A
Authority
CN
China
Prior art keywords
smoke
particles
particle
powder material
solid
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Pending
Application number
CN201911181144.0A
Other languages
English (en)
Inventor
高阳
郝爱民
徐英浩
李帅
秦洪
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Beihang University
Original Assignee
Beihang University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Beihang University filed Critical Beihang University
Priority to CN201911181144.0A priority Critical patent/CN110909472A/zh
Publication of CN110909472A publication Critical patent/CN110909472A/zh
Pending legal-status Critical Current

Links

Images

Classifications

    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T13/00Animation
    • G06T13/203D [Three Dimensional] animation
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T15/003D [Three Dimensional] image rendering
    • G06T15/005General purpose rendering architectures
    • GPHYSICS
    • G06COMPUTING; CALCULATING OR COUNTING
    • G06TIMAGE DATA PROCESSING OR GENERATION, IN GENERAL
    • G06T17/00Three dimensional [3D] modelling, e.g. data description of 3D objects
    • G06T17/20Finite element generation, e.g. wire-frame surface description, tesselation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Graphics (AREA)
  • Geometry (AREA)
  • Software Systems (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种基于混合模型的粉末材料仿真方法,涉及基于物理的流体仿真技术领域,以实现粉末材料的逼真建模为目的,将粉末材料分解为由颗粒状固体粒子和烟雾状烟雾两种材质构成,基于Affine Particle in Cell方法(APIC)和欧拉网格方法的混合模型,围绕粉末材料仿真展开研究。首先利用基于APIC方法的流体模型进行颗粒状粒子材质的建模,仿真粉末材料中颗粒固体的运动,然后利用Mohr‑Coulomb屈服准则和外摩擦条件修改APIC流体模型,实现颗粒状固体堆积、滑动的状态仿真;接着通过改进的欧拉方法对粉末材料的烟雾状态进行模拟,并通过提出的双向转化算法仿真粉末材料中固体与烟雾的相互转化,进而增强粉末状材料的真实感,实现对诸如粉末飘洒、烟雾爆炸等基于物理的粉末材料运动过程模拟。

Description

一种基于混合模型的粉末材料仿真方法
技术领域
本发明涉及基于计算机图形学的真实类流体材料仿真的技术领域,具体涉及一种基于混合模型的粉末材料仿真方法。
背景技术
近几十年来,关注实现流体各种性质模拟与仿真的工作和文章层出不穷,流体仿真的效果和效率也越来越好,但是其中一个分支——类流体及其衍生行为得到的关注远不如流体行为仿真的关注度。粉末材料的组成复杂,由微米级到毫米级不同尺度的微小粒子组成,其中毫米级的颗粒表现为类似沙粒的固体属性,而微米级的细小粒子则在运动过程中极易像烟雾一样飘散,同时在运动过程中,大颗粒固体粒子有可能破碎形成烟雾,烟雾也可能聚集生成粒子团。因此,基于计算机图形学的真实物理现象仿真方面,粉末材料的模拟复杂且困难。
在近年来的粉末或颗粒材质物体的仿真研究中,一种简单有效的模拟方法是把类似粉末的材质作为流体模拟仿真,而在真实的场景中,例如沙漠中一旦有风流动,地面的固体粒子中会有许多较小的颗粒在风力的影响下,变成浮动在空气中的烟尘。当风力达到一定程度时,空气中浮动的沙尘变得足够多,这时空气中的浮尘达到了一定的密度,沙尘暴就形成了。换言之,粉末状材质处在不同外力环境中时,会呈现出不同的状态,从固态(聚集成堆静止在地面)到液态(像流体一样在地面流动),甚至还有气态(在足够外力影响下浮到空中形成烟雾)。
为实现基于计算机图形学的真实感粉末材料仿真,可将粉末材料的状态简化为两种:粒子状态和烟雾状态。烟雾是指悬浮在空气中的固体微粒,根据国际标准化组织的规定,粒子直径小于75μm的固体悬浮物定义为烟雾,习惯上烟雾又可以称为灰尘、尘埃、烟尘、沙尘等。烟雾在生活中几乎到处可见,土壤和岩石风化后裂解成许多细小的颗粒。而粉末材料的粒子状态,则是具有固体颗粒材料(Granular Material)(通常指粒径大于100μm的固体颗粒)性质的状态。这两种状态的粒径差异对其运动性质产生了巨大影响,使它们具有截然不同的受力和性质:考虑到颗粒材料较小但不可忽视的粒子直径,当粉末材料处于固体颗粒状态时,它的个体颗粒间的相互作用主要取决于颗粒间的摩擦力,而空气对颗粒产生的作用和其相互作用则可以忽略不计;对于粉末材料的烟雾状态,由于其粒子直径较小,空气与颗粒之间的相互作用不但不可以忽略,还占据对其运动状态影响的主导地位,,因此颗粒的间距比较大,接触力和干摩擦与空气产生的影响相比较几乎可以忽略。因此当粉末材料处于烟雾状态时,其会呈现出雾状材料(Fog Material)的物理特性。
现有的流体和类流体模拟仿真技术具有一定的基础和研究,包括对流体材料(Fluid Material)、颗粒材料(Granular Material)和雾状材料(Fog Material)等的模拟仿真,但是没有专注于粉末材料(Powdered Material)的仿真研究,难以逼真模拟粉末材料的运动过程。而粉末材料运动又广泛应用于影视制作、游戏场景建模等方面。因此针对市场和应用的大量需求,考虑到当前技术缺失的现状,对于粉末材料的研究分析与建模就具有着重要的意义和价值。
发明内容
本发明要解决的技术问题是:针对真实感流体环境建模中粉末材料难以逼真模拟的问题,提出一种基于APIC方法和欧拉方法结合的混合建模方法,设计并实现一个将粉末材料分为可相互转换的颗粒材质和烟雾材质的仿真模型,实现粉末材料的逼真建模和真实感仿真(如烟尘弥漫、粒子团洒落等)。
本发明采用的技术方案为:一种基于混合模型的粉末材料仿真方法,包括以下四个步骤:
步骤(1)、根据粉末状材料由颗粒不同大小的固体粒子组成的特性,将粉末状材料分为固体粒子和烟雾两种材质分别进行建模,以增加仿真过程中的真实感。其中,改进基于混合模型表示的仿射单元格粒子方法(Affine Particle In Cell,APIC)用来计算固体粒子运动,基于网格表示的欧拉方法计算烟雾运动;
步骤(2)、改进APIC流体仿真模型,引入摩尔-科伦坡(Mohr-Coulomb)屈服条件,计算固体粒子所受的应力,根据固体粒子间受力,判断固体粒子集合之间受力情况:保持满足屈服条件的受静摩擦力的粒子间相对静止,同时保证固体粒子集合具有统一的线速度和角速度;对不满足屈服条件的固体粒子作为滑动粒子,使用中心差分法计算滑动速度,同时执行摩擦边界条件,根据滑动粒子的切向速度与法向速度的比值,给切向速度施加负增量。通过上述手段将APIC流体模型进行改进,使之可求解颗粒状固体的静止、滑动过程,用于模拟粉末材料中的固体粒子材质;
步骤(3)、在进行粉末材料的烟雾粒子仿真时,将重力与空气阻力作为动力学方程的外力项,然后根据APIC模型计算烟雾粒子速度场,并通过求解Navier-Stokes(N-S)方程以获得基于网格的密度场,再利用半拉格朗日法更新密度场,最终建立一个可模拟烟雾粒子运动的烟雾求解器;
步骤(4)、设计一种双向相态转化机制,并在固体粒子仿真过程的每个时间步内监测固体粒子的速度和位置变化,设置转化阈值,当固体粒子速度超过阈值时,固体粒子在质量守恒条件的约束下转化成为烟雾。同时,烟雾的运动速度小于阈值时,在质量守恒的前提下转换回固体粒子,从而实现烟雾与固体粒子的双向转换。
本发明的原理在于:
本发明提出了一种基于混合建模的粉末材料仿真建模方法,用于进行诸如粉末飘洒现象的真实感仿真。将粉末材料的模拟简化为两种状态:基于粒子表示的固体颗粒状态和基于密度场表示的烟雾状态。首先对基于APIC方法的流体求解器进行了基于GPU的算法实现,然后根据Mohr-Coulomb屈服准则将流体求解器改进为可模拟粒子运动及堆积、滑动过程的颗粒材质求解器,用于模拟粉末材料的固体颗粒状态;之后,本发明修改了基于网格方法的欧拉模型,用于模拟粉末材料的烟雾状态;最后,为了实现两种材质状态的相互转化,本发明还提出了一种双向转化机制。
本发明与现有技术相比的优点在于:
(1)本发明提出的基于混合模型的粉末材料仿真方法,应用于计算机动画和计算机可视化仿真领域,相对于现有的类流体仿真方法物理更加真实,算法更加严谨、精确。
(2)本发明的方法设计了粉末材料的两种材质,用固体粒子和烟雾这两种材质共同表示粉末材料,相较于现有的单一表示的仿真方法模型,在满足物理真实的同时,能够更好地表现粉末材料的临界特性,提升了仿真效果的真实感。
(3)相较于现有计算流体动力学领域的单一表示的类流体仿真方法,本发明的方法设计了粉末材料粒子和烟雾两种材质之间的双向转化算法,保证了提出方法的物理正确性和真实感效果,保证了仿真中物质和能量转化时的物质守恒条件,进而可对诸如烟尘弥漫、烟雾爆炸等场景进行逼真模拟。
附图说明
图1为本发明提出的混合模型建模框架;
图2为本发明提出的仿真方法运算示意图;
图3为本发明提出的粉末材料仿真方法运算步骤示意图;
图4为本发明提出的混合模型三维仿真结果;
图5为本发明提出的混合模型三维仿真渲染结果。
具体实施方式
图1给出了本发明提出的混合模型的建模框架,图2是本发明提出的粉末材料仿真方法的运算示意图,图3是本发明提出的粉末材料仿真方法的总体运算步骤示意图,下面结合其他附图及具体实施方式进一步说明本发明。
如图1所示,本发明提出的方法框架主要包含两个可交互的模型:固体粒子仿真模型和烟雾仿真模型,并且在耦合算法的控制下,粒子可以汽化生成烟雾,烟雾也能够沉淀重新生成粒子。
如图2所示,本发明提出一种基于混合模型的粉末材料仿真方法,具体实施为耦合APIC流体模型和欧拉网格模型进行共同仿真,主要步骤介绍如下:
1、粉末材料固体粒子状态仿真
使用传统的欧拉方法求解N-S方程,对粒子进行插值时使用速度而不是速度改变量。在每个时间步长中将粒子属性传递到网格中,进行动力学运算,使其具有固体粒子的特性,然后再将网格属性插值回粒子,并通过额外携带的矩阵运算,在数值传递的同时保证角动量守恒,以满足旋转、剪切等运动的要求。该方法中使用交错网格(MAC-Grid)来进行数值存储,使用三个方向上的3×1向量,分别代表三维中仿射速度场的三个矢量分量存储在每个粒子上。因此在进行粒子到网格的转换时,使用如下公式:
Figure BDA0002291290550000041
Figure BDA0002291290550000042
其中上标n表示一次仿真中的第n个时间步;下标a代表x,y或z方向的网格面,ea表示a方向的单位向量;下标p和下标i分别代表粒子p和网格面节点i处的数值。此外,
Figure BDA0002291290550000043
表示插值权重,通常N(x)为三线性插值的核函数;xai是与方向a相关联的MAC网格面所在的位置,xp是粒子所在的位置;
Figure BDA0002291290550000044
Figure BDA0002291290550000045
分别是MAC网格面上存储的质量和速度分量,
Figure BDA0002291290550000046
是粒子携带的仿射分量。
为了模拟粉末材料的固体粒子状态,在动力学的压强投影步骤之后增加一个用于对固体粒子摩擦力和塑性进行控制的部分,执行如下步骤:
首先,使用标准的中心差分法(Standard Central Differences)来估算网格中每个单元格的3×3张量应变率D,使用如下公式:
Figure BDA0002291290550000047
其中
Figure BDA0002291290550000048
表示梯度算子,
Figure BDA0002291290550000049
是网格上速度场的梯度,而
Figure BDA00022912905500000410
代表它的转置。然后对使用如下公式对固体粒子流动时的摩擦应力σf进行计算:
Figure BDA00022912905500000411
在公式(4)中,
Figure BDA00022912905500000412
是摩擦角度,代表固体粒子材料静止时沙堆的最大坡度,
Figure BDA00022912905500000413
越小意味着静止的沙堆状态越扁平;p是压强投影步骤中计算得到的压力,其梯度表示对速度的影响,也就意味着固体粒子内部的流动方向;|D|F是应变率张量D的Frobenius范数。同时,通过如下公式计算σrigid
Figure BDA0002291290550000051
通过公式(5)计算得到的σrigid的意义是,在两个单元格之间,其表示的固体粒子相对凝固(即非塑性)所需的最小应力,其中ρ是该单元格的密度,Δx是单元格的边长,Δt是每个循环进行的时间步长。通过公式(4)和(5)得到了固体粒子滑动的摩擦力分量σf和固体粒子静止的摩擦力分量σrigid,也可以将其理解为对于每个单元格,在其中心处的摩擦力法向分量σrigid和切向分量σf(沿着压强投影的梯度方向)。通过比较摩擦力在这两个方向上的分量大小,就可以获得每个单元格用于判断的屈服准则:
σfrigid+c (6)
公式(6)中的c是凝固系数(Cohesion Coefficient),是对于基础Mohr-Coulumb屈服准则的一个常用的增强项。凝聚系数c适用于沙土和其它支持屈服前张力的微粘性材质,即使是对于固体粒子这种无粘性材料,给予c一个微小的值也能够避免数值误差引发的不稳定小滑坡现象,而对于固体粒子的滑动不产生视觉可见的影响。
随后,根据公式(6)的判断结果决定如何对单元格进行处理。如果单元格满足屈服准则,则将该单元格标记为刚性,并在该单元格中存储σrigid用于后续计算;否则,存储滑动摩擦力σf。然后对两个区域进行不同的处理:对于刚性单元格,在全部的模拟空间中找到标记为刚性的单元格的所有联通分量组,并将每个联通组中的速度场投影到刚体运动的空间:
Figure BDA0002291290550000052
公式(7)中
Figure BDA0002291290550000053
表示刚性联通分量组的质心,
Figure BDA0002291290550000054
表示组内的平均线速度,
Figure BDA0002291290550000055
表示组内的平均角速度。这个过程称为相对固化(Rigidification)操作,实质就是保证同一个联通分量组内的格子的角速度和速度保持一致,这个过程的计算通过动量守恒和角动量守恒来完成。
对于所有非刚性单元格(即不满足屈服准则的格子),使用如下公式对速度进行更新:
Figure BDA0002291290550000056
公式(8)中的σf正是之前存储在单元格的滑动摩擦力σf
Figure BDA0002291290550000057
是用中心差分法计算的σf的散度。由于u是存储在交错网格节点的速度场,所以摩擦力的散度可以很方便地通过中心差分法应用到速度场上,避免散度计算的空间零问题。
网格到粒子的转换使用如下公式:
Figure BDA0002291290550000061
Figure BDA0002291290550000062
其中上标n+1表示这个时间步结束后的下一个时间步中使用的存储物理量;
Figure BDA0002291290550000063
表示与之前相同的插值权重的梯度;
Figure BDA0002291290550000064
表示网格上的
Figure BDA0002291290550000065
在经过动力学运算后得到的新速度值。通过这个公式将权重梯度和速度的乘积求和,得到的就是这个粒子由其周围网格节点(三维网格中通常是8个)速度共同组成和影响的仿射状态向量。
2、粉末材料烟雾状态仿真
粉末材料烟雾状态建模过程,使用基于网格表示的欧拉求解器进行体积建模,同时改良了控制方程中的外力项和场源,求解器中的烟雾速度场ud的控制方程如下:
Figure BDA0002291290550000066
这个控制方程由基础的Navier-Stokes方程变形衍生而来,其中ud为烟雾的速度场。对于烟雾速度场ud的求解,Navier-Stokes方程中的外力项被阻力项代替,求解方式与普通流体求解器相同。通常的雾状材料建模仿真中,会认为重力和浮力相互抵消而忽略,而本方法为实现守恒与最终的稳态,以受到阻力的重力作为外力项。对于满足转化条件的固体粒子状态粒子,将其转化为烟雾(dust)状态,即删除该粒子,并在该粒子所处位置产生烟雾(dust)的密度场和速度场的值,并将其叠加进已有的场。密度场由如下方程决定:
Figure BDA0002291290550000067
其中,μ表示扩散系数,d表示耗散系数,S表示产生密度场的源,即固体粒子状态的粒子转化而来的密度场。扩散系数μ越高,代表当固体粒子转化成烟雾时,会对更大范围的场产生影响;由于粉末材料需要满足质量守恒,而不会像普通雾状材料一样耗散在空气中,本发明的耗散系数d通常设置为0,即在这个方程中忽略烟雾的耗散。
3、粉末材料两种状态的双向转化
在求解粉末材料运动函数的过程中,存在粉末材料两种状态之间的相互转化,包括固体粒子状态向烟雾状态的转化和烟雾状态向固体粒子状态的转化。这种转化实现物质和能量转移,满足物质和能量的守恒,使用密度和速度值作为阈值设置转化条件。
首先是固体粒子状态向烟雾状态的转化,即粒子表示转化成网格场值的表示。经过分析后选定的转化条件是粒子速度与其周围一定半径内粒子的平均速度的差值,即粒子与周围的相对速度,计算公式如下:
Figure BDA0002291290550000071
其中S是粒子p相邻单元格中的周围粒子,ωi=N(dij,h)是一个平滑的加权核函数,其中dij是第i个粒子和第j个粒子之间的距离,h是核函数的半径,控制该粒子受影响的范围。通过该控制方程计算得到的
Figure BDA0002291290550000072
就是粒子p与其周围半径h内所有粒子的平均速度的差值,即该粒子的相对速度,当
Figure BDA0002291290550000073
满足本发明设置的相对速度阈值
Figure BDA0002291290550000074
时,将其标记为待转化粒子并在之后一起进行转化过程。同时,当粒子绝对速度达到阈值
Figure BDA0002291290550000075
时,也发生转化。另一方面,烟雾状态向固体粒子状态的转化实质上是场的表示转化成粒子表示。在类似于静电吸引力的驱动下,烟雾密集区域将发生烟雾状态向固体粒子状态的转化,生成一个或多个固体粒子状态的粒子,即当一个网格单元的密度达到阈值
Figure BDA0002291290550000076
此时,该单元格中的烟雾聚集从而转化成固体粒子。其中,N是固体粒子转化成烟雾状态时受其影响的单元格数量(8或27);
Figure BDA0002291290550000077
是单个粒子发生转换时每个单元格受其影响增加的密度;α是一个转化系数,用于控制转化阈值的大小,代表α个粒子全部转化到一个网格获得的密度值,选择设置α=1方便质量守恒和生成粒子。在密度达到ρthres的前提下,烟雾单元格速度还需要小于速度阈值
Figure BDA0002291290550000078
转化才会发生。同时还设置一个静止速度阈值
Figure BDA0002291290550000079
作为转化条件,当场中单元格的速度低于最小速度阈值
Figure BDA00022912905500000710
时,则这个格子的烟雾停止扩散,将其转化为固体粒子状态并生成固体粒子,然后在场中减掉相应的数值。在特殊情况下,如果某些单元格的速度低于最小速度阈值
Figure BDA00022912905500000711
但其密度也特别低,不足以生成一个(甚至半个)固体粒子,则将其视为肉眼无法观察到离散微量烟雾,这种微量的耗散是合理的。这里使用的密度下限被称为耗散密度ρdis,其意义是密度低于耗散密度ρdis的单元格被认为是烟雾含量可忽略的耗散单元格。耗散密度ρdis应该被设置为极低,通常设置的耗散密度
Figure BDA00022912905500000712
Figure BDA00022912905500000713
简言之,发生烟雾状态向固体粒子状态转化的条件有两个:一是单元格的密度大于ρthres且单元格速度小于
Figure BDA00022912905500000714
即减缓速度的密集烟雾;二是单元格的速度小于
Figure BDA00022912905500000715
且单元格密度大于耗散密度ρdis,即不可忽视的低速烟雾。
由此,以上过程建立起完整的粉末材料两种状态之间相互转化的物理机制。
4、求解粉末材料运动方程步骤
求解粉末材料运动方程的步骤如图3所示,在一个时间步长Δt内,执行如下具体运算步骤如下:
首先对粒子进行迁移(advect)操作,根据粒子速度移动粒子的位置;然后对粒子施加外力(重力、风场风力等);然后将粒子属性通过仿射变换映射到网格上;执行摩擦边界条件,根据边界单元的切向速度与法向速度的比值,来给边界切向速度一个负增量;计算每个单元网格的应变率张量,并将结果代入预设的屈服准则进行判断,对于满足屈服准则的网格单元执行相对凝固处理,对不满足条件的单元更新其速度;将网格属性通过仿射变换映射回粒子,对粒子进行属性的更新;
然后对预设的标准进行计算和判断,根据结果执行固体粒子状态向烟雾状态的转化,在转化发生的位置生成烟雾场的值;
之后是基于欧拉求解器的烟雾模拟部分,计算密度场,包括之前步骤中根据质量和动量守恒转化产生的密度;然后通过计算获得其速度场,并更新密度场和速度场的值;
对预设的标准进行计算和判断,然后根据结果执行烟雾状态向固体粒子状态的转化,在转化发生的位置生成具有相同物理属性的粒子;
最后,时间步长Δt+1,进入下一个仿真的时间步内,重复以上步骤,实现包含两种状态的粉末材料的连续仿真。
如图4和图5所示,为证明本发明在计算机动画和计算机仿真领域的正确性,设计了一个封闭边界的三维流场场景,中间正中心有一个固定的实体模型,场景上方偏右有一个类似水龙头的材料出口作为粉末材料的输入,不停倾倒粉末在场景中间。图4中兔子模型表示不动的固体模型粒子,自上而下的离散粒子表示处于固体粒子状态的粉末材料,烟雾的表示处于烟雾状态的粉末材料密度场。通过这两个图可以明显看出,本发明能够真实有效的用两种状态的共同表示对处于运动临界状态的粉末材料进行仿真。图5是本发明涉及方法的渲染效果。

Claims (4)

1.一种基于混合模型的粉末材料仿真方法,其特征在于包括以下步骤:
步骤(1)、根据粉末状材料由颗粒不同大小的固体粒子组成的特性,将粉末状材料分为固体粒子和烟雾两种材质分别进行表示,以增加仿真过程中的真实感,固体粒子的运动过程采用基于混合模型表示的仿射单元格粒子模型(Affine Particle In Cell,APIC)进行仿真,烟雾粒子的运动过程采用基于网格表示的欧拉模型进行仿真;
步骤(2)、基于步骤(1)中混合模型表示的APIC模型,引入摩尔-科伦坡(Mohr-Coulomb)屈服条件,计算固体粒子所受的应力,根据计算得到的固体粒子间应力,判断固体粒子集合之间受力情况:保持满足屈服条件的受静摩擦力的粒子间相对静止,同时保证固体粒子集合具有统一的线速度和角速度;对不满足屈服条件的固体粒子作为滑动粒子,使用中心差分法计算滑动速度,同时执行摩擦边界条件,根据滑动粒子的切向速度与法向速度的比值,给切向速度施加负增量;通过将APIC流体模型进行改进,使之能够计算颗粒状固体粒子的运动,从而模拟粉末材料中的固体粒子的运动;
步骤(3)、基于步骤(1)中网格表示的欧拉模型,在进行粉末材料的烟雾粒子仿真时,将重力与空气阻力作为动力学方程的外力项,然后根据APIC模型计算烟雾粒子速度场,并通过求解Navier-Stokes(N-S)方程以获得基于网格的密度场,再利用半拉格朗日法更新密度场,最终建立一个可模拟烟雾粒子运动的烟雾求解器;
步骤(4)、设计一种双向相态转化机制,并在固体粒子运动的每个时间步内监测固体粒子的速度和位置变化,设置转化阈值,当固体粒子速度超过所述阈值时,固体粒子在质量守恒条件的约束下转化成为烟雾;同时,烟雾的运动速度小于或等于阈值时,在质量守恒的前提下转换回固体粒子,实现烟雾粒子与固体粒子的双向转换,最终达到对烟雾与固体粒子混合的粉末状材料仿真。
2.根据权利要求1所述的基于混合模型的粉末材料仿真方法,其特征在于:步骤(2)中,所述用APIC方法仿真颗粒状固体的运动,是在APIC的压强投影步骤之后增加一个用于对固体粒子摩擦力和塑性进行控制的部分,执行如下步骤:
首先,使用标准的中心差分法来估算网格中每个单元格的3×3张量应变率D,使用如下公式:
Figure FDA0002291290540000011
其中
Figure FDA0002291290540000021
是网格上速度u的梯度,而
Figure FDA0002291290540000022
代表它的转置,然后对使用如下公式对固体粒子流动时的摩擦应力σf进行计算:
Figure FDA0002291290540000023
在以上公式中,
Figure FDA0002291290540000024
是摩擦角度,代表固体粒子材料静止时沙堆的最大坡度,
Figure FDA0002291290540000025
越小意味着静止的沙堆状态越扁平,p是压强投影步骤中计算得到的压力,其梯度表示对速度的影响,也就意味着固体粒子内部的流动方向;|D|F是应变率张量D的Frobenius范数。
3.根据权利要求1所述的基于混合模型的粉末材料仿真方法,其特征在于:步骤(3)中,在进行粉末材料的烟雾粒子仿真时,使用基于网格表示的欧拉求解器进行体积建模,烟雾速度场ud的控制方程如下:
Figure FDA0002291290540000026
这个控制方程由基础的Navier-Stokes方程变形衍生而来,其中ud为烟雾的速度场,
Figure FDA0002291290540000027
为梯度算子,ρ是烟雾密度,p是压强相,f是外力项,对于满足转化条件的固体粒子,将其转化为烟雾状态后即删除该粒子,并在该粒子所处位置产生烟雾的密度场和速度场的值,烟雾体积与密度的乘积需要和转化前的粒子质量相等以保证质量守恒,密度场由如下方程决定:
Figure FDA0002291290540000028
其中,μ表示扩散系数,d表示耗散系数,S表示产生密度场的源,即固体粒子状态的粒子转化而来的密度场。
4.根据权利要求1所述的基于混合模型的粉末材料仿真方法,其特征在于:在步骤(4)中,双向相态转化机制是固体粒子向烟雾状态的转化和烟雾向固体粒子状态的转化;
固体状态向烟雾状态的转化,即粒子属性转化成密度场,转化条件是粒子速度与其周围一定半径内粒子的平均速度的差值,即当前粒子与周围粒子的相对速度,计算公式如下:
Figure FDA0002291290540000029
其中S是粒子p相邻单元格中的周围粒子,ωi=N(dij,h)是一个平滑的加权核函数,dij是第i个粒子和第j个粒子之间的距离,h是核函数的半径,控制该粒子受影响的范围,当
Figure FDA00022912905400000210
满足相对速度阈值
Figure FDA00022912905400000211
时,该粒子将转化为烟雾状态;同时,当粒子绝对速度达到阈值
Figure FDA00022912905400000212
时,也发生转化,即发生固体粒子向烟雾状态的转化要满足两个条件之一:一是当前粒子与周围粒子或空气发生剧烈摩擦产生较大的相对速度;二是该粒子的绝对速度达到绝对速度阈值
Figure FDA00022912905400000213
另一方面,烟雾向固体粒子的转化是基于场的表示转化成粒子表示,当一个网格单元的密度达到阈值ρthres的时候,该单元格中的烟雾聚集从而转化成固体粒子,ρthres表示如下:
Figure FDA0002291290540000031
其中,N是固体粒子粒子转化成烟雾状态时受其影响的单元格数量;
Figure FDA0002291290540000032
是单个粒子发生转换时每个单元格受其影响增加的密度;α是一个转化系数,用于控制转化阈值的大小,代表α个粒子全部转化到一个网格单元格获得的密度值,选择设置α=1保证质量守恒条件;在密度达到ρthres的前提下,烟雾网格速度还需要小于速度阈值
Figure FDA0002291290540000033
转化才会发生,同时设置一个静止速度阈值
Figure FDA0002291290540000034
作为转化条件,当网格的速度低于最小速度阈值
Figure FDA0002291290540000035
时,则认为这个网格中的烟雾停止扩散,将其转化为固体粒子状态并生成固体粒子,然后在场中减掉相应的密度值;
由此,建立起完整的粉末材料两种状态之间相互转化的物理机制。
CN201911181144.0A 2019-11-27 2019-11-27 一种基于混合模型的粉末材料仿真方法 Pending CN110909472A (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201911181144.0A CN110909472A (zh) 2019-11-27 2019-11-27 一种基于混合模型的粉末材料仿真方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201911181144.0A CN110909472A (zh) 2019-11-27 2019-11-27 一种基于混合模型的粉末材料仿真方法

Publications (1)

Publication Number Publication Date
CN110909472A true CN110909472A (zh) 2020-03-24

Family

ID=69819831

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201911181144.0A Pending CN110909472A (zh) 2019-11-27 2019-11-27 一种基于混合模型的粉末材料仿真方法

Country Status (1)

Country Link
CN (1) CN110909472A (zh)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102147932A (zh) * 2011-03-30 2011-08-10 北京航空航天大学 一种基于可移动欧拉网格的模型驱动的烟雾模拟方法
CN102831275A (zh) * 2012-08-31 2012-12-19 中国科学技术大学 一种3d流体的仿真方法及系统
CN103699715A (zh) * 2013-12-01 2014-04-02 北京航空航天大学 一种基于光滑粒子流体动力学和非线性有限元的流固耦合方法
CN104268943A (zh) * 2014-09-28 2015-01-07 北京航空航天大学 一种基于欧拉-拉格朗日耦合方法的流体仿真方法
CN105389839A (zh) * 2015-11-06 2016-03-09 北京航空航天大学 基于流体分析的流体参数估计方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102147932A (zh) * 2011-03-30 2011-08-10 北京航空航天大学 一种基于可移动欧拉网格的模型驱动的烟雾模拟方法
CN102831275A (zh) * 2012-08-31 2012-12-19 中国科学技术大学 一种3d流体的仿真方法及系统
CN103699715A (zh) * 2013-12-01 2014-04-02 北京航空航天大学 一种基于光滑粒子流体动力学和非线性有限元的流固耦合方法
CN104268943A (zh) * 2014-09-28 2015-01-07 北京航空航天大学 一种基于欧拉-拉格朗日耦合方法的流体仿真方法
CN105389839A (zh) * 2015-11-06 2016-03-09 北京航空航天大学 基于流体分析的流体参数估计方法

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
G. I. KELBALIYEV 等: ""Modeling the granulation of powdered materials by rolling"", 《THEORETICAL FOUNDATIONS OF CHEMICAL ENGINEERING》 *
YANG GAO 等: ""A Hybrid Method for Powdered Materials Modeling"", 《25TH ACM SYMPOSIUM ON VIRTUAL REALITY SOFTWARE AND TECHNOLOGY》 *
闫澍旺 等: ""桶形基础液压下沉过程的耦合欧拉-拉格朗日有限元法分析"", 《岩士力学》 *

Similar Documents

Publication Publication Date Title
Macklin et al. Unified particle physics for real-time applications
Guendelman et al. Coupling water and smoke to thin deformable and rigid shells
US7647214B2 (en) Method for simulating stable but non-dissipative water
Feldman et al. Animating gases with hybrid meshes
Génevaux et al. Simulating Fluid-Solid Interaction.
CN109657322B (zh) 一种固液多相适用于泥石流的动力学数值模拟方法
CN108269299A (zh) 一种基于sph方法近似求解的粘性流体建模方法
CN108491619A (zh) 基于物理与非物理混合的复杂场景流固耦合高效模拟方法
CN109992830B (zh) 基于物质点方法的山体滑坡灾害场景模拟方法
Zhu et al. Animating Sand as a Surface Flow.
Chiba et al. Visual simulation of water currents using a particle‐based behavioural model
Wang et al. Real-time snowing simulation
CN107798198B (zh) 一种基于物理的融化现象逼真模拟方法
Nixon et al. A fluid-based soft-object model
Habte et al. Particle sedimentation using hybrid Lattice Boltzmann-immersed boundary method scheme
Chen et al. A heuristic approach to the simulation of water drops and flows on glass panes
Chang et al. A particle-based method for viscoelastic fluids animation
CN110909472A (zh) 一种基于混合模型的粉末材料仿真方法
Im et al. Visual simulation of rapidly freezing water based on crystallization
Chang et al. A particle-based method for granular flow simulation
Wong et al. Hybrid‐based snow simulation and snow rendering with shell textures
Chen et al. A hybrid method for water droplet simulation
Fujisawa et al. Particle-based shallow water simulation with splashes and breaking waves
Yang et al. Interactive coupling between a tree and raindrops
Haumann et al. An application of motion design and control for physically-based animation

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
RJ01 Rejection of invention patent application after publication
RJ01 Rejection of invention patent application after publication

Application publication date: 20200324