一种SPH多相流触觉交互研究的方法

摘要

本发明公开了一种SPH多相流触觉交互研究的方法，该方法包括：步骤1、SPH混合流模拟；步骤2、实现SPH流体的固/液耦合；步骤3、实现触觉交互操控与触觉力渲染。与现有技术相比，与现有技术相比，本发明可以为用户提供在虚拟场景中流体交互与触觉力真实反馈，可以用于虚拟现实教育领域与游戏互动等，在视觉和触觉上为用提供更真实的虚拟现实交互体验；通过本发明获得了自由度控制方式和3自由度力反馈的设备性能与虚拟场景中的流体交互及感受流体反馈而来的力。

权利要求

1.一种SPH多相流触觉交互研究的方法，其特征在于，该方法包括以下步骤：

步骤(1)、实现SPH混合流模拟，具体包括以下处理：

将各相流体混合后得到的混合流体模拟求解的问题转化为每个时间步长对于SPH粒子的各相流体积分数αk和混合速度um的求解：其中，pm为混合流体的压强，ρm为混合流体的密度，g为重力加速度，t为模拟的时间步长，▽·Tm为混合流体的粘性力，▽·TDm为混合流体中的各相流体之间由于漂移速度导致的动量变化，umk为混合流体中的各相流体的漂移速度，k为流体的相数,um·▽表示混合速度与梯度算子相乘，▽表示梯度算子，▽pm表示混合流体压强的梯度，fhaptic表示SPH粒子与动态刚体交互所受到的触觉力；其中umk需要满足以下方程组：公式(3)表示混合流体中每一点对应的各相流体漂移速度之和为零；公式(4)则是该漂移速度的求解公式；其中，σdif为各相流体的扩散系数，k'为其他相流体，▽αk表示体积分数的梯度，ρk′表示其他相流体的静态密度；

压力项▽pm的计算方式为：

其中，Pmi为混合流体中刚体粒子i的压强，Pmj为混合流体中刚体粒子j的压强；Pj为刚体粒子j的压强；

对各相流之间由于漂移速度导致的混合流体的动量变化▽·TDm的计算公式为：

其中，W(xji,h)为高斯核函数，h为高斯核函数的作用半径，xij＝(xi-xj)是SPH粒子i与刚体粒子j的距离，umki为混合流体中SPH粒子i的流体漂移速度，umkj为混合流体中刚体粒子j的流体漂移速度，αki为流体SPH粒子i的各相流体积分数αk，αkj为刚体粒子j的各相流体积分数αk，mj表示其他粒子j的质量；流体粒子的粘性力由下式求解：其中，μi为SPH粒子i的粘性系数，μj为刚体粒子j的粘性系数，umi为SPH粒子i混合后的速度，umj为刚体粒子j混合后的速度，W表示核函数，mj表示其他粒子j的质量，ρj表示刚性粒子j的密度；

SPH粒子与动态刚体交互所受到的触觉力fhaptic分解为fpres与fvis两项，fipres为刚体粒子所受到的压力，fivis为与刚体粒子交互的粘性力，在多相流模拟中一旦发生了流体与动态刚体的交互，则对其进行求解：在公式(7)、(8)中，i为SPH粒子，j为刚体粒子，pj为刚体粒子j的压强，Pmi为混合流体中刚体粒子i的压强，μj为刚体粒子j的粘性系数，μmi为混合流体中SPH粒子i的粘性系数，uj为刚体粒子的速度，umi为SPH粒子i混合后的速度，ρj为刚体粒子j的密度，Pmi为刚体粒子的压强，mj为刚体粒子的质量，▽表示梯度算子；在模拟步骤中，通过对上述公式求解，模拟触觉交互时多相流的混合扩散效果，并具有实时性，实现交互控制的任务；

步骤(2)、实现SPH流体的固/液耦合，具体包括以下处理：

通过计算所有刚体粒子所受的合力Fhaptic作为触觉交互的反馈力，其中：

刚体粒子受到的触觉力公式表示为压力与粘性力：

fihaptic＝fipres+fivis    (11)

其中，压力fipres与粘性力fivis的计算公式如下：

最后计算合力，反馈给触觉设备作为反馈力输出：

步骤(3)、实现触觉交互操控与触觉力渲染，具体包括以下处理：

将刚体的位移与旋转表现为每一个刚体粒子的空间位置相对质心的修正过程；通过触觉设备的位移更新质心的位置，通过触觉设备的旋转更新刚体的旋转角速度；使用带汉明窗的FIR(Finite Impulse Response)滤波器将离散的低频率的触觉力数据平滑化为1000赫兹的高频率信息。

