CN103869983B - 一种用于力触觉人机交互的柔性物体变形仿真方法 - Google Patents

Abstract

本发明提出一种用于力触觉人机交互的柔性物体变形仿真方法，基于等边三角形分割的弹簧质点模型，将柔性物体的变形视为由于作用力从局部区域内的表面质点逐层地、顺序地传播到组织内的各个质点，直至到达某预设的边界为止；并不需要对变形物体作整体考虑，而是有效地将计算资源集中投入到局部区域里的质点。本方法的计算只涉及柔性物体表面的局部区域内的质点，变形区域由接触力决定，分割精度可以根据准确性要求调整，不仅大大简化了计算量，加快了变形计算的速度，且调整灵活，适用范围广泛，可应用于具有各向同性的柔性物体在平稳接触下的变形情况，例如虚拟外科手术仿真、遥操作机器人控制等虚拟现实人机交互领域。

Description

一种用于力触觉人机交互的柔性物体变形仿真方法

技术领域

本发明属于柔性触觉交互技术领域，涉及一种柔性物体变形仿真方法，尤其是涉及一种用于力触觉人机交互的，基于等边三角形分割的弹簧质点模型的柔性物体变形仿真方法。

背景技术

柔性触觉交互技术是虚拟现实技术在触觉再现上的一个重要应用，通过计算机构成虚拟的三维空间，将现实环境复制到计算机中去产生逼真的“虚拟环境”，从而使得用户在多种感官上产生一种沉浸于虚拟环境的感觉。某些专用的交互设备可以将人的操作动作转换成为虚拟现实系统能够识别的指令序列，同时该设备可以使人能够感受到虚拟环境提供的反馈。当人与虚拟环境交互时，将图形显示和触觉显示结合起来，使操作者不但能通过眼睛观察到虚拟场景的视觉信息，还可以通过手或者相关工具对真实环境进行操作，而且还可以感受得到来自真实环境的反作用力，可以极大增强操作者身临其境的感受。

在柔性触觉交互系统中，操作者通过该系统可以触摸虚拟环境中的虚拟物体，感知虚拟物体的柔性、刚度、表面纹理等触觉特性。该技术的应用范围非常广泛，目前在娱乐、教育、艺术、军事、航空、医学、机器人等多方面有着相当广泛的应用。

模拟物体变形技术主要有两个方面，基于几何模型和基于物理模型技术。由于基于几何模型方法很难模拟复杂的物体受力变形，所以基于物理模型的模拟物体变形技术受到越来越多的研究者的关注。目前应用于物体的物理变形模型方法主要有质点-弹簧模型和有限单元模型。

在经典的质点-弹簧模型中，整个物体被模拟成一个由大量质点和弹簧组成的系统，一般每个质点都和其邻域中的至少8个邻点用弹簧相连。进行变形计算时，要对系统中所有的质点建立动力学偏微分方程。对这些偏微分方程，一般采用经典的欧拉方法、龙格-库塔或共轭梯度法求解。求解过程需要迭代多次才能收敛到平衡状态。尽管和有限元方法相比计算量小的多，但是这种方法也存在如下缺点：

首先，由于模型是由不连续的拉格朗日方程得来，这个决定系统动态特性的微分方程必须满足一些特殊条件，才能避免在求数值解时失稳。

其次，存在力的到达深度不受限制的问题。力的到达深度是指在质点一弹簧模型中，由应力引起的变形过程：力首先由中心质点通过弹簧传播到邻近的质点，再由这些邻近的质点进一步传播至更远的和他们相邻的质点。一般的算法没有对弹簧的弹性力的传播深度做限制，即认为接触力影响的范围是所有的质点，因而要对所有的质点建立拉格朗日力学方程。但是实际作用力能够影响受力变形范围是有限的。所以对于那些势能无限趋近于零的弹簧质点就没有计算的必要。

再者，由于网格的拓扑结构很大程度上决定了系统的行为，若拓扑结构设计不合理就会导致失真。具体建立某个物体模型时，需要根据物体的材料来确定所采用的弹簧参数。但这些参数一般很难从物体材料的特性参数中直接得到，如果参数值选取不合适就会导致较大的误差，而且一些材料的属性也无法在模型中得到自然的表达；要达到相同的计算精度，弹簧质点模型需要比其它模型更细密的网格分割，这使得计算量大大增加，目前的趋势是采用细节层次技术来减少要处理的质点数，但不同层次的模型还存在动态衔接问题；在动态弹簧质点系统中，刚度K与计算所需的时间步长成反比，即若物体的刚度大则计算所需的时间步长就要减小。

发明内容

针对弹簧质点模型存在缺乏力的作用深度限制以及精确性和实时性难于控制的问题，本发明提出一种用于力触觉人机交互的柔性物体变形仿真方法，该仿真方法基于等边三角形分割的弹簧质点模型，能准确快速的计算变形量和力反馈，实现对柔性物体的实时变形仿真，从而提高虚拟力触觉交互的逼真度。

为了达到上述目的，本发明提供如下技术方案：

一种用于力触觉人机交互的柔性物体变形仿真方法，包括如下步骤：

步骤1对虚拟场景进行初始化；

步骤2检测虚拟手和虚拟物体是否接触，当检测到虚拟手和虚拟物体接触式，从力触觉人机交互装置获得接触力；

