一种三维流固单向耦合的实现方法及系统的制作方法
【专利摘要】本发明公开了一种三维流固单向耦合的实现方法及系统，其中，该方法包括：将固体边界的三角形网格离散成粒子，并进行边界粒子均匀化采样；计算该固体中流体粒子的状态，包括：流体与固体边界处流体粒子的密度、流体粒子的压力及粘力，以及流体粒子表面张力；并利用一数值积分模式对流体粒子的速度和位置信息进行更新；当所述流体粒子穿透所述固体边界时，对所述流体粒子的速度和位置进行校正，从而实现三维流固单向耦合。本发明公开的方法及系统增强了边界处理时的稳定性，成功模拟出压力、粘力和阻止穿透效果的同时大大减小了计算的复杂度，并且有效模拟出表面张力的宏观效果，而且并不会出现直接模拟聚合力情况下的网状失真现象。
【专利说明】一种三维流固单向耦合的实现方法及系统
【技术领域】
[0001] 本发明涉及基于物理的计算机实时动画【技术领域】，尤其涉及一种三维流固单向耦 合的实现方法及系统
【背景技术】
[0002] 拉格朗日流体模拟在计算机动画领域一直是个热门的研究议题，其基于粒子的流 体表示天然地支持模拟一些微小尺度的流体现象，同时无条件地保证质量守恒。光滑粒 子动力学SPH(SmoothedParticleHydrodynamics)方法是一种典型的拉格朗日方法，由 Lucy、Monaghan&Gingold等人于1977年提出用来解决天体物理学中的行星运动问题，后来 被引入到计算流体力学领域。在模拟流体时它是将连续的流体用相互作用的质点组描述， 各个质点上承载各种物理量，包括质量、速度、密度等，通过求解质点组的动力学偏微分方 程组和跟踪每个质点的运动，求得整个系统的力学行为。SPH方法使用核函数对密度、压力、 粘力项进行离散化，得到Navier-Stokes方程组的离散化计算形式，从而在每步迭代过程 近似解出各个物理量，模拟流体运动。
[0003] 关于SPH方法模拟流固耦合现象属于流固两相交界处边界处理的范畴，由于SPH 方法的粒子属性导致其在处理多相流问题时需要花费的计算代价很大，目前不存在一个统 一的模型和方法能完美的解决固液耦合问题，几种具有代表性的方法如下：
[0004]I)Muller等人于2004年提出的基于Lennard-Jones势能的固液交互模型，其为流 体粒子和固体粒子间施加一个距离相关的力，该方法在宏观上可以近似模拟，但在微观尺 度上不满足流体运动的物理规律，而且计算所需的时间步长较大。
[0005] 2)SeungtaikOh等人于2009年提出的基于冲量边界力（IBF)方法，其为固液交互 所施加的边界力的计算基于粒子间碰撞过程中产生的冲量大小和方向，该方法需要计算虚 拟的冲量边界力。
[0006] 3)Becker等人于2009年引入directforce来为边界力建模，该方法采用预测校 正模式对与边界碰撞后粒子的位置和速度进行校正，达到没有穿透现象的发生并模拟了多 种边界条件，该方法对于流固双向耦合计算的复杂高，每个迭代中需要两次周围邻居粒子 查询。
【发明内容】
[0007] 本发明的目的是提供一种三维流固单向耦合的实现方法及系统，增强了边界处理 时的稳定性，成功模拟出压力、粘力和阻止穿透效果的同时大大减小了计算的复杂度，并且 有效模拟出表面张力的宏观效果，而且并不会出现直接模拟聚合力情况下的网状失真现 象。
[0008] 本发明的目的是通过以下技术方案实现的：
[0009] -种三维流固单向耦合的实现方法，该方法包括：
[0010] 将固体边界的三角形网格离散成粒子，并进行边界粒子均匀化采样；
[0011] 计算该固体中流体粒子的状态，包括：流体与固体边界处流体粒子的密度、流体粒 子的压力及粘力，以及流体粒子表面张力；
[0012] 并利用一数值积分模式对流体粒子的速度和位置信息进行更新；
[0013] 当所述流体粒子穿透所述固体边界时，对所述流体粒子的速度和位置进行校正， 从而实现三维流固单向I禹合。
[0014] 进一步的，所述将固体边界的三角形网格离散成粒子，并进行边界粒子均匀化采 样包括：
[0015] 读取固体模型的三角形网格，并建立固体模型的周围符号距离场，具体的：将固体 模型中每个三角形网格建立一个包围盒，并为每个包围盒以统一的方式划分网格；计算每 个网格点到固体模型的符号距离，符号距离场以一个哈希表来存储，形成一个存在于固体 模型周围的窄带宽的符号距离场；
[0016] 抽取符号距离场的等值面，并在等值面的每个三角形网格上初始化边界粒子，每 个三角形网格上布置gA/Jir2个粒子；其中，g表示控制粒子密度的参数，A表示每个三角形 网格的面积，r表示粒子的半径；
[0017] 根据粒子的分布来决定每个粒子在固体表面法向和切向^的速度；其中，
【权利要求】
1. 一种三维流固单向耦合的实现方法，其特征在于，该方法包括： 将固体边界的三角形网格离散成粒子，并进行边界粒子均匀化采样； 计算该固体中流体粒子的状态，包括：流体与固体边界处流体粒子的密度、流体粒子的 压力及粘力，以及流体粒子表面张力； 并利用一数值积分模式对流体粒子的速度和位置信息进行更新； 当所述流体粒子穿透所述固体边界时，对所述流体粒子的速度和位置进行校正，从而 实现三维流固单向耦合。
2. 根据权利要求1所述的方法，其特征在于，所述将固体边界的三角形网格离散成粒 子，并进行边界粒子均匀化采样包括： 读取固体模型的三角形网格，并建立固体模型的周围符号距离场，具体的：将固体模型 中每个三角形网格建立一个包围盒，并为每个包围盒以统一的方式划分网格；计算每个网 格点到固体模型的符号距离，符号距离场以一个哈希表来存储，形成一个存在于固体模型 周围的窄带宽的符号距离场； 抽取符号距离场的等值面，并在等值面的每个三角形网格上初始化边界粒子，每个三 角形网格上布置gA/πr2个粒子；其中，g表示控制粒子密度的参数，A表示每个三角形网格 的面积，r表示粒子的半径； 根据粒子的分布来决定每个粒子在固体表面法向f和切向^的速度；其中， ..t1ζ___γ-J 巧=-0(1)及，巧=ΣΛ謂创|1 -确;X表示粒子的位置，0(；π和及分别表示 位置X处的符号距离和法向，Q表示一个类高斯核函数，S是整个粒子集合，Π是粒子X周 围第1粒子的位置；通过所述固体表面法向保证了粒子只能在边界表面移动，通过所述 切向€保证了粒子由密集向稀疏的区域移动； 对粒子的速度进行积分来更新粒子位置； 通过给整体粒子集合S的总位移设置一个粒子数量相关的阈值，来判断是否收敛；当 总位移低于阈值时，表明此时迭代过程已经趋于稳定，均匀化的粒子采样已经形成。
3. 根据权利要求1所述的方法，其特征在于，计算该固体中流体粒子状态的公式包括： 流体与固体边界处流体粒子的密度的计算公式为： Pi= XjmjWij+X^bi(P0) ffik； ψμ(ρ〇) =P0Vbi ; 其中，Pi表示流体与固体边界处流体粒子i的密度；IIij表示流体粒子i周围第j个流 体粒子的质量，Wij是一个支持域有限的核函数，PC1表示流体的静止密度，Vbi表示固体边界 粒子bi的体积，k表示固体粒子bi周围的固体粒子； 流体粒子的压力及粘力的计算公式为：
其中，表示固体边界粒子bi对流体粒子i所施加的压力,Iiii表示流体粒子i的质 量，Pi表示流体粒子i的压强，Pi表示流体粒子i的密度，▽Wibi表示核函数Wibi的梯度向 量函数表示固体边界粒子bi对流体粒子i的粘力，μ表示流体的粘度系数，Vfi表 示流体粒子i的速度，▽ 2Wibi表示核函数Wibi的拉普拉斯Laplacian函数； 流体粒子表面张力的计算公式为： Fi^jcur = -ymifci-rij)； 其中，Y表示控制力大小的参数，IHi表示流体粒子i的质量，Iii和rij分别表示流体粒 子i与流体粒子j的法向。
4. 根据权利要求1所述的方法，其特征在于，所述当所述流体粒子穿透所述固体边界 时，对所述流体粒子的速度和位置进行校正包括： 对固体边界区域进行标记，当流体粒子进入标记过的区域时判断流体粒子位置处的符 号距离从而得到流体粒子是否穿透边界； 若穿透，则按照下述公式分别对粒子的速度和位置进行校正： Vi = Vi-(1+Cr) (Vi · n)n； Xi = Xi+(1+Dr) I s |n； 其中，Vi表示流体粒子i的速度，(^表示控制动量损失的参数，值介于O?1之间；Xi 表示流体粒子i的位置，W表示调节回弹距离的一个随机数，s表示Xi处的符号距离，η表 示边界的法向。
5. -种三维流固单向耦合的实现系统，其特征在于，该系统包括： 固体边界粒子采样模块，用于将固体边界的三角形网格离散成粒子，并进行边界粒子 均匀化采样； 流体粒子状态计算模块，用于计算该固体中流体粒子的状态，包括：流体与固体边界处 流体粒子的密度、流体粒子的压力及粘力，以及流体粒子表面张力； 流体粒子信息更新模块，用于并利用一数值积分模式对流体粒子的速度和位置信息进 行更新； 穿透粒子的速度和位置校正模块，用于当所述流体粒子穿透所述固体边界时，对所述 流体粒子的速度和位置进行校正，从而实现三维流固单向耦合。
6. 根据权利要求5所述的系统，其特征在于，所述固体边界粒子采样模块包括： 符号距离场建立模块，用于读取固体模型的三角形网格，并建立固体模型的周围符号 距离场，具体的：将固体模型中每个三角形网格建立一个包围盒，并为每个包围盒以统一的 方式划分网格；计算每个网格点到固体模型的符号距离，符号距离场以一个哈希表来存储， 形成一个存在于固体模型周围的窄带宽的符号距离场； 边界粒子初始化模块，用于抽取符号距离场的等值面，并在等值面的每个三角形网格 上初始化边界粒子，每个三角形网格上布置gA/πr2个粒子；其中，g表示控制粒子密度的 参数，A表示每个三角形网格的面积，r表示粒子的半径； 法向与切向速度决定模块，用于根据粒子的分布来决定每个粒子在固体表面法向f ................t. X-XI 和切向f的速度；其中，|Γ= _0认_K=Σ.mS..p0(||l-詞)IJ-圳I;X表示粒子的位 置，0(Z)和#分别表示位置X处的符号距离和法向，Q表示一个类高斯核函数，S是整个粒 子集合，Π是粒子X周围第1粒子的位置；通过所述固体表面法向f保证了粒子只能在边 界表面移动，通过所述切向保证了粒子由密集向稀疏的区域移动； 位置更新模块，用于对粒子的速度进行积分来更新粒子位置； 均匀化采样模块，用于通过给整体粒子集合S的总位移设置一个粒子数量相关的阈 值，来判断是否收敛；当总位移低于阈值时，表明此时迭代过程已经趋于稳定，均匀化的粒 子采样已经形成。
7. 根据权利要求5所述的系统，其特征在于，所述流体粒子状态计算模块包括： 密度计算模块，流体与固体边界处流体粒子的密度的计算公式为： Pi= XjmjWij+X^bi(P0) ffik； ψμ(ρ〇) =P0Vbi ; 其中，Pi表示流体与固体边界处流体粒子i的密度；IIij表示流体粒子i周围第j个粒 子的质量，Wij是一个支持域有限的核函数，PC1表示流体的静止密度，Vbi表示固体边界粒子 bi的体积，k表示固体粒子bi周围的固体粒子； 压力及粘力计算模块，流体粒子的压力及粘力的计算公式为：
其中，fIw表示固体边界粒子bi对流体粒子i所施加的压力,IIii表示流体粒子i的质 量，Pi表示流体粒子i的压强，Pi表示流体粒子i的密度，▽Wibi表示核函数Wibi的梯度向 量函数；表示固体边界粒子bi对流体粒子i的粘力，μ表示流体的粘度系数，Vfi表 示流体粒子i的速度，▽ 2Wibi表示核函数Wibi的拉普拉斯Laplacian函数； 表面张力计算模块，流体粒子表面张力的计算公式为： Fi^jcur = -ymifci-rij)； 其中，Y表示控制力大小的参数，IHi表示流体粒子i的质量，Iii和rij分别表示流体粒 子i与流体粒子j的法向。
8. 根据权利要求5所述的系统，其特征在于，所述穿透粒子的速度和位置校正模块包 括： 穿透判断模块，用于对固体边界区域进行标记，当流体粒子进入标记过的区域时判断 流体粒子位置处的符号距离从而得到流体粒子是否穿透边界； 校正模块，用于穿透时，则按照下述公式分别对粒子的速度和位置进行校正： Vi = Vi-(1+Cr) (Vi · n)n ; Xi = Xi+(1+Dr) I s |n； 其中，Vi表示流体粒子i的速度，(^表示控制动量损失的参数，值介于O?I之间；Xi 表示流体粒子i的位置，W表示调节回弹距离的一个随机数，s表示Xi处的符号距离，η表 示边界的法向。
【文档编号】G06T13/20GK104318598SQ201410553159
【公开日】2015年1月28日 申请日期:2014年10月17日 优先权日:2014年10月17日
【发明者】董兰芳, 章恒 申请人:中国科学技术大学