说明书

一种SPH多相流触觉交互研究的方法

技术领域

[0001] 本发明属于计算机虚拟现实应用中关于触觉交互的流体仿真技术，特别是涉及一种SPH多相流触觉交互研究的方法。

背景技术

[0002] 在计算机图形学领域中，对基于物理的流体控制与交互的研究在游戏动画、数字娱乐、虚拟仿真等领域具有较强的应用价值。在虚拟现实教育和游戏中，常常需要在流体交互时进行搅拌，如游戏中的炼金，炼药动作。在基于物理的流体动画制作中，基于触觉设备的流体交互模拟也是一个研究热点。

[0003] 触觉交互设备的发展与成熟，让流体控制与交互能够在虚拟场景交互的过程中得到真实的力反馈效果。触觉交互设备不仅可以控制选择的对象在三维场景中移动，还能够让用户通过设备来感受真实的触觉力，因此被广泛应用于虚拟现实领域如手术模拟[1]、3D雕刻[2]、虚拟场景搭建[3]等方面。在虚拟现实教育和游戏中，常常需要对于虚拟的流体进行交互。相对于单相流，对多相流的交互在实现上的困难在于触觉交互要能够让用户分辨出各种流体性质的不同，并且在加入了混溶的物理模拟方案后仍能有很好的实时性。

[0004] Cirio等人[4，5]使用使用HAPTION Virtuose 6D35-45实现了基于SPH方法的固体、液体及固/液混合态的物体的交互力反馈。Menelas等人[6]采用该设备实现了非稳态流中的力的模拟。Yang等人[7]通过Novint Falcon控制球形刚体在流体中运动，可以实现水流中的反馈力的计算；该设备为三自由度输入，在实际应用中的局限性较大。Pier等人[8]对SPH方法在流体交互中的计算进行了并行化加速，提升了运行效率。相对于前人实时计算反馈力的策略，Dobashi等人[3]使用自制的设备Spidar-G完成了虚拟划艇场景搭建，能够通过划桨的反馈力真实体验划艇运动；该平台可以通过离线计算的方式来提高反馈力的频率获得更好的体验，但相应的该方法在可拓展性上存在一定的问题。Geomagic系列是现在使用最广泛的触觉交互设备，有各种精度级别的设备，其中Geomagic touch(原phantom omni)作为入门级6自由度触觉交互设备被广为使用。

[0005] 目前基于触觉设备的流体人机交互方法主要分为网格和SPH方法两种。Yang等人[7]使用网格的方法，并针对实时性做了优化和评价；该研究指出了实时性对触觉交互的重要性。文献[9,10]同样使用网格的方法，基于Geomagic系列的触觉交互设备]模拟了在容器中搅动流体的反馈效果，但该方法对于流体混合的处理并非真正三维空间下的模拟。而文献[4,5]中使用基本SPH方法模拟各种状态的流体和固体，通过SPH haptic force来反馈触觉力给输入设备，之后通过CUDA并行计算来加速模拟效率，最终达到较好的实时性。该研究成功模拟了使用锅等厨具互相倾倒流体的真实触觉效果。Ma等人[11]中将设备用于肝脏模型的学习。

[0006] 相对于网格方法，基于SPH(光滑粒子动力学)方法的触觉交互上具有独特的优势。SPH方法不仅可以模拟流体，也可以模拟固体、粒子类型扩展非常方便，可以无缝完成流固耦合平台，能够很好的解决模拟平滑稳定的触觉力、多触点处理、计算的并行性和可扩展性等触觉模拟研究中很重要的问题。现有的流体交互技术在多相流领域，尤其是对基于SPH流体的多相流触觉交互缺少研究。

[0007] [1]曾妍文、许天春、岳龙旺等；《带有力感觉的显微外科手术血管的仿真研究》计算机仿真，2006，23(7)：214-218；

[0008] [2]陈飞《；基于触觉反馈的虚拟变形与切割研究》，哈尔滨工业大学，2009；

[0009] [3]Dobashi Y,Yamamoto T,Sato M,et al.Aprecomputed approach for re al-time haptic interaction with fluids.Computer Graphics and Applications,2007,27(3):90-92；

[0010] [4]Cirio G,Marchal M,Otaduy M A,et al.Six-Dof haptic interaction with fluids,solids,and their transitions.Proceedings of World Haptics Conference(WHC),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,Robotics and Automotive Mechanics Conference(CERMA),2011:391-396；

[0015] [9]Mora J,Lee W  S.Real-time fluid interaction with a haptic device.Proceedings of Haptic,Audio and Visual Environments and Games,2007:160-165；

