一种sph多相流触觉交互研究的方法
【专利摘要】本发明公开了一种SPH多相流触觉交互研究的方法，该方法包括：步骤1、SPH混合流模拟；步骤2、实现SPH流体的固/液耦合；步骤3、实现触觉交互操控与触觉力渲染。与现有技术相比，与现有技术相比，本发明可以为用户提供在虚拟场景中流体交互与触觉力真实反馈，可以用于虚拟现实教育领域与游戏互动等，在视觉和触觉上为用提供更真实的虚拟现实交互体验；通过本发明获得了自由度控制方式和3自由度力反馈的设备性能与虚拟场景中的流体交互及感受流体反馈而来的力。
【专利说明】
-种SPH多相流触觉交互研究的方法
技术领域
[0001] 本发明属于计算机虚拟现实应用中关于触觉交互的流体仿真技术，特别是设及一 种SPH多相流触觉交互研究的方法。
【背景技术】
[0002] 在计算机图形学领域中，对基于物理的流体控制与交互的研究在游戏动画、数字 娱乐、虚拟仿真等领域具有较强的应用价值。在虚拟现实教育和游戏中，常常需要在流体交 互时进行揽拌，如游戏中的炼金，炼药动作。在基于物理的流体动画制作中，基于触觉设备 的流体交互模拟也是一个研究热点。
[0003] 触觉交互设备的发展与成熟，让流体控制与交互能够在虚拟场景交互的过程中得 到真实的力反馈效果。触觉交互设备不仅可W控制选择的对象在Ξ维场景中移动，还能够 让用户通过设备来感受真实的触觉力，因此被广泛应用于虚拟现实领域如手术模拟EU、3D 雕刻W、虚拟场景搭建W等方面。在虚拟现实教育和游戏中，常常需要对于虚拟的流体进行 交互。相对于单相流，对多相流的交互在实现上的困难在于触觉交互要能够让用户分辨出 各种流体性质的不同，并且在加入了混溶的物理模拟方案后仍能有很好的实时性。
[0004] Cirio等人[4'5]使用使用HAPTI0N Vbtuose 6D35-45实现了基于SPH方法的固体、 液体及固/液混合态的物体的交互力反馈。Menelas等人W采用该设备实现了非稳态流中的 力的模拟。Yang等人W通过Novint化Icon控制球形刚体在流体中运动，可W实现水流中的 反馈力的计算;该设备为Ξ自由度输入，在实际应用中的局限性较大。Pier等人W对SW1方 法在流体交互中的计算进行了并行化加速，提升了运行效率。相对于前人实时计算反馈力 的策略，Dobashi等人W使用自制的设备Spidar-G完成了虚拟划艇场景搭建，能够通过划奖 的反馈力真实体验划艇运动;该平台可W通过离线计算的方式来提高反馈力的频率获得更 好的体验，但相应的该方法在可拓展性上存在一定的问题。Geomagic系列是现在使用最广 泛的触觉交互设备，有各种精度级别的设备，其中Geomagic touch(原地antom 0皿i)作为 入口级6自由度触觉交互设备被广为使用。
[0005] 目前基于触觉设备的流体人机交互方法主要分为网格和SPH方法两种。Yang等 人W使用网格的方法，并针对实时性做了优化和评价;该研究指出了实时性对触觉交互的 重要性。文献W'W同样使用网格的方法，基于Geomagic系列的触觉交互设备3模拟了在容器 中揽动流体的反馈效果，但该方法对于流体混合的处理并非真正Ξ维空间下的模拟。而文 献中使用基本SPH方法模拟各种状态的流体和固体，通过SPH haptic force来反馈触觉 力给输入设备，之后通过CUDA并行计算来加速模拟效率，最终达到较好的实时性。该研究成 功模拟了使用锅等厨具互相倾倒流体的真实触觉效果。Ma等人中将设备用于肝脏模型 的学习。
[0006] 相对于网格方法，基于SPH(光滑粒子动力学)方法的触觉交互上具有独特的优势。 SPH方法不仅可W模拟流体，也可W模拟固体、粒子类型扩展非常方便，可W无缝完成流固 禪合平台，能够很好的解决模拟平滑稳定的触觉力、多触点处理、计算的并行性和可扩展性 等触觉模拟研究中很重要的问题。现有的流体交互技术在多相流领域，尤其是对基于SKI流 体的多相流触觉交互缺少研究。
[0007] [1]曾妍文、许天春、岳龙旺等;《带有力感觉的显微外科手术血管的仿真研究》计 算机仿真，2006,23(7) :214-218;
[0008] [ 2]陈飞;《基于触觉反馈的虚拟变形与切割研究》，哈尔滨工业大学，2009;
[0009] [3]Dobashi Y,Yamamoto Τ,Sato M,et al.A precomputed approach for re al-time haptic interaction with fluids.Computer Graphics and Applications, 2007,27(3):90-92;
[0010] [4]Cirio G，Marchal M，0taduy M A,et al.Six-Dof haptic interaction with fluids，solids， and their transitions.Proceedings of World Haptics Conference (WHO,2013:157-162；
[0011] [5]Cirio G，Marchal M，Hillaire S，et al. Six degrees-of-freedom haptic interaction with fluids.Visualization and Computer Graphics,2011,17(11):1714- 1727;
[0012] [6]Menelas B，Ammi M,Pastur L，et al.Haptical exploration of an unsteady flow.Proceedings of EuroHaptics conference，2009:232-237;
[0013] [7]Yang M，Lu J,Safonova A，et al.GPU methods for real-time haptic interaction with 3D fluids. Proceedings of Haptic Audio visual Environments and Games,2009:24-29;
[0014] [8]Pier J M，Figueroa I，Huegel J.CUDA-enabled particle-based 3D fluid haptic simulation . Proceedings of Electronics，民obotics and Automotive Mechanics Conference(CERMA),2011:391-396；
[001 日] [9]Mora J, Lee W S.民eal-time fluid interaction with a haptic device.Proceedings of Haptic,Audio and Visual Environments and Games,2007: 160-165;
[0016] [lOjMora J,Lee W S.Real-time 3D fluid interaction with a haptic user inter face.Proceedings of 3D User Interfaces，2008:75-81;
[0017] [ll]Ma N,Liu Y,Qiao A，et al.Design of three-dimensional interactive visualization system based on force feedback device.Proceedings of Bioinformatics and Biomedical Engineering,2008:1780-1783〇
【发明内容】
[0018] 基于上述现有技术和存在的问题，本发明提出了一种SPH多相流触觉交互研究的 方法，针对现在基于SPH方法的触觉交互中仅有流固禪合^及固体与不可溶流体交互而没 有基于多相流的触觉模拟的问题进行研究，实现基于SKI的多相流触觉反馈的真实模拟效 果D
[0019] 本发明提出了一种SPH多相流触觉交互研究的方法，该方法包括W下步骤：
[0020] 步骤1、实现SPH纔合流模拟，具体包括W下处理：
[0021] 将各相流体纔合后得到的纔合流体模拟求解的问题转化为每个时间步长对于SPH 粒子的各相流体积分数Qk和混合速度Um的求解：
[0024] 其中，Pm为混合流体的压强，Pm为混合流体的密度，g为重力加速度，t为模拟的时间 步长，▽· Tm为混合流体的粘性力，▽· Tom为混合流体中的各相流体之间由于漂移速度导 致的动量变化，Umk为混合流体中的各相流体的漂移速度，k为流体的相数，Um · ▽表示混合 速度与梯度算子相乘，▽表示梯度算子，Vpm表示混合流体压强的梯度;其中Umk需要满足W 下方程组：
[0025]
[0026] 公式(3)表示混合流体中每一点对应的各相流体漂移速度之和为零;公式(4)则是 该漂移速度的求解公式;其中，Odif为各相流体的扩散系数，k '为其他相流体，Vak表示体积 分数的梯度，表示其他相流体的静态密度在本发明中扩散系数取0.1有较好的效果。
[0027] 压力项Vpm的计算方式为：
[002引
[0029] 其中，为混合流体中刚体粒子i的压强，刚体粒子j的压强。
[0030] 对各相流之间由于漂移速度导致的混合流体的动量变化▽· TDm的计算公式为：
[0031]
[0032] 其中，WU扣h)为高斯核函数，h为高斯核函数的作用半径，w = (xi-xj)是SPH粒子 i与刚体粒子j的距离，Umki为混合流体中SPH粒子i的流体漂移速度，Umw为混合流体中刚体 粒子j的流体漂移速度，：?功流体SPH粒子i的各相流体积分数ak，aw为刚体粒子j的各相流 体积分数Qk，mj表示其他粒子j的质量；
[0033] 对于粘性力的求解，由于各个流体粒子的粘性力随着各相流所占的百分比的变化 而变化，流体粒子的粘性力由下式求解：
[0034]
[00对其中，μι为SPH粒子i的粘性系数，yj为刚体粒子j的粘性系数，Umi为SPH粒子i混合 后的速度，Umj为刚体粒子j混合后的速度，W表示核函数，1?表示其他粒子j的质量。
[0036] SPH粒子与动态刚体交互所受到的触觉力fhaPtie在本发明中可分解为fPTes与fvis两 项，fiPtes为刚体粒子所受到的压力，为与刚体粒子交互的粘性力，在多相流模拟中一旦 发生了流体与动态刚体的交互，则对其进行求解：
[0039] 在公式(7)、（8)中，i为SPH粒子，j为刚体粒子，p功刚体粒子j的压强，Pm功混合流 体中刚体粒子i的压强，叫为刚体粒子j的粘性系数，μ?ι为混合流体中SPH粒子i的粘性系数， Uj为刚体粒子的速度，Umi为SPH粒子i混合后的速度，Pj为刚体粒子j的密度，Pmi为刚体粒子 的压强，1?表示其他粒子j的质量，▽表示梯度算子。在本发明的模拟步骤中，通过对上述公 式求解，能够较好的模拟触觉交互时多相流的混合扩散效果，并且有较好的实时性，可W实 现实时交互控制的任务。
[0040] 步骤2、实现SPH流体的固/液禪合，具体包括W下处理：
[0041] 通过计算所有刚体粒子所受的合力FhaPtK作为触觉交互的反馈力，其中：
[0042] 刚体粒子受到的触觉力公式表示为压力与粘性力：
[0049] 步骤3、实现触觉交互操控与触觉力擅染，具体包括W下处理：
[0050] 将刚体的位移与旋转表现为每一个刚体粒子的空间位置相对质屯、的修正过程。通 过触觉设备的位移更新质屯、的位置，通过触觉设备的旋转更新刚体的旋转角速度;使用带 汉明窗的FIR(Finite Impulse Response)滤波器将离散的低频率的触觉力数据平滑化为 1000赫兹的高频率信息。
[0051] 与现有技术相比，本发明可W为用户提供在虚拟场景中流体交互与触觉力真实反 馈，可W用于虚拟现实教育领域与游戏互动等，在视觉和触觉上为用提供更真实的虚拟现 实交互体验;通过本发明获得了自由度控制方式和3自由度力反馈的设备性能与虚拟场景 中的流体交互及感受流体反馈而来的力。
【附图说明】
[0052] 图1为实现本发明的SK1多相流触觉交互研究的方法的整体框架图。
[0053] 图2为使用触觉设备控制虚拟场景中物体与不可混溶流体交互的效果；
[0054] 图3为使用触觉设备控制虚拟场景中物体与可混溶流体交互的效果。
【具体实施方式】
[0055] W下结合附图及【具体实施方式】，进一步详述本发明的技术方案。
[0056] 本发明使用的触觉交互设备是Geomagic touch,该设备具有六自由度的输入与Ξ 自由度的反馈力输出。在本发明的流体模拟中，W应用CUDA框架的SPH多相流作为基础，并 且加入与动态刚体的固/液禪合的方法和触觉擅染的方法和交互的方法，实现该方法的整 体框架如图1所示。
[0057] 本发明的具体实施流程，包括W下步骤：
[0化引1、实现SP取昆合流模拟
[0059] 在流体模拟中，对于可混溶的液体，需要考虑速度场不均匀分布所带来的体积分 数变化。本发明改进SP取昆合流模拟的方法用于触觉交互，W实现混合流体揽拌的效果。
[0060] 本发明将各相流体混合后得到的混合流体模拟求解的问题转化为每个时间步长 对于SPH粒子的各相流体积分数Qk和混合速度Um的求解：
[0063] 其中，Pm为混合流体的压强，Pm为混合流体的密度，g为重力加速度，t为模拟的时间 步长，▽· Tm为混合流体的粘性力，▽· Tom为混合流体中的各相流体之间由于漂移速度导 致的动量变化，Umk为混合流体中的各相流体的漂移速度，k为流体的相数，，Um · ▽表示混合 速度与梯度算子相乘，▽表示梯度算子，Vpm表示混合流体压强的梯度。公式（1)与（2)分别 转化自质量守恒与动量守恒。其中，公式(2)加入了触觉力项fhaPtK，表示动态刚体对混合流 体的作用力，该项在产生固/液禪合时才做计算。
[0064] 求解各相流体体积分数ak的关键在于对混合流体中的各相流体的漂移速度Umk的 求解，考虑到运算的效率W及前两项相对于扩散效果来说较为不明显，因此在计算漂移速 度时只考虑了扩散效应的计算，本模拟方法中U血需要满足W下方程组：
[00 化]
[0066] 公式(3)表示混合流体中每一点对应的各相流体漂移速度之和为零;公式(4)则是 该漂移速度的求解公式;其中，Odif为各相流体的扩散系数，k '为其他相流体，Vak表示体积 分数的梯度，表示其他相流体的静态密度，在本发明中扩散系数取0.1有较好的效果。
[0067] 压力项Vpm的计算方式与参照基本的SPH方法中单相流压力项的计算方式
[006引
[0069] 其中，为混合流体中刚体粒子i的压强，刚体粒子j的压强。
[0070] 对各相流之间由于漂移速度导致的混合流体的动量变化▽· TDm的计算公式为：
[0071]
[0072] 其中，W(x扣h)为高斯核函数，h为高斯核函数的作用半径，w = (xi-xj)是SPH粒子 i与刚体粒子j的距离，Umki为混合流体中SPH粒子i的流体漂移速度，Umw为混合流体中刚体 粒子j的流体漂移速度，：?功流体SPH粒子i的各相流体积分数ak，aw为刚体粒子j的各相流 体积分数ak，ny表示其他粒子j的质量。对于粘性力的求解，由于各个流体粒子的粘性力随着 各相流所占的百分比的变化而变化，流体粒子的粘性力由下式求解：
[0073]
[0074] 其中，μι为SPH粒子i的粘性系数，μ功刚体粒子j的粘性系数，Umi为SPH粒子i混合 后的速度，Umj为刚体粒子j混合后的速度，W表示核函数，mj表示其他粒子j的质量。
[0075] SPH粒子与动态刚体交互所受到的触觉力fhaPtie在本发明中可分解为fPTes与fvis两 项，fiPtes为刚体粒子所受到的压力，为与刚体粒子交互的粘性力，在多相流模拟中一旦 发生了流体与动态刚体的交互，则对其进行求解：
[0078]在公式(7)、（8)中，i为SPH粒子，j为刚体粒子，PJ为刚体粒子j的压强，Pm功混合流 体中刚体粒子i的压强，叫为刚体粒子j的粘性系数，μ?ι为混合流体中SPH粒子i的粘性系数， UJ为刚体粒子的速度，Umi为SPH粒子i混合后的速度，P功刚体粒子j的密度，Pmi为刚体粒子 的压强，1?表示其他粒子j的质量，▽为梯度算子。在本发明的模拟步骤中，通过对上述公式 求解，能够较好的模拟触觉交互时多相流的混合扩散效果，并且有较好的实时性，可W实现 实时交互控制的任务。在该步骤中，本发明基于混合流模型，将应用于SK1方法的多相流模 拟公式简化，保留最重要的扩散效应计算，在提高运行效率的同时保持逼真的多相流模拟 现象，并应用触觉交互的方法。
[00巧]2、实现SPH流体的固/液禪合
[0080] 在流体模拟中，需要与流体发生交互作用的实体大致包括两类，静态实体和动态 实体。相对于静态实体，动态刚体与流体粒子的禪合的问题相对来说更为复杂。本发明利用 触觉设备来操控动态刚体与流体交互，并将反馈力实时传回给操控者。因此，反馈力的精确 计算显得尤为重要。
[0081] 为保证交互与计算的精确性，对于该步骤，本发明使用统一粒子模型的方法，即： 将刚体看作结构稳定的刚体粒子的结合，使用体素化的方法将刚体转化为集合粒子，并且 使用高斯核函数来计算各种粒子之间的相互作用。
[0082] 将通过计算所有刚体粒子所受的合力FhaptK作为触觉交互的反馈力。对刚体粒子 压强的pH十算使用W下公式W防止渗流现象产生：
[0083]
[0084] 而流体粒子的压强pf也同样做运样的处理：
[0085]
[0086] 其中，p为动态密度，ρτ和pf为上述刚体粒子与流体的静态密度，kt与kf为压强系 数;运里kf取1.5，kt取3.0。运两个公式使得压力的计算始终保持正值形成排斥力防止刚体 与流体粒子相互渗流。产生的流体粒子进入刚体粒子，不符合实际。
[0087] 刚体粒子受到的触觉力公式表示为压力与粘性力：
[008引
（U)
[0089]其中，压力fiPns与粘性力fiVis的计算公式如下：
[0092]最后计算合力，可W反馈给触觉设备作为反馈力输出：
[oow]
片4)。
[0094]在该步骤中，本发明使用体素化的方法将网格模型转化为刚体粒子的集合，使用 统一粒子模型的方法对运些刚体粒子赋予初始属性和状态W后让其与流体粒子同样进行 状态更新，并在与流体交互的过程中保持刚体结构的稳定性且防止渗流效果;并计算所有 刚体粒子所受的合力作为触觉交互的反馈力。
[00M] 3、实现触觉交互操控与触觉力擅染
[0096] 动态刚体的运动可W看作是刚体围绕质屯、进行平移和旋转，并且不会产生形变。 本发明中操作触觉设备控制虚拟场景中粒子表示的动态刚体运动可W转化为另一种思路： 将刚体的位移与旋转表现为每一个刚体粒子的空间位置相对质屯、的修正过程。
[0097] 假设所有粒子的质量相同，对于一个刚体粒子，速度公式可W表示为：
[009引 Ui = Ur+C0rXqi (15)
[0099] 其中，Ur是刚体质屯、的速度，ωτ是刚体的角速度，Ui是单个刚体粒子j的速度，qi是 粒子ij与质屯、的相对位置。在程序框架中通过每一个时间步长更新交互设备的位置与旋转 角度的输入参数来实时的控制虚拟场景中动态刚体的运动。
[0100] 在触觉力擅染部分，本发明将触觉力的更新从模拟步骤分离开并且通过滤波来填 充模拟步骤之间过渡触觉力实现平滑模拟的效果。本发明使用带汉明窗的FIR(Finite Impulse Response)滤波器将离散的低频率的触觉力数据平滑化为1000赫兹的高频率信 息，运样可W防止使用者在操作过程中感受到震动。虽然FIR滤波器的特性会带来一定的反 馈力延迟，但是由于延迟时间很短，在实际使用时很难发现。经过实验发现，使用滤波的方 法将触觉力平滑化效果好于未平滑的初始效果。
[0101] 本发明实验机器配置如下：Intel Core i7 4700mq，GTX765m，16G RAM,Windows 8.1，CUDA 7.0，VS2013。并支持Windows 7，Windows8.1，Windows 10环境，IG内存W上， nVidia GT430W上显卡配置的微机。
【主权项】
1. 一种SPH多相流触觉交互研究的方法，其特征在于，该方法包括以下步骤： 步骤1、实现SPH混合流模拟，具体包括以下处理： 将各相流体混合后得到的混合流体模拟求解的问题转化为每个时间步长对于SPH粒子 的各相流体积分数〇k和混合速度um的求解：其中，Pm为混合流体的压强，Pm为混合流体的密度，g为重力加速度，t为模拟的时间步 长，V ?；为混合流体的粘性力，力^为混合流体中的各相流体之间由于漂移速度导致的 动量变化，Umk为混合流体中的各相流体的漂移速度，k为流体的相数，& · V表示混合速度与 梯度算子相乘，▽表示梯度算子，V/?,,表示混合流体压强的梯度;其中umk需要满足以下方程 组：公式(3)表示混合流体中每一点对应的各相流体漂移速度之和为零;公式(4)则是该漂 移速度的求解公式;其中，〇dif为各相流体的扩散系数，k '为其他相流体，▽%表示体积分数 的梯度，表示其他相流体的静态密度，在本发明中扩散系数取0.1有较好的效果。 压力项V/7.,,的计算方式为：其中，为混合流体中刚体粒子i的压强，刚体粒子j的压强； 对各相流之间由于漂移速度导致的混合流体的动量变化V · 的计算公式为：其中，W(Xji，h)为高斯核函数，h为高斯核函数的作用半径，Xij= (Xi-Xj)是SPH粒子i与 刚体粒子j的距离，umkl为混合流体中SPH粒子i的流体漂移速度，umkj为混合流体中刚体粒子 j的流体漂移速度，：aki为流体SPH粒子i的各相流体积分数ak，akj为刚体粒子j的各相流体积 分数ak，mj表示其他粒子j的质量;流体粒子的粘性力由下式求解：其中，Pi为SPH粒子i的粘性系数，μ」为刚体粒子j的粘性系数，umi为SPH粒子i混合后的速 度，叫为刚体粒子j混合后的速度，W表示核函数，m谦示其他粒子j的质量。 SPH粒子与动态刚体交互所受到的触觉力fhaptle在本发明中可分解为fpi:es与f vls两项， 为刚体粒子所受到的压力，为与刚体粒子交互的粘性力，在多相流模拟中一旦发 生了流体与动态刚体的交互，则对其进行求解：在公式(7)、（8)中，i为SPH粒子，j为刚体粒子，p伪刚体粒子j的压强，Pmi为混合流体中 刚体粒子i的压强，μ伪刚体粒子j的粘性系数，为混合流体中SPH粒子i的粘性系数，1!伪 刚体粒子的速度，um^SPH粒子i混合后的速度，P伪刚体粒子j的密度，p mi为刚体粒子的压 强，πυ为刚体粒子的质量，V表示梯度算子。在本发明的模拟步骤中，通过对上述公式求解， 能够较好的模拟触觉交互时多相流的混合扩散效果，并且有较好的实时性，可以实现实时 交互控制的任务。 步骤2、实现SPH流体的固/液耦合，具体包括以下处理： 通过计算所有刚体粒子所受的合力Fhaptlc;作为触觉交互的反馈力，其中： 刚体粒子受到的触觉力公式表示为压力与粘性力：其中，压力与粘性力的计算公式如下：最后计算合力，可以反馈给触觉设备作为反馈力输出：步骤3、实现触觉交互操控与触觉力渲染，具体包括以下处理： 将刚体的位移与旋转表现为每一个刚体粒子的空间位置相对质心的修正过程。通过触 觉设备的位移更新质心的位置，通过触觉设备的旋转更新刚体的旋转角速度;使用带汉明 窗的FIR(Finite Impulse Response)滤波器将离散的低频率的触觉力数据平滑化为1000 赫兹的高频率信息。
【文档编号】G06F19/00GK105825059SQ201610152586
【公开日】2016年8月3日
【申请日】2016年3月17日
【发明人】刘世光, 张枭勇
【申请人】天津大学