步骤3根据步骤2中获得的接触力，基于等边三角形分割的弹簧质点模型计算虚拟物体的变形和力反馈，所述基于等边三角形分割的弹簧质点模型建立方法为：

（1）以接触点为中心，把整个柔性物体分割成一系列均匀的同心等边三角形，边长依次为a，2a，3a……；每圈均匀分布有6个离散的由弹簧相连的质点；

（2）假设对x正方向的弹簧质点，f_z是作用力f在垂直方向上作用力分量，对于各向同性的情况，f_z＝f/6；

在作用点中心，根据受力平衡关系有：

f_z＝k_N1Δz₁+k_T1Δa_isinθ₁＝f_N1+f_T1sinθ₁ （1）

对于第n个质点有：

$\left\{\begin{array}{c}{f}_{T_i}\cos {\mathrm{\θ}}_i=f_{T_{(i-1)}}\cos {\mathrm{\θ}}_{i-1}\\ f_{N_i}=k_{N_i}{\mathrm{\Δz}}_i=f_{T_{i-1}}\sin {\mathrm{\θ}}_{i-1}-f_{T_i}\sin {\mathrm{\θ}}_i\end{array}\right.i=\mathrm{2,3},...,m,...---\left(2\right)$

其中表示切向力，表示法向力，表示第N_i个质点法线方向上弹簧的弹性模量，Δa_i表示切向上第i个质点与第i+1个质点间弹簧的变形量，△z_i表示第i个质点在法向上的变形量，θ_i为第i个质点与第i+1个质点间连线与水平方向的夹角；

假设作用力和形变是缓慢均匀的，得到

$f_{T_i}\cos {\mathrm{\θ}}_1=f_{T_2}\cos {\mathrm{\θ}}_2=...=f_{T_{i-1}}\cos {\mathrm{\θ}}_{i-1}=f_{T_i}\cos {\mathrm{\θ}}_i=...\left(3\right)$

由式(2)、式(3)推出：

$f_{N_1}=f_z-f_{T_1}\sin {\mathrm{\θ}}_1$

$f_{N_2}=f_{T_1}\sin {\mathrm{\θ}}_1-f_{T_1}\sin {\mathrm{\θ}}_2$

……

$f_{N_{i-1}}=f_{T_{i-1}}\sin {\mathrm{\θ}}_{i-1}-f_{T_{i-1}}\sin {\mathrm{\θ}}_{i-1}$

$f_{N_i}=f_{T_{i-1}}\sin {\mathrm{\θ}}_{i-1}-f_{T_i}\sin {\mathrm{\θ}}_i$

合并以上方程有：

$f_{N_i}=f_{T_{i-1}}\sin {\mathrm{\θ}}_{i-1}-f_{T_i}\sin {\mathrm{\θ}}_i---\left(4\right)$

设当i=m，sinθ_m≈0，则

$f_{N_1}+f_{N_2}+...+f_{N_{m-1}}+f_{N_m}=f_z---\left(5\right)$

$k_{N_1}{\mathrm{\Δz}}_1+k_{N_2}{\mathrm{\Δz}}_2+...+k_{N_{m-1}}{\mathrm{\Δz}}_{m-1}+k_{N_m}{\mathrm{\Δz}}_m=f_z---\left(6\right)$

由弹性力学可知，根据等效关系，法线方向上第i个弹簧的弹性模量正比于所对应的三角形或梯形的面积s_i，表示如下：

$k_{N_i}\∝s_i=\frac{\sqrt{3}}{4}\{{\left(\mathrm{ia}\right)}^2-{\left[(i-1)a\right]}^2\}=\frac{\sqrt{3}}{4}(2i-1)a^2---\left(7\right)$

由上式可得到：

$k_{N_1}:k_{N_2}:k_{N_3}:...:k_{N_m}=1:3:5:...:(2m-1)---\left(8\right)$

切线方向上相邻的弹簧弹性模量之比为常数，即第i个切向弹簧的弹性模量和第一个切向弹簧弹性模量是指数倍关系，表示如下：

其中i＝2,3...m，β为系数,λ_i≥1 （9）

则切线方向上形变量为：

${\mathrm{\Δa}}_i=\frac{f_{T_i}}{k_{T_i}}=\frac{f_{T_1}\cos {\mathrm{\θ}}_1}{\cos {\mathrm{\θ}}_i}\×\frac{1}{{\mathrm{\λ}}_i{k}_{T_1}}=\frac{1}{{\mathrm{\λ}}_i}\frac{\cos {\mathrm{\θ}}_1}{\cos {\mathrm{\θ}}_i}{\mathrm{\Δa}}_1,(i=\mathrm{1,2}...m)---\left(10\right)$

$\frac{{\mathrm{\Δa}}_i}{{\mathrm{\Δa}}_1}=\frac{1}{{\mathrm{\λ}}_i}\frac{\cos {\mathrm{\θ}}_1}{\cos {\mathrm{\θ}}_i}=\frac{1}{{\mathrm{\λ}}_i}\frac{\frac{a}{a+{\mathrm{\Δa}}_1}}{\frac{a}{a+{\mathrm{\Δa}}_i}}=\frac{1}{{\mathrm{\λ}}_i}\frac{1+\frac{{\mathrm{\Δa}}_i}{a}}{1+\frac{{\mathrm{\Δa}}_1}{a}}---\left(11\right)$

$\frac{\frac{{\mathrm{\Δa}}_i}{a}}{\frac{{\mathrm{\Δa}}_1}{a}}=\frac{1}{{\mathrm{\λ}}_i}\frac{1+\frac{{\mathrm{\Δa}}_i}{a}}{1+\frac{{\mathrm{\Δa}}_1}{a}}---\left(12\right)$

令是第i个切向弹簧的变形率，则有

${\mathrm{\ϵ}}_i=\frac{{\mathrm{\ϵ}}_1}{{\mathrm{\λ}}_i+({\mathrm{\λ}}_i-1){\mathrm{\ϵ}}_1}---\left(13\right)$

法线方向上第i和第i+1个点之间的相对形变：

$l_i=\sqrt{{(a+{\mathrm{\Δa}}_i)}^2-a^2}=\sqrt{({2\mathrm{\ϵ}}_i+{\mathrm{\ϵ}}_i^2)}a---\left(14\right)$

将公式（13）代入（14），可得：

$\frac{l_i}{a}=\frac{1}{{\mathrm{\λ}}_i+({\mathrm{\λ}}_i-1){\mathrm{\ϵ}}_1}\sqrt{[{2\mathrm{\λ}}_i+(2{\mathrm{\λ}}_1-1){\mathrm{\ϵ}}_1]{\mathrm{\ϵ}}_1}---\left(15\right)$

当i＝1时，有： $l_1=a\sqrt{(2+{\mathrm{\ϵ}}_1){\mathrm{\ϵ}}_1}$

所以： $\frac{l_i}{l_1}=\frac{1}{(1+{\mathrm{\ϵ}}_i){\mathrm{\λ}}_i-{\mathrm{\ϵ}}_1}\sqrt{\frac{2(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_i-{\mathrm{\ϵ}}_1}{2+{\mathrm{\ϵ}}_1}}---\left(16\right)$

令ε_i＜0.5，当i≥2时有λ_i＞3，则

(1+ε_i)λ_i>>ε_i

则： $\frac{l_i}{l_1}\≈\sqrt{\frac{2}{(1+{\mathrm{\ϵ}}_1)(2+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_i}}=\frac{1}{\mathrm{\ρ}}\·\frac{1}{\sqrt{{\mathrm{\λ}}_i}},i=\mathrm{2,3}....,m---\left(17\right)$ 其中 $\mathrm{\η}=\sqrt{(1+{\mathrm{\ϵ}}_1)(1+\frac{{\mathrm{\ϵ}}_1}{2})},$ 所以有以下重要关系式：

$l_i=\frac{1}{\mathrm{\ρ}}\·\frac{1}{\sqrt{{\mathrm{\λ}}_i}}\·L_1,i=\mathrm{2,3}....,m---\left(18\right)$

每个质点的总变形量为：

$\left\{\begin{array}{c}{\mathrm{\Δz}}_1=l_1+l_2+l_3+...+l_m\\ {\mathrm{\Δz}}_2=l_2+l_3+...+l_m\\ ...\\ {\mathrm{\Δz}}_m=l_m\end{array}\right.---\left(19\right)$

将式（19）、（9）代入式（6）有

$F_z=k_{N_1}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}(2i-1){\mathrm{\Δz}}_i---\left(20\right)$

$F_z=k_{N_1}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{i}^2{l}_i---\left(21\right)$

其中，m表示在虚拟物体上以受力点为中心，划分的同心的等边三角形的个数，由方程（15）得到

$\frac{l_m}{a}=\frac{1}{(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_m-{\mathrm{\ϵ}}_1}\sqrt{[2(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_m-{\mathrm{\ϵ}}_1]{\mathrm{\ϵ}}_1}---\left(22\right)$

因为(1+ε₁)λ_m＞＞ε₁，所以

$\frac{l_m}{a}\≈\frac{1}{(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_m}\sqrt{\left[2(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_m\right]{\mathrm{\ϵ}}_1}=\sqrt{\frac{{2\mathrm{\ϵ}}_1}{(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_m}}---\left(23\right)$

由上式得到： $m=\mathrm{\β}\ln \left(\frac{{2\mathrm{\ϵ}}_1}{1+{\mathrm{\ϵ}}_1}{\left(\frac{a}{l_m}\right)}^2\right)+1---\left(24\right)$

其中，l_m/a代表了计算区域边界的变形量，而ε₁表示受力中心点的切向弹簧的变形率；

（3）计算虚拟物体发生变形的区域，获得发生变形的弹簧层数，计算该局部区域内各层弹簧变形量之和为虚拟物体的变形，各层弹簧弹力的合力为力反馈。

优选的，所述β取值为4.0，受力中心点的切向弹簧的变形率ε₁取值为0.3，l_m/a取值为1/57.3。

有益效果：

与现有技术相比，本发明提供的柔性物体变形仿真方法，在虚拟现实环境中的虚拟手和虚拟柔性物体的接触过程中，将柔性物体的变形视为由于作用力从局部区域内的表面质点逐层地、顺序地传播到组织内的各个质点，直至到达某预设的边界为止；基于给定的接触力，计算将发生变形的局部区域，并不需要对变形物体作整体考虑，而是有效地将计算资源集中投入到局部区域里的质点；输出的力反馈信号等于逐层弹簧弹性力的叠加。本方法的计算只涉及柔性物体表面的局部区域内的质点，变形区域由接触力决定，分割精度可以根据准确性要求调整，不仅大大简化了计算量，加快了变形计算的速度，且调整灵活，在力触觉人机交互过程中，自然舒适、力触觉感觉平稳、模拟效果逼真。本发明适用范围广泛，可应用于具有各向同性的柔性物体在平稳接触下的变形情况，例如虚拟外科手术仿真、遥操作机器人控制等虚拟现实人机交互领域。

附图说明

图1为本发明提供的柔性变形仿真流程图；

图2为基于等边三角形分割的弹簧-质点形变模型；

图3为法向弹簧弹性模量与半径关系图；

图4为相邻质点法向变形位置示意图；

图5为虚拟长方体受压后受力变形渲染效果图。

具体实施方式

以下将结合具体实施例对本发明提供的技术方案进行详细说明，应理解下述具体实施方式仅用于说明本发明而不用于限制本发明的范围。

本例采用Phantom Omini手控器作为力触觉人机交互装置，以6058个质点，11328个等边三角网格构成的各向同性的长方体（如图5所示）为例来进行仿真模拟，各向同性的虚拟长方体进行按压变形仿真，并将交互过程中产生的力触觉信息实时反馈给操作者，软件用VC++6.0、MFC、OpenGL编程，在Pentium42.80GHz，1GB内存，NVIDIA GeForce2显卡的PC机上实现了变形仿真。具体地说，仿真过程如图1所示，包括如下步骤：

步骤1力触觉人机交互装置Phantom Omini手控器和虚拟场景初始化；

步骤2当检测到虚拟手和虚拟物体发生接触时，在接触力作用下，计算虚拟物体发生形变的区域，在该区域内基于等边三角形分割的弹簧质点模型进行仿真计算，将柔性物体的变形视为由于作用力从局部区域内的表面质点逐层地、顺序地传播到组织内的各个质点，直至到达预设的边界为止。系统中各个质点的相对位移的叠加对外等效为物体表面变形，与之相连的弹簧弹性力的合力等于物体表面的接触力；在人机交互过程中输出力反馈信号，等效于各层弹簧弹性力的合力，具体操作如下：

（1）以接触点为中心，把整个柔性物体分割成一系列均匀的同心等边三角形，边长依次为a，2a，3a……；每圈均匀分布有6个离散的由弹簧相连的质点，如图2所示。

以上模型的主要优点在于使用三角形面片来构造虚拟物体表面模型是最常用的方法，能更好的和虚拟物体的几何模型相结合；物理意义和几何建模的结合也更加紧密，能确保两者在受力变化时的一致性。选择等边三角形的原因是以一个顶点为中心的六个等边三角形是中心对称的，而且其构成的正六边形扩展后是类似蜂巢格的单元结构，可以和任意多这样的单元结构拼合扩展构图；且对于任意顶点来说周围的图形都是对称的，所以可以简化对虚拟物体变形的分析。

假设在未受力时，物体表面和内部的质点由弹簧相连，均匀分布。先考虑物体的形变各向同性的情况，即在水平平面上的所有的弹簧—质点都具有相同的特性。整个变形是各个分割单元变形迭代累加的结果。

（2）由于受力变形情况以受力点为中心对称分布，x方向和y方向的变形相同，为简化分析过程，可以只分析x正方向的受力和形变情况。假设对x正方向的弹簧质点，f_z是作用力f在垂直方向上作用力分量，对于各向同性的情况，f_z＝f/6。

在作用点中心，根据受力平衡关系有：

f_z＝k_N1Δz₁+k_T1Δa_isinθ₁＝f_N1+f_T1sinθ₁ （1）

对于第n个质点有：

$\left\{\begin{array}{c}{f}_{T_i}\cos {\mathrm{\θ}}_i=f_{T_{(i-1)}}\cos {\mathrm{\θ}}_{i-1}\\ f_{N_i}=k_{N_i}{\mathrm{\Δz}}_i=f_{T_{i-1}}\sin {\mathrm{\θ}}_{i-1}-f_{T_i}\sin {\mathrm{\θ}}_i\end{array}\right.i=\mathrm{2,3},...,m,...---\left(2\right)$

其中表示切向力，表示法向力，表示第N_i个质点法线方向上弹簧的弹性模量，Δa_i表示切向上第i个质点与第i+1个质点间弹簧的变形量，△z_i表示第i个质点在法向上的变形量，θ_i为第i个质点与第i+1个质点间连线与水平方向的夹角。各个方向上的作用力和形变是均匀的，即物体形变是各向同性的。

$f_{T_i}\cos {\mathrm{\θ}}_1=f_{T_2}\cos {\mathrm{\θ}}_2=...=f_{T_{i-1}}\cos {\mathrm{\θ}}_{i-1}=f_{T_i}\cos {\mathrm{\θ}}_i=...\left(3\right)$

由式(2)、式(3)推出：

$f_{N_1}=f_z-f_{T_1}\sin {\mathrm{\θ}}_1$

$f_{N_2}=f_{T_1}\sin {\mathrm{\θ}}_1-f_{T_1}\sin {\mathrm{\θ}}_2$

……

$f_{N_{i-1}}=f_{T_{i-1}}\sin {\mathrm{\θ}}_{i-1}-f_{T_{i-1}}\sin {\mathrm{\θ}}_{i-1}$

$f_{N_i}=f_{T_{i-1}}\sin {\mathrm{\θ}}_{i-1}-f_{T_i}\sin {\mathrm{\θ}}_i$

合并以上方程有：

$f_{N_1}+f_{N_2}+...+f_{N_{i-1}}+f_{N_i}=f_z-f_{T_i}\sin {\mathrm{\θ}}_1---\left(4\right)$

在sinθ_m≈0时，令i=m，m是发生形变的三角形面片的边界，则

$f_{N_1}+f_{N_2}+...+f_{N_{m-1}}+f_{N_m}=f_z---\left(5\right)$

$k_{N_1}{\mathrm{\Δz}}_1+k_{N_2}{\mathrm{\Δz}}_2+...+k_{N_{m-1}}{\mathrm{\Δz}}_{m-1}+k_{N_m}{\mathrm{\Δz}}_m=f_z---\left(6\right)$

由弹性力学可知，根据等效关系，法线方向上第i个弹簧的弹性模量正比于所对应的三角形或梯形的面积s_i，如图3所示。

$k_{N_i}\∝s_i=\frac{\sqrt{3}}{4}\{{\left(\mathrm{ia}\right)}^2-{\left[(i-1)a\right]}^2\}=\frac{\sqrt{3}}{4}(2i-1)a^2---\left(7\right)$

可得到：

$k_{N_1}:k_{N_2}:k_{N_3}:...:k_{N_m}=1:3:5:...:(2m-1)---\left(8\right)$

切线方向上相邻的弹簧弹性模量之比为常数，即第i个切向弹簧的弹性模量和第一个切向弹簧弹性模量是指数倍关系：

其中i＝2,3...m，β为系数,λ_i≥1 （9）

切线方向上形变量为：

${\mathrm{\Δa}}_i=\frac{f_{T_i}}{k_{T_i}}=\frac{f_{T_1}\cos {\mathrm{\θ}}_1}{\cos {\mathrm{\θ}}_i}\×\frac{1}{{\mathrm{\λ}}_i{k}_{T_1}}=\frac{1}{{\mathrm{\λ}}_i}\frac{\cos {\mathrm{\θ}}_1}{\cos {\mathrm{\θ}}_i}{\mathrm{\Δa}}_1,(i=\mathrm{1,2}...m)---\left(10\right)$

$\frac{{\mathrm{\Δa}}_i}{{\mathrm{\Δa}}_1}=\frac{1}{{\mathrm{\λ}}_i}\frac{\cos {\mathrm{\θ}}_1}{\cos {\mathrm{\θ}}_i}=\frac{1}{{\mathrm{\λ}}_i}\frac{\frac{a}{a+{\mathrm{\Δa}}_1}}{\frac{a}{a+{\mathrm{\Δa}}_i}}=\frac{1}{{\mathrm{\λ}}_i}\frac{1+\frac{{\mathrm{\Δa}}_i}{a}}{1+\frac{{\mathrm{\Δa}}_1}{a}}---\left(11\right)$

$\frac{\frac{{\mathrm{\Δa}}_i}{a}}{\frac{{\mathrm{\Δa}}_1}{a}}=\frac{1}{{\mathrm{\λ}}_i}\frac{1+\frac{{\mathrm{\Δa}}_i}{a}}{1+\frac{{\mathrm{\Δa}}_1}{a}}---\left(12\right)$

令是第i个切向弹簧的变形率，则有

${\mathrm{\ϵ}}_i=\frac{{\mathrm{\ϵ}}_1}{{\mathrm{\λ}}_i+({\mathrm{\λ}}_i-1){\mathrm{\ϵ}}_1}---\left(13\right)$

如图4所示，法线方向上第i和第i+1个点之间的相对形变为：

$l_i=\sqrt{{(a+{\mathrm{\Δa}}_i)}^2-a^2}=\sqrt{({2\mathrm{\ϵ}}_i+{\mathrm{\ϵ}}_i^2)}a---\left(14\right)$

将公式（13）代入（14），可得：

$\frac{l_i}{a}=\frac{1}{{\mathrm{\λ}}_i+({\mathrm{\λ}}_i-1){\mathrm{\ϵ}}_1}\sqrt{[{2\mathrm{\λ}}_i+(2{\mathrm{\λ}}_1-1){\mathrm{\ϵ}}_1]{\mathrm{\ϵ}}_1}---\left(15\right)$

当i＝1时，有： $l_1=a\sqrt{(2+{\mathrm{\ϵ}}_1){\mathrm{\ϵ}}_1}$

所以： $\frac{l_i}{l_1}=\frac{1}{(1+{\mathrm{\ϵ}}_i){\mathrm{\λ}}_i-{\mathrm{\ϵ}}_1}\sqrt{\frac{2(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_i-{\mathrm{\ϵ}}_1}{2+{\mathrm{\ϵ}}_1}}---\left(16\right)$

一般情况下：ε_i＜0.5，当i≥2时有λ_i＞3，则

(1+ε_i)λ_i＞＞ε_i

则： $\frac{l_i}{l_1}\≈\sqrt{\frac{2}{(1+{\mathrm{\ϵ}}_1)(2+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_i}}=\frac{1}{\mathrm{\ρ}}\·\frac{1}{\sqrt{{\mathrm{\λ}}_i}},i=\mathrm{2,3}....,m---\left(17\right)$

其中所以有以下重要关系式

$l_i=\frac{1}{\mathrm{\ρ}}\·\frac{1}{\sqrt{{\mathrm{\λ}}_i}}\·L_1,i=\mathrm{2,3}....,m---\left(18\right)$

每个质点的总变形量为：

$\left\{\begin{array}{c}{\mathrm{\Δz}}_1=l_1+l_2+l_3+...+l_m\\ {\mathrm{\Δz}}_2=l_2+l_3+...+l_m\\ ...\\ {\mathrm{\Δz}}_m=l_m\end{array}\right.---\left(19\right)$

将式（19）、（9）代入式（6）有

$F_z=k_{N_1}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}(2i-1){\mathrm{\Δz}}_i---\left(20\right)$

$F_z=k_{N_1}\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{i}^2{l}_i---\left(21\right)$

这里，m是模型的一个重要边界参数，它表示在虚拟物体上以受力点为中心，划分的同心的等边三角形的个数，由方程（15）得到

$\frac{l_m}{a}=\frac{1}{(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_m-{\mathrm{\ϵ}}_1}\sqrt{[2(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_m-{\mathrm{\ϵ}}_1]{\mathrm{\ϵ}}_1}---\left(22\right)$

因为(1+ε₁)λ_m＞＞ε₁，所以

$\frac{l_m}{a}\≈\frac{1}{(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_m}\sqrt{\left[2(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_m\right]{\mathrm{\ϵ}}_1}=\sqrt{\frac{{2\mathrm{\ϵ}}_1}{(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_m}}---\left(23\right)$

$m=\mathrm{\β}\ln \left(\frac{{2\mathrm{\ϵ}}_1}{1+{\mathrm{\ϵ}}_1}{\left(\frac{a}{l_m}\right)}^2\right)+1---\left(24\right)$

其中，l_m/a代表了计算区域边界的变形量，是边界参数，而ε₁则代表了受力中心点的切向弹簧的变形率。

在上述模型中，计算量的大小取决于同心的等边的三角形的分割数目m的大小，由式（22）和（24）可知，m的大小取决于边界条件l_ma和接触点的变形率ε₁，因此整个模型的计算量就取决于ε₁和边界参数l_m/a的选取：当边界条件l_m/a固定时，ε₁增大，计算点数m迅速增加，计算量会随之增大。当a固定时，不同的物体，受相同的力时，ε₁越大表示物体的越柔软；同一物体，受力越大，受力点的变形率ε₁越大。而变形区域大小由参数m和a决定。参数m取决于变形率ε₁，而ε₁又和受力大小有关。对于同一种材料，若a固定，力越大，计算量就越大，变形区域也会相应增大，但变形区域也不会无限扩大，随着力的增大，变形区域增大的速度越来越缓慢。

经实验验证，本例中取β=4.0，受力中心点的切向弹簧的变形率ε₁=0.3，在上述取值下变形效果较好。

考虑到计算精度，取l_m/a=1/57.3，即模型边界上的弹簧变形角度为1°，认为边界外无变形量。

计算结果按四舍五入取整，得到m=8。

即计算8层弹簧质点模型的变形量和力反馈，计算结果按四舍五入取小数点后四位。

受力中心点的切向弹簧的变形率ε₁=0.3，则当i＝1时，有：

$l_1=a\sqrt{(2+{\mathrm{\ϵ}}_1){\mathrm{\ϵ}}_1},$ 即得到l₁=0.8307a

l_{i（i=2,3,…m）}和l₁的关系为：

$l_i=\frac{1}{\mathrm{\ρ}}\·\frac{1}{\sqrt{{\mathrm{\λ}}_i}}\·l_1,$ i＝2,3....,m，其中 $\mathrm{\ρ}=\sqrt{(1+{\mathrm{\ϵ}}_1)(1+\frac{{\mathrm{\ϵ}}_1}{2})}$

由于切线方向上相邻的弹簧弹性模量之比为常数，即第i个切向弹簧的弹性模量和第一个切向弹簧弹性模量是指数倍关系，即：

其中i＝2,3...m，β为系数,λ_i≥1

取e=2.71828，β=4.0，λ₂=1.2840，

计算出：l₂=0.5996a

重复上述过程，可以计算出l₃…l_m。

从而得到每个质点的变形量ΔZ₁，ΔZ₂…ΔZ_m。

计算虚拟物体发生变形的区域，获得发生变形的弹簧层数m，该局部区域内各层弹簧变形量之和为虚拟物体的变形，各层弹簧弹力的合力为力反馈。

各层弹簧的弹力的合力为力反馈：

$f_z=f_{N_1}+f_{N_2}+...+f_{N_{m-1}}+f_{N_m}$

每层弹簧的弹力可以根据弹性模量和变形量计算，从而得到力反馈如下：

$f_z=k_{N_1}{\mathrm{\Δz}}_1+k_{N_1}+k_{N_2}{\mathrm{\Δz}}_2+...+k_{N_{m-1}}{\mathrm{\Δz}}_{m-1}+k_{N_m}{\mathrm{\Δz}}_m$

步骤3，根据步骤2中计算出的虚拟物体的变形更新绘制变形后的虚拟物体（如图5所示，），并将力反馈信息实时反馈给操作者。

根据上述步骤，本方法完成一次变形计算的时间开销为30ms，1秒内变形计算的刷新率达33次，变形效果如图5所示。实验结果表明：本发明提供的建模方法不仅计算简单，而且能够保证触觉接触力和变形计算具有较高精度，模拟效果逼真。在交互过程中，操作者可以实时、真实地感知到变形仿真过程中虚拟手与虚拟长方体之间的力触觉信息，交互自然，舒适，力触觉感觉平稳。

本发明方案所公开的技术手段不仅限于上述实施方式所公开的技术手段，还包括由以上技术特征任意组合所组成的技术方案。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也视为本发明的保护范围。

Claims (2)

1.一种用于力触觉人机交互的柔性物体变形仿真方法，其特征在于，包括如下步骤：

步骤1对虚拟场景进行初始化；

步骤2检测虚拟手和虚拟物体是否接触，当检测到虚拟手和虚拟物体接触式，从力触觉人机交互装置获得接触力；

步骤3根据步骤2中获得的接触力，基于等边三角形分割的弹簧质点模型计算虚拟物体的变形和力反馈，所述基于等边三角形分割的弹簧质点模型建立方法为：

(1)以接触点为中心，把整个柔性物体分割成一系列均匀的同心等边三角形，边长依次为a，2a，3a……；每圈均匀分布有6个离散的由弹簧相连的质点；

(2)假设对x正方向的弹簧质点，f_z是作用力f在垂直方向上作用力分量，对于各向同性的情况，f_z＝f/6；

在作用点中心，根据受力平衡关系有：

f_z＝k_N1Δz₁+k_T1Δa_isinθ₁＝f_N1+f_T1sinθ₁ (1)

对于第n个质点有：

$\left\{\begin{array}{c}{f}_{T_i}{\mathrm{cos\θ}}_i=f_{T_{\left(i-1\right)}}{\mathrm{cos\θ}}_{i-1}\\ f_{N_i}=k_{N_i}{\mathrm{\Δz}}_i=f_{T_{i-1}}{\mathrm{sin\θ}}_{i-1}-f_{T_i}{\mathrm{sin\θ}}_i\end{array}\right.,i=2,3,...,m,...---\left(2\right)$

其中表示切向力，表示法向力，表示第N_i个质点法线方向上弹簧的弹性模量，Δa_i表示切向上第i个质点与第i+1个质点间弹簧的变形量，△z_i表示第i个质点在法向上的变形量，θ_i为第i个质点与第i+1个质点间连线与水平方向的夹角；

假设作用力和形变是缓慢均匀的，得到

$f_{T_1}{\mathrm{cos\θ}}_1=f_{T_2}{\mathrm{cos\θ}}_2=...=f_{T_{i-1}}{\mathrm{cos\θ}}_{i-1}=f_{T_i}{\mathrm{cos\θ}}_i=...---\left(3\right)$

由式(2)、式(3)推出：

$f_{N_1}=f_z-f_{T_1}{\mathrm{sin\θ}}_1$

$f_{N_2}=f_{T_1}{\mathrm{sin\θ}}_1-f_{T_1}{\mathrm{sin\θ}}_2$

……

$f_{N_{i-1}}=f_{T_{i-2}}{\mathrm{sin\θ}}_{i-2}-f_{T_{i-1}}{\mathrm{sin\θ}}_{i-1}$

$f_{N_i}=f_{T_{i-1}}{\mathrm{sin\θ}}_{i-1}-f_{T_i}{\mathrm{sin\θ}}_i$

合并以上方程有：

$f_{N_1}+f_{N_2}+...+f_{N_{i-1}}+f_{N_i}=f_z-f_{T_i}{\mathrm{sin\θ}}_i---\left(4\right)$

设当i＝m，sinθ_m≈0，则

$f_{N_1}+f_{N_2}+...+f_{N_{m-1}}+f_{N_m}=f_z---\left(5\right)$

$k_{N_1}{\mathrm{\Δz}}_1+k_{N_2}{\mathrm{\Δz}}_2+...+k_{N_{m-1}}{\mathrm{\Δz}}_{m-1}+k_{Nm}{\mathrm{\Δz}}_m=f_z---\left(6\right)$

由弹性力学可知，根据等效关系，法线方向上第i个弹簧的弹性模量k_Ni正比于所对应的三角形或梯形的面积s_i，表示如下：

$k_{N_i}\∝s_i=\frac{\sqrt{3}}{4}\{{\left(ia\right)}^2-{\[(i-1)a\]}^2\}=\frac{\sqrt{3}}{4}(2i-1)a^2---\left(7\right)$

由上式可得到：

$k_{N_1}:k_{N_2}:k_{N_3}:...:k_{N_m}=1:3:5:...:(2m-1)---\left(8\right)$

切线方向上相邻的弹簧弹性模量之比为常数，即第i个切向弹簧的弹性模量和第一个切向弹簧弹性模量是指数倍关系，表示如下：

其中i＝2,3...m，β为系数,λ_i≥1 (9)

则切线方向上形变量为：

${\mathrm{\Δa}}_i=\frac{f_{T_i}}{k_{T_i}}=\frac{f_{T_1}{\mathrm{cos\θ}}_1}{{\mathrm{cos\θ}}_i}\×\frac{1}{{\mathrm{\λ}}_i{k}_{T_1}}=\frac{1}{{\mathrm{\λ}}_i}\frac{{\mathrm{cos\θ}}_1}{{\mathrm{cos\θ}}_i}{\mathrm{\Δa}}_1,\left(i=1,2...m\right)---\left(10\right)$

$\frac{{\mathrm{\Δa}}_i}{{\mathrm{\Δa}}_1}=\frac{1}{{\mathrm{\λ}}_i}\frac{{\mathrm{cos\θ}}_1}{{\mathrm{cos\θ}}_i}=\frac{1}{{\mathrm{\λ}}_i}\frac{\frac{a}{a+{\mathrm{\Δa}}_1}}{\frac{a}{a+{\mathrm{\Δa}}_i}}=\frac{1}{{\mathrm{\λ}}_i}\frac{1+\frac{{\mathrm{\Δa}}_i}{a}}{1+\frac{{\mathrm{\Δa}}_1}{a}}---\left(11\right)$

$\frac{\frac{{\mathrm{\Δa}}_i}{a}}{\frac{{\mathrm{\Δa}}_1}{a}}=\frac{1}{{\mathrm{\λ}}_i}\frac{1+\frac{{\mathrm{\Δa}}_i}{a}}{1+\frac{{\mathrm{\Δa}}_1}{a}}---\left(12\right)$

令是第i个切向弹簧的变形率，则有

${\mathrm{\ϵ}}_i=\frac{{\mathrm{\ϵ}}_1}{{\mathrm{\λ}}_i+({\mathrm{\λ}}_i-1){\mathrm{\ϵ}}_1}---\left(13\right)$

法线方向上第i和第i+1个点之间的相对形变：

$l_i=\sqrt{{(a+{\mathrm{\Δa}}_i)}^2-a^2}=\sqrt{(2{\mathrm{\ϵ}}_i+{\mathrm{\ϵ}}_i^2)}a---\left(14\right)$

将公式(13)代入(14)，可得：

$\frac{l_i}{a}=\frac{1}{{\mathrm{\λ}}_i+({\mathrm{\λ}}_i-1){\mathrm{\ϵ}}_1}\sqrt{\[2{\mathrm{\λ}}_i+(2{\mathrm{\λ}}_i-1){\mathrm{\ϵ}}_1\]{\mathrm{\ϵ}}_1}---\left(15\right)$

当i＝1时，有：

所以： $\frac{l_i}{l_1}=\frac{1}{(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_i-{\mathrm{\ϵ}}_1}\sqrt{\frac{2(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_i-{\mathrm{\ϵ}}_1}{2+{\mathrm{\ϵ}}_1}}---\left(16\right)$

令ε_i＜0.5，当i≥2时有λ_i＞3，则

(1+ε_i)λ_i＞＞ε_i

则： $\frac{l_i}{l_1}\≈\sqrt{\frac{2}{(1+{\mathrm{\ϵ}}_1)(2+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_i}}=\frac{1}{\mathrm{\ρ}}\·\frac{1}{\sqrt{{\mathrm{\λ}}_i}},i=2,3....,m---\left(17\right)$

其中所以有以下重要关系式：

$l_i=\frac{1}{\mathrm{\ρ}}\·\frac{1}{\sqrt{{\mathrm{\λ}}_i}}\·l_1,i=2,3....,m---\left(18\right)$

每个质点的总变形量为：

$\left\{\begin{array}{c}{\mathrm{\Δz}}_1=l_1+l_2+l_3+...+l_m\\ {\mathrm{\Δz}}_2=l_2+l_3+...+l_m\\ ...\\ {\mathrm{\Δz}}_m=l_m\end{array}\right.---\left(19\right)$

将式(19)、(9)代入式(6)有

$F_z=k_{N_1}\underset{i=1}{\overset{m}{\Σ}}(2i-1){\mathrm{\Δz}}_i---\left(20\right)$

$F_z=k_{N_1}\underset{i=1}{\overset{m}{\Σ}}{i}^2{l}_i---\left(21\right)$

其中，m表示在虚拟物体上以受力点为中心，划分的同心的等边三角形的个数，由方程(15)得到

$\frac{l_m}{a}=\frac{1}{(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_m-{\mathrm{\ϵ}}_1}\sqrt{\[2(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_m-{\mathrm{\ϵ}}_1\]{\mathrm{\ϵ}}_1}---\left(22\right)$

因为(1+ε₁)λ_m＞＞ε₁，所以

$\frac{l_m}{a}\≈\frac{1}{(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_m}\sqrt{\[2(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_m\]{\mathrm{\ϵ}}_1}=\sqrt{\frac{2{\mathrm{\ϵ}}_1}{(1+{\mathrm{\ϵ}}_1){\mathrm{\λ}}_m}}---\left(23\right)$

由上式得到： $m=\mathrm{\β}ln\left(\frac{2{\mathrm{\ϵ}}_1}{1+{\mathrm{\ϵ}}_1}{\left(\frac{a}{l_m}\right)}^2\right)+1---\left(24\right)$

其中，l_m/a代表了计算区域边界的变形量，而ε₁表示受力中心点的切向弹簧的变形率；

(3)计算虚拟物体发生变形的区域，获得发生变形的弹簧层数，计算发生变形的区域内各层弹簧变形量之和为虚拟物体的变形，各层弹簧弹力的合力为力反馈。

2.根据权利要求1所述的用于力触觉人机交互的柔性物体变形仿真方法，其特征在于：所述β取值为4.0，受力中心点的切向弹簧的变形率ε₁取值为0.3，l_m/a取值为1/57.3。