[0016] [10]Mora J,Lee W S.Real-time 3D fluid interaction with a haptic user inter face.Proceedings of 3D User Interfaces,2008:75-81；

[0017] [11]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方法的触觉交互中仅有流固耦合以及固体与不可溶流体交互而没有基于多相流的触觉模拟的问题进行研究，实现基于SPH的多相流触觉反馈的真实模拟效果。

[0019] 本发明提出了一种SPH多相流触觉交互研究的方法，该方法包括以下步骤：

[0020] 步骤1、实现SPH混合流模拟，具体包括以下处理：

[0021] 将各相流体混合后得到的混合流体模拟求解的问题转化为每个时间步长对于SPH粒子的各相流体积分数αk和混合速度um的求解：

[0022]

[0023]

[0024] 其中，pm为混合流体的压强，ρm为混合流体的密度，g为重力加速度，t为模拟的时间步长，▽·Tm为混合流体的粘性力，▽·TDm为混合流体中的各相流体之间由于漂移速度导致的动量变化，umk为混合流体中的各相流体的漂移速度，k为流体的相数,um·▽表示混合速度与梯度算子相乘，▽表示梯度算子，▽pm表示混合流体压强的梯度，fhaptic表示SPH粒子与动态刚体交互所受到的触觉力；其中umk需要满足以下方程组：

[0025]

[0026] 公式(3)表示混合流体中每一点对应的各相流体漂移速度之和为零；公式(4)则是该漂移速度的求解公式；其中，σdif为各相流体的扩散系数，k'为其他相流体，▽αk表示体积分数的梯度，ρk′表示其他相流体的静态密度；

[0027] 压力项▽pm的计算方式为：

[0028]

[0029] 其中，Pmi为混合流体中刚体粒子i的压强，Pmj为混合流体中刚体粒子j的压强；Pj为刚体粒子j的压强；

[0030] 对各相流之间由于漂移速度导致的混合流体的动量变化▽·TDm的计算公式为：

[0031]

[0032] 其中，W(xji,h)为高斯核函数，h为高斯核函数的作用半径，xij＝(xi-xj)是SPH粒子i与刚体粒子j的距离，umki为混合流体中SPH粒子i的流体漂移速度，umkj为混合流体中刚体粒子j的流体漂移速度， 为流体SPH粒子i的各相流体积分数αk，αkj为刚体粒子j的各相流体积分数αk，mj表示其他粒子j的质量；流体粒子的粘性力由下式求解：

[0033]

[0034] 其中，μi为SPH粒子i的粘性系数，μj为刚体粒子j的粘性系数，umi为SPH粒子i混合后的速度，umj为刚体粒子j混合后的速度，W表示核函数，mj表示其他粒子j的质量，ρj表示刚性粒子j的密度；

[0035] SPH粒子与动态刚体交互所受到的触觉力fhaptic在本发明中可分解为fpres与fvis两项，fipres为刚体粒子所受到的压力，fivis为与刚体粒子交互的粘性力，在多相流模拟中一旦发生了流体与动态刚体的交互，则对其进行求解：

[0036]

[0037]

[0038] 在公式(7)、(8)中，i为SPH粒子，j为刚体粒子，pj为刚体粒子j的压强，Pmi为混合流体中刚体粒子i的压强，μj为刚体粒子j的粘性系数，μmi为混合流体中SPH粒子i的粘性系数，uj为刚体粒子的速度，umi为SPH粒子i混合后的速度，ρj为刚体粒子j的密度，Pmi为刚体粒子的压强，mj为刚体粒子的质量，▽表示梯度算子；在本发明的模拟步骤中，通过对上述公式求解，模拟触觉交互时多相流的混合扩散效果，并具有实时性，实现交互控制的任务；

[0039] 步骤2、实现SPH流体的固/液耦合，具体包括以下处理：

[0040] 通过计算所有刚体粒子所受的合力Fhaptic作为触觉交互的反馈力，其中：

[0041] 刚体粒子受到的触觉力公式表示为压力与粘性力：

[0042] fihaptic＝fipres+fivis       (11)

[0043] 其中，压力fipres与粘性力fivis的计算公式如下：

[0044]

[0045]

[0046] 最后计算合力，反馈给触觉设备作为反馈力输出：

[0047]

[0048] 步骤3、实现触觉交互操控与触觉力渲染，具体包括以下处理：

[0049] 将刚体的位移与旋转表现为每一个刚体粒子的空间位置相对质心的修正过程；通过触觉设备的位移更新质心的位置，通过触觉设备的旋转更新刚体的旋转角速度；使用带汉明窗的FIR(Finite Impulse Response)滤波器将离散的低频率的触觉力数据平滑化为1000赫兹的高频率信息。

[0050] 与现有技术相比，本发明可以为用户提供在虚拟场景中流体交互与触觉力真实反馈，可以用于虚拟现实教育领域与游戏互动等，在视觉和触觉上为用提供更真实的虚拟现实交互体验；通过本发明获得了自由度控制方式和3自由度力反馈的设备性能与虚拟场景中的流体交互及感受流体反馈而来的力。

附图说明

[0051] 图1为实现本发明的SPH多相流触觉交互研究的方法的整体框架图。

[0052] 图2为使用触觉设备控制虚拟场景中物体与不可混溶流体交互的效果；

[0053] 图3为使用触觉设备控制虚拟场景中物体与可混溶流体交互的效果。

具体实施方式

[0054] 以下结合附图及具体实施方式，进一步详述本发明的技术方案。

[0055] 本发明使用的触觉交互设备是Geomagic touch，该设备具有六自由度的输入与三自由度的反馈力输出。在本发明的流体模拟中，以应用CUDA框架的SPH多相流作为基础，并且加入与动态刚体的固/液耦合的方法和触觉渲染的方法和交互的方法，实现该方法的整体框架如图1所示。

[0056] 本发明的具体实施流程，包括以下步骤：

[0057] 1、实现SPH混合流模拟

[0058] 在流体模拟中，对于可混溶的液体，需要考虑速度场不均匀分布所带来的体积分数变化。本发明改进SPH混合流模拟的方法用于触觉交互，以实现混合流体搅拌的效果。

[0059] 本发明将各相流体混合后得到的混合流体模拟求解的问题转化为每个时间步长对于SPH粒子的各相流体积分数αk和混合速度um的求解：

[0060]

[0061]

[0062] 其中，pm为混合流体的压强，ρm为混合流体的密度，g为重力加速度，t为模拟的时间步长，▽·Tm为混合流体的粘性力，▽·TDm为混合流体中的各相流体之间由于漂移速度导致的动量变化，umk为混合流体中的各相流体的漂移速度，k为流体的相数，,um·▽表示混合速度与梯度算子相乘，▽表示梯度算子，▽pm表示混合流体压强的梯度。公式(1)与(2)分别转化自质量守恒与动量守恒。其中，公式(2)加入了触觉力项fhaptic，表示动态刚体对混合流体的作用力，该项在产生固/液耦合时才做计算。

[0063] 求解各相流体体积分数αk的关键在于对混合流体中的各相流体的漂移速度umk的求解，考虑到运算的效率以及前两项相对于扩散效果来说较为不明显，因此在计算漂移速度时只考虑了扩散效应的计算，本模拟方法中umk需要满足以下方程组：

[0064]

[0065] 公式(3)表示混合流体中每一点对应的各相流体漂移速度之和为零；公式(4)则是该漂移速度的求解公式；其中，σdif为各相流体的扩散系数，k'为其他相流体，▽αk表示体积分数的梯度，表示其他相流体的静态密度，在本发明中扩散系数取0.1有较好的效果。

[0066] 压力项▽pm的计算方式与参照基本的SPH方法中单相流压力项的计算方式[0067]

[0068] 其中，为混合流体中刚体粒子i的压强，刚体粒子j的压强。

[0069] 对各相流之间由于漂移速度导致的混合流体的动量变化▽·TDm的计算公式为：

[0070]

[0071] 其中，W(xji,h)为高斯核函数，h为高斯核函数的作用半径，xij＝(xi-xj)是SPH粒子i与刚体粒子j的距离，umki为混合流体中SPH粒子i的流体漂移速度，umkj为混合流体中刚体粒子j的流体漂移速度，：αki为流体SPH粒子i的各相流体积分数αk，αkj为刚体粒子j的各相流体积分数αk，mj表示其他粒子j的质量。对于粘性力的求解，由于各个流体粒子的粘性力随着各相流所占的百分比的变化而变化，流体粒子的粘性力由下式求解：

[0072]

[0073] 其中，μi为SPH粒子i的粘性系数，μj为刚体粒子j的粘性系数，umi为SPH粒子i混合后的速度，umj为刚体粒子j混合后的速度，W表示核函数，mj表示其他粒子j的质量。

[0074] SPH粒子与动态刚体交互所受到的触觉力fhaptic在本发明中可分解为fpres与fvis两项，fipres为刚体粒子所受到的压力，fivis为与刚体粒子交互的粘性力，在多相流模拟中一旦发生了流体与动态刚体的交互，则对其进行求解：

[0075]

[0076]

[0077] 在公式(7)、(8)中，i为SPH粒子，j为刚体粒子，pj为刚体粒子j的压强，Pmi为混合流体中刚体粒子i的压强，μj为刚体粒子j的粘性系数，μmi为混合流体中SPH粒子i的粘性系数，uj为刚体粒子的速度，umi为SPH粒子i混合后的速度，ρj为刚体粒子j的密度，Pmi为刚体粒子的压强，mj表示其他粒子j的质量，▽为梯度算子。在本发明的模拟步骤中，通过对上述公式求解，能够较好的模拟触觉交互时多相流的混合扩散效果，并且有较好的实时性，可以实现实时交互控制的任务。在该步骤中，本发明基于混合流模型，将应用于SPH方法的多相流模拟公式简化，保留最重要的扩散效应计算，在提高运行效率的同时保持逼真的多相流模拟现象，并应用触觉交互的方法。

[0078] 2、实现SPH流体的固/液耦合

[0079] 在流体模拟中，需要与流体发生交互作用的实体大致包括两类，静态实体和动态实体。相对于静态实体，动态刚体与流体粒子的耦合的问题相对来说更为复杂。本发明利用触觉设备来操控动态刚体与流体交互，并将反馈力实时传回给操控者。因此，反馈力的精确计算显得尤为重要。

[0080] 为保证交互与计算的精确性，对于该步骤，本发明使用统一粒子模型的方法，即：将刚体看作结构稳定的刚体粒子的结合，使用体素化的方法将刚体转化为集合粒子，并且使用高斯核函数来计算各种粒子之间的相互作用。

[0081] 将通过计算所有刚体粒子所受的合力Fhaptic作为触觉交互的反馈力。对刚体粒子压强的pr计算使用以下公式以防止渗流现象产生：

[0082]

[0083] 而流体粒子的压强pf也同样做这样的处理：

[0084]

[0085] 其中，ρ为动态密度，ρr和ρf为上述刚体粒子与流体的静态密度，kr与kf为压强系数；这里kf取1.5，kr取3.0。这两个公式使得压力的计算始终保持正值形成排斥力防止刚体与流体粒子相互渗流。产生的流体粒子进入刚体粒子，不符合实际。

[0086] 刚体粒子受到的触觉力公式表示为压力与粘性力：

[0087] fihaptic＝fipres+fivis       (11)

[0088] 其中，压力fipres与粘性力fivis的计算公式如下：

[0089]

[0090]

[0091] 最后计算合力，可以反馈给触觉设备作为反馈力输出：

[0092]

[0093] 在该步骤中，本发明使用体素化的方法将网格模型转化为刚体粒子的集合，使用统一粒子模型的方法对这些刚体粒子赋予初始属性和状态以后让其与流体粒子同样进行状态更新，并在与流体交互的过程中保持刚体结构的稳定性且防止渗流效果；并计算所有刚体粒子所受的合力作为触觉交互的反馈力。

[0094] 3、实现触觉交互操控与触觉力渲染

[0095] 动态刚体的运动可以看作是刚体围绕质心进行平移和旋转，并且不会产生形变。本发明中操作触觉设备控制虚拟场景中粒子表示的动态刚体运动可以转化为另一种思路：
将刚体的位移与旋转表现为每一个刚体粒子的空间位置相对质心的修正过程。

[0096] 假设所有粒子的质量相同，对于一个刚体粒子，速度公式可以表示为：

[0097] ui＝ur+ωr×qi      (15)

[0098] 其中，ur是刚体质心的速度，ωr是刚体的角速度，ui是单个刚体粒子j的速度，qi是粒子ij与质心的相对位置。在程序框架中通过每一个时间步长更新交互设备的位置与旋转角度的输入参数来实时的控制虚拟场景中动态刚体的运动。

[0099] 在触觉力渲染部分，本发明将触觉力的更新从模拟步骤分离开并且通过滤波来填充模拟步骤之间过渡触觉力实现平滑模拟的效果。本发明使用带汉明窗的FIR(Finite Impulse Response)滤波器将离散的低频率的触觉力数据平滑化为1000赫兹的高频率信息，这样可以防止使用者在操作过程中感受到震动。虽然FIR滤波器的特性会带来一定的反馈力延迟，但是由于延迟时间很短，在实际使用时很难发现。经过实验发现，使用滤波的方法将触觉力平滑化效果好于未平滑的初始效果。

[0100] 本发明实验机器配置如下：Intel Core i7 4700mq，GTX765m，16G RAM，Windows 8.1，CUDA 7.0，VS2013。并支持Windows 7，Windows8.1，Windows 10环境，1G内存以上，nVidia GT430以上显卡配置的微机。