CN101635063A - 一种椭圆截面弹性柱体自由扭转/变形的模拟方法 - Google Patents

Abstract

本发明公开了一种椭圆截面弹性柱体自由扭转/变形的模拟方法，其特征是以椭圆截面弹性柱体的棱线为经线，以截面线为纬线，通过相邻的经线、相邻的纬线将椭圆截面弹性柱体侧面分割成一系列块状区域；在椭圆截面弹性柱体自由端施加扭矩作用下，使其产生自由扭转变形，利用扭转应力函数与单位长度的相对扭转角及扭矩之间的对应关系可求出任一关键点在扭矩作用下的变形量；采用插值计算计算任意块状区域内任意一点的变形量，进而更新绘制扭转后的弹性柱体模型；该模拟方法计算简单，能准确快速的计算扭转变形，实现对弹性柱体的实时变形仿真。在交互过程中，自然舒适、力触觉感觉平稳、模拟效果逼真。

Description

一种椭圆截面弹性柱体自由扭转/变形的模拟方法

技术领域

本发明涉及一种扭转/变形的模拟方法，尤其涉及一种用于虚拟现实人机交互的基于物理意义的椭圆截面弹性柱体自由扭转/变形的模拟方法。

背景技术

虚拟现实中的力/触觉再现是指用户通过力/触觉设备进行触摸、感知和操纵虚拟物体等一系列相互作用来获得虚拟物体特性，表征力/触觉感知信息的过程。用户能够以自然的方式与虚拟环境交互，从而产生与真实环境一致的沉浸感，在远程医疗、遥操作、虚拟制造、虚拟手术等领域有着广阔的应用前景。

虚拟弹性体的力/触觉变形模拟方法是力/触觉再现技术中最为重要的环节。目前，虚拟物体变形按照变形成因可分力/变形和力矩/变形。其中，虚拟力矩/变形研究在虚拟手术的力/触觉人机交互已成为愈来愈受关注的前沿性研究课题。在虚拟手术中，当医生进入虚拟世界，通过身体(尤其是手和手臂)的运动，与虚拟模型和环境进行交互，获取真实力/触觉和动感反馈，从而形成对虚拟器官的一个完整的认识。尤其对于肠道、血管等柱状体，切割、拉伸、搅动等操作型任务引起的扭转变形的实时性和准确性是手术成功的关键。因此，弹性柱体自由扭转/变形技术无疑对推到虚拟手术的发展具有重要意义。

发明内容

本发明提出一种椭圆截面弹性柱体自由扭转/变形的模拟方法，并将其用于虚拟现实人机交互的弹性柱体组织的变形仿真。该方法能准确快速的计算扭转变形量，实现对弹性组织的实时变形仿真，从而提高虚拟力触觉交互的逼真度。

本发明采用如下技术方案：

一种椭圆截面弹性柱体自由扭转/变形的模拟方法，其特征在于具体步骤如下：

步骤1将椭圆截面弹性柱体的一端固定，并以固定端面所在的平面为XY平面，以椭圆截面弹性柱体的中轴线为Z轴，建立XYZ空间坐标系，XYZ空间坐标系的原点为固定端面的形心O，以椭圆截面弹性柱体的另一端为自由端；在椭圆截面弹性柱体的侧面上选取出i条棱线且相邻棱线之间的距离相等，其中i的取值范围为2～360，取j个与XY平面平行的截面且相邻截面之间的距离相等，其中j的取值范围为2～360，且截面中包含固定端面和自由端的端面；以各个截面与各条棱线的交点作为关键点，以截面与棱线的交线为截面线，为棱线为经线，以截面线为纬线，相邻的经线、相邻的纬线将弹性柱体侧面分割成一系列块状区域，块状区域个数为：

s＝i·(j-1)   i＝2，3，......，360   j＝2，3，......，360(1)

其中s为椭圆截面弹性柱体侧面被分割成的块状区域个数，i为侧面上选取的棱线条数，j为选取的截面个数；

步骤2假定任一块状区域K由A_k、B_k、C_k、D_k四个关键点围成，且对应各个关键点在扭矩作用下的位移变形量分别为

$\mathrm{\Δ}{B}_k(\mathrm{\Δ}{x}_{B_k},\mathrm{\Δ}{y}_{B_k},\mathrm{\Δ}{z}_{B_k}),\mathrm{\Δ}{C}_k(\mathrm{\Δ}{x}_{C_k},\mathrm{\Δ}{y}_{C_k},\mathrm{\Δ}{z}_{C_k}),\mathrm{\Δ}{D}_k(\mathrm{\Δ}{x}_{D_k},\mathrm{\Δ}{y}_{D_k},{\mathrm{\Δz}}_{D_k});$ 块状区域内存在任一点Q_k，且Q_k离A_k、B_k、C_k、D_k四个关键点的直线距离分别为

直线距离可用空间任意两点的距离公式求出；这四个关键点对Q_k位移变形量影响系数分别记为且其影响系数之和为1；

$\left[\begin{array}{c}{\mathrm{\ρ}}_{A_k}^{Q_k}\\ {\mathrm{\ρ}}_{B_k}^{Q_k}\\ {\mathrm{\ρ}}_{C_k}^{Q_k}\\ {\mathrm{\ρ}}_{D_k}^{Q_k}\end{array}\right]=\left[\begin{array}{c}{d}_{B_k}^{Q_k}{d}_{C_k}^{Q_k}{d}_{D_k}^{Q_k}\\ d_{A_k}^{Q_k}{d}_{C_k}^{Q_k}{d}_{D_k}^{Q_k}\\ d_{A_k}^{Q_k}{d}_{B_k}^{Q_k}{d}_{D_k}^{Q_k}\\ d_{A_k}^{Q_k}{d}_{B_k}^{Q_k}{d}_{C_k}^{Q_k}\end{array}\right]/(d_{B_k}^{Q_k}{d}_{C_k}^{Q_k}{d}_{D_k}^{Q_k}+d_{A_k}^{Q_k}{d}_{C_k}^{Q_k}{d}_{D_k}^{Q_k}+d_{A_k}^{Q_k}{d}_{B_k}^{Q_k}{d}_{D_k}^{Q_k}+d_{A_k}^{Q_k}{d}_{B_k}^{Q_k}{d}_{C_k}^{Q_k})---\left(2\right)$

则Q_k点的位移变形量

可由A_k、B_k、C_k、D_k四个关键点的位移变形量及影响系数的线性插值求出：

$\mathrm{\Δ}{Q}_k={\mathrm{\ρ}}_{A_k}^{Q_k}\·\mathrm{\Δ}{A}_k+{\mathrm{\ρ}}_{B_k}^{Q_k}\·\mathrm{\Δ}{B}_k+{\mathrm{\ρ}}_{C_k}^{Q_k}\·\mathrm{\Δ}{C}_k+{\mathrm{\ρ}}_{D_k}^{Q_k}\·\mathrm{\Δ}{D}_k---\left(3\right)$

上述任一关键点的位移变形量的计算方法相同，可由步骤3求得；步骤3在椭圆截面弹性柱体的自由端施加一扭矩M_T，使椭圆截面弹性柱体产生自由扭转变形，在选取的任意水平截面XY方向上存在一个扭转应力函数

所述扭转应力函数与沿XY方向的变量x、y满足：

其中，a、b为椭圆的长半轴、短半轴，G为剪切模量，其与弹性材料的性质有关，θ为单位长度的相对扭转角，这里单位长度指自由端面与固定端面之间设定的单位垂直距离；

扭转应力函数

与扭转力矩M_T存在以下函数关系：

其中

为扭转应力函数

在截面 $\frac{x^2}{a^2}+\frac{y^2}{b^2}\≤1$ 上沿x、y方向的二重积分；

椭圆截面弹性柱体扭转后，任意截面的扭转角与该截面到固定端面的垂直距离成正比，则离固定端面垂直距离为h处截面的扭转角为h·θ；假定初始位置时，在离截面形心O′为距离r的任一点P_o(x_o，y_o，z_o)，其中r的取值范围介于a、b之间。在扭矩M_T的作用下，扭转后到达一个新的位置P_e(x_e，y_e，z_e)，则该点沿x、y方向的位移分量变形量Δx、Δy可分别表示为：

Δx＝x_e-x_o＝-r·sina·h·θ＝-y·h·θ(6)

Δy＝y_e-y_o＝r·cosa·h·θ＝x·h·θ(7)

椭圆截面弹性柱体中任意截面的翘曲变形与其离固定端面的垂直距离h无关，这里的翘曲指在扭矩作用下，截面上各点沿Z轴方向产生的位移分量，取任意横截面，对于椭圆截面弹性柱体其对应的翘曲函数Δz是该截平面上自变量x，y坐标的函数，即：

$\mathrm{\Δz}=f(x,y)=-\frac{a^2-b^2}{{\mathrm{\πa}}^3{b}^{3}G}{M}_T\mathrm{xy}---\left(8\right)$

步骤4根据步骤2、3，获得各个关键点及各个块状区域内任意一点在三维空间里的坐标，更新绘制扭转后的弹性柱体模型。

本发明的优点：

(1)和以往常用的基于物理意义的扭转/变形仿真力触觉模拟方法相比，该模拟方法以弹性力学为理论基础，将椭圆截面弹性柱体侧面分割成一系列块状区域，任意一块状区域的任意一个关键点的位移变形量计算方法相同，从而简化计算，加快了变形计算的速度。

(2)对于任意一块状区域内的任一点的位移变形量，通过组成该块状区域的四个关键点的位移变形量及影响系数的线性插值求出，线性插值计算简单方便，使得此模拟方法具有快速灵活的特点。

(3)通过调节划分块状区域的个数、弹性材料的剪切模量、单位长度的相对扭转角等参数，就可模拟不同类型的弹性柱体模型，具有广泛的适用性。

(4)可将其应用于虚拟外科手术仿真、遥操作机器人控制、远程医疗等领域。

附图说明

图1椭圆截面弹性柱体任意关键点的分布及块状区域内任意点插值变形计算的思想；

图2椭圆截面弹性柱体侧面上任意关键点的扭转变形；

图3椭圆截面弹性柱体力矩/扭转模拟方法流程图；

图4该模拟方法的系统图；

图5椭圆截面弹性柱体扭转变形实例图，(a)是椭圆截面弹性柱体未受扭矩示意图，(b)是椭圆截面弹性柱体受扭矩作用下的示意图。

具体实施方式：

实施例1：

一种椭圆截面弹性柱体自由扭转/变形的模拟方法，其特征在于具体步骤如下：

步骤1将椭圆截面弹性柱体的一端固定，并以固定端面所在的平面为XY平面，以椭圆截面弹性柱体的中轴线为Z轴，建立XYZ空间坐标系，XYZ空间坐标系的原点为固定端面的形心O，以椭圆截面弹性柱体的另一端为自由端；在椭圆截面弹性柱体的侧面上选取出i条棱线且相邻棱线之间的距离相等，其中i的取值范围为2～360，取j个与XY平面平行的截面且相邻截面之间的距离相等，其中j的取值范围为2～360，且截面中包含固定端面和自由端的端面；以各个截面与各条棱线的交点作为关键点，以截面与棱线的交线为截面线，为棱线为经线，以截面线为纬线，相邻的经线、相邻的纬线将弹性柱体侧面分割成一系列块状区域，块状区域个数为：

s＝i·(j-1)  i＝2，3，......，360  j＝2，3，......，360(1)

其中s为椭圆截面弹性柱体侧面被分割成的块状区域个数，i为侧面上选取的棱线条数，j为选取的截面个数；

步骤2假定任一块状区域K由A_k、B_k、C_k、D_k四个关键点围成，且对应各个关键点在扭矩作用下的位移变形量分别为

$\mathrm{\Δ}{B}_k(\mathrm{\Δ}{x}_{B_k},\mathrm{\Δ}{y}_{B_k},\mathrm{\Δ}{z}_{B_k}),\mathrm{\Δ}{C}_k(\mathrm{\Δ}{x}_{C_k},\mathrm{\Δ}{y}_{C_k},\mathrm{\Δ}{z}_{C_k}),\mathrm{\Δ}{D}_k(\mathrm{\Δ}{x}_{D_k},\mathrm{\Δ}{y}_{D_k},{\mathrm{\Δz}}_{D_k});$ 块状区域内存在任一点Q_k，且Q_k离A_k、B_k、C_k、D_k四个关键点的直线距离分别为

直线距离可用空间任意两点的距离公式求出；这四个关键点对Q_k位移变形量影响系数分别记为

且其影响系数之和为1；

$\left[\begin{array}{c}{\mathrm{\ρ}}_{A_k}^{Q_k}\\ {\mathrm{\ρ}}_{B_k}^{Q_k}\\ {\mathrm{\ρ}}_{C_k}^{Q_k}\\ {\mathrm{\ρ}}_{D_k}^{Q_k}\end{array}\right]=\left[\begin{array}{c}{d}_{B_k}^{Q_k}{d}_{C_k}^{Q_k}{d}_{D_k}^{Q_k}\\ d_{A_k}^{Q_k}{d}_{C_k}^{Q_k}{d}_{D_k}^{Q_k}\\ d_{A_k}^{Q_k}{d}_{B_k}^{Q_k}{d}_{D_k}^{Q_k}\\ d_{A_k}^{Q_k}{d}_{B_k}^{Q_k}{d}_{C_k}^{Q_k}\end{array}\right]/(d_{B_k}^{Q_k}{d}_{C_k}^{Q_k}{d}_{D_k}^{Q_k}+d_{A_k}^{Q_k}{d}_{C_k}^{Q_k}{d}_{D_k}^{Q_k}+d_{A_k}^{Q_k}{d}_{B_k}^{Q_k}{d}_{D_k}^{Q_k}+d_{A_k}^{Q_k}{d}_{B_k}^{Q_k}{d}_{C_k}^{Q_k})---\left(2\right)$

则Q_k点的位移变形量

可由A_k、B_k、C_k、D_k四个关键点的位移变形量及影响系数的线性插值求出：

$\mathrm{\Δ}{Q}_k={\mathrm{\ρ}}_{A_k}^{Q_k}\·\mathrm{\Δ}{A}_k+{\mathrm{\ρ}}_{B_k}^{Q_k}\·\mathrm{\Δ}{B}_k+{\mathrm{\ρ}}_{C_k}^{Q_k}\·\mathrm{\Δ}{C}_k+{\mathrm{\ρ}}_{D_k}^{Q_k}\·\mathrm{\Δ}{D}_k---\left(3\right)$

上述任一关键点的位移变形量的计算方法相同，可由步骤3求得；

步骤3在椭圆截面弹性柱体的自由端施加一扭矩M_T，使椭圆截面弹性柱体产生自由扭转变形，在选取的任意水平截面XY方向上存在一个扭转应力函数

所述扭转应力函数与沿XY方向的变量x、y满足：

其中，a、b为椭圆的长半轴、短半轴，G为剪切模量，其与弹性材料的性质有关，θ为单位长度的相对扭转角，这里单位长度指自由端面与固定端面之间设定的单位垂直距离；

扭转应力函数

与扭转力矩M_T存在以下函数关系：

其中

为扭转应力函数

在截面 $\frac{x^2}{a^2}+\frac{y^2}{b^2}\≤1$ 上沿x、y方向的二重积分；

椭圆截面弹性柱体扭转后，任意截面的扭转角与该截面到固定端面的垂直距离成正比，则离固定端面垂直距离为h处截面的扭转角为h·θ；假定初始位置时，在离截面形心O′为距离r的任一点P_o(x_o，y_o，z_o)，其中r的取值范围介于a、b之间。在扭矩M_T的作用下，扭转后到达一个新的位置P_e(x_e，y_e，z_e)，则该点沿x、y方向的位移分量变形量Δx、Δy可分别表示为：

Δx＝x_e-x_o＝-r·sin a·h·θ＝-y·h·θ(6)

Δy＝y_e-y_o＝r·cosa·h·θ＝x·h·θ(7)

椭圆截面弹性柱体中任意截面的翘曲变形与其离固定端面的垂直距离h无关，这里的翘曲指在扭矩作用下，截面上各点沿Z轴方向产生的位移分量，取任意横截面，对于椭圆截面弹性柱体其对应的翘曲函数Δz是该截平面上自变量x，y坐标的函数，即：

$\mathrm{\Δz}=f(x,y)=-\frac{a^2-b^2}{{\mathrm{\πa}}^3{b}^{3}G}{M}_T\mathrm{xy}---\left(8\right)$

步骤4根据步骤2、3，获得各个关键点及各个块状区域内任意一点在三维空间里的坐标，更新绘制扭转后的弹性柱体模型。

具体实施例2：

1、对椭圆截面弹性柱体进行量化，并对其侧面进行块状区域划分。

本实例中椭圆截面弹性柱体长半轴a为0.04m、短半轴b为0.03m、高为0.2m，椭圆截面弹性柱体的一端固定，并以长半轴和短半轴所在的平面为XY平面，以其高所在方向中轴线为Z轴，建立XYZ空间坐标系，XYZ空间坐标系的原点为固定端面的形心O，以椭圆截面弹性柱体的另一端为自由端；在椭圆截面弹性柱体的侧面上选取出6条棱线且相邻棱线之间的距离相等，取5个与XY平面平行的截面且相邻截面之间的距离相等，且截面中包含固定端面和自由端的端面；以各个截面与各条棱线的交点作为关键点，以截面与棱线的交线为截面线，为棱线为经线，以截面线为纬线，相邻的经线、相邻的纬线将弹性柱体侧面分割成一系列块状区域，块状区域个数为：

s＝i·(j-1)＝6·(5-1)＝24(1)

其中s为椭圆截面弹性柱体侧面被分割成的块状区域个数，i为侧面上选取的棱线条数，j为选取的截面个数；

2、选取椭圆侧面上任一块状区域K内的A_k、B_k、C_k、D_k四个关键点及块状区域K内的任意一点Q_k，并计算Q_k离A_k、B_k、C_k、D_k四个关键点的直线距离分别为以下计算的数据按四舍五入法保留三位有效数字；例如，关键点及区域内一点得坐标为：A_k(0.0159，-0.0275，0.15)、B_k(0.0159，-0.0275，0.1)、C_k(-0.0159，-0.0275，0.15)、D_k(-0.0159，-0.0275，0.1)、Q_k(0，-0.03，0.12)

Q_k到四个关键点A_k、B_k、C_k、D_k的距离，用空间任意两点的直线公式分别求得， $d_{A_k}^{Q_k}=0.0162m,$ $d_{B_k}^{Q_k}=0.0160m,$ $d_{C_k}^{Q_k}=0.0162m,$ $d_{D_k}^{Q_k}=0.0160m,$ 这四个关键点对Q_k位移变形量影响系数分别记为

且其影响系数之和为1；

$\left[\begin{array}{c}{\mathrm{\ρ}}_{A_k}^{Q_k}\\ {\mathrm{\ρ}}_{B_k}^{Q_k}\\ {\mathrm{\ρ}}_{C_k}^{Q_k}\\ {\mathrm{\ρ}}_{D_k}^{Q_k}\end{array}\right]=\left[\begin{array}{c}{d}_{B_k}^{Q_k}{d}_{C_k}^{Q_k}{d}_{D_k}^{Q_k}\\ d_{A_k}^{Q_k}{d}_{C_k}^{Q_k}{d}_{D_k}^{Q_k}\\ d_{A_k}^{Q_k}{d}_{B_k}^{Q_k}{d}_{D_k}^{Q_k}\\ d_{A_k}^{Q_k}{d}_{B_k}^{Q_k}{d}_{C_k}^{Q_k}\end{array}\right]/(d_{B_k}^{Q_k}{d}_{C_k}^{Q_k}{d}_{D_k}^{Q_k}+d_{A_k}^{Q_k}{d}_{C_k}^{Q_k}{d}_{D_k}^{Q_k}+d_{A_k}^{Q_k}{d}_{B_k}^{Q_k}{d}_{D_k}^{Q_k}+d_{A_k}^{Q_k}{d}_{B_k}^{Q_k}{d}_{C_k}^{Q_k})$

$=\left[\begin{array}{c}0.0160\×0.0162\×0.0160\\ 0.0162\×0.0162\×0.0160\\ 0.0162\×0.0160\×0.0160\\ 0.0162\×0.0160\×0.0162\end{array}\right]/(0.0160\×0.0162\×0.0160+0.0162\×0.0162\×0.0160$

$+0.0162\×0.0160\×0.0160+0.0162\×0.0160\×0.0162)$

$=\left[\begin{array}{c}0.248\\ 0.252\\ 0.248\\ 0.252\end{array}\right]---\left(2\right)$

3、在椭圆弹性柱体的自由端施加一逆时针扭矩M_T，值为5N·m，使弹性柱体产生自由扭转变形。

在选取的任意水平截面XY方向上存在一个扭转应力函数

对于椭圆截面所述扭转应力函数与沿XY方向的变量x、y满足：

其中，a、b为椭圆的长半轴、短半轴，剪切模量G为2.93MPa，θ为单位长度的相对扭转角，这里单位长度指自由端面与固定端面之间设定的单位垂直距离，这里设定为1米；

扭转应力函数

与扭转力矩M_T存在以下函数关系：

即 $\mathrm{\θ}=\frac{a^2+b^2}{\mathrm{\π}{a}^3{b}^{3}G}{M}_T---\left(5\right)$

其中为扭转应力函数

在截面 $\frac{x^2}{a^2}+\frac{y^2}{b^2}\≤1$ 上沿x、y方向的二重积分；

椭圆弹性柱体扭转后，任意截面的扭转角与该截面到固定端面的垂直距离成正比，则离固定端面垂直距离为h处截面的扭转角为h·θ；分别计算四个关键点A_k、B_k、C_k、D_k沿x，y方向的位移变形量分别为；

A_k在x方向的位移变形量为：

$\mathrm{\Δ}{x}_{A_k}=-y_{A_k}\·\mathrm{\θ}\·h_{A_k}=-\frac{(a^2+b^2)}{\mathrm{\π}{a}^3{b}^{3}G}{M}_T{y}_{A_k}{h}_{A_k}$

$=-\frac{[{\left(0.04\right)}^2+{\left(0.03\right)}^2]}{3.14\×{\left(0.04\right)}^3\×{\left(0.03\right)}^3\×2.93\×{10}^6}\×5\×(-0.0275)\×0.15$

$=0.00324$

A_k在y方向的位移变形量为：

$\mathrm{\Δ}{y}_{A_k}=x_{A_k}\mathrm{\θ}{h}_{A_k}=\frac{(a^2+b^2)}{\mathrm{\π}{a}^3{b}^{3}G}{M}_T{x}_{A_k}{h}_{A_k}$

$=\frac{[{\left(0.04\right)}^2+{\left(0.03\right)}^2]}{3.14\×{\left(0.04\right)}^3\×{\left(0.03\right)}^3\×2.93\×{10}^6}\×5\×0.0159\×0.15$

$=0.00188$

同理，对于区域K内其它的三个关键点B_k(0.0159，-0.0275，0.1)、C_k(-0.0159，-0.0275，0.015)、D_k(-0.0159，-0.0275，0.1)沿x，y方向的位移分量分别为：

$\mathrm{\Δ}{x}_{B_k}=-y_{B_k}\mathrm{\θ}{h}_{B_k}=-\frac{(a^2+b^2)}{\mathrm{\π}{a}^3{b}^{3}G}{M}_T{y}_{B_k}{h}_{B_k}$

$=-\frac{[{\left(0.04\right)}^2+{\left(0.03\right)}^2]}{3.14\×{\left(0.04\right)}^3\×{\left(0.03\right)}^3\×2.93\×{10}^6}\×5\×(-0.0275)\×0.1$

$=0.00216$

$\mathrm{\Δ}{y}_{B_k}=x_{B_k}\mathrm{\θ}{h}_{B_k}=\frac{(a^2+b^2)}{\mathrm{\π}{a}^3{b}^{3}G}{M}_T{x}_{B_k}{h}_{B_k}$

$=\frac{[{\left(0.04\right)}^2+{\left(0.03\right)}^2]}{3.14\×{\left(0.04\right)}^3\×{\left(0.03\right)}^3\×2.93\×{10}^6}\×5\×0.0159\×0.1$

$=0.00125$

$\mathrm{\Δ}{x}_{C_k}=-y_{C_k}\mathrm{\θ}{h}_{C_k}=-\frac{(a^2+b^2)}{\mathrm{\π}{a}^3{b}^{3}G}{M}_T{y}_{C_k}{h}_{C_k}$

$=-\frac{[{\left(0.04\right)}^2+{\left(0.03\right)}^2]}{3.14\×{\left(0.04\right)}^3\×{\left(0.03\right)}^3\×2.93\×{10}^6}\×5\×(-0.0275)\×0.15$

$=0.00324$

$\mathrm{\Δ}{y}_{C_k}=x_{C_k}\mathrm{\θ}{h}_{C_k}=\frac{(a^2+b^2)}{\mathrm{\π}{a}^3{b}^{3}G}{M}_T{x}_{C_k}{h}_{C_k}$

$\begin{array}{c}=\frac{[{\left(0.04\right)}^2+{\left(0.03\right)}^2]}{3.14\×{\left(0.04\right)}^3\×{\left(0.03\right)}^3\×2.93\×{10}^6}\×5\×(-0.0159)\×0.15\\ =-0.00188\end{array}$

$\mathrm{\Δ}{x}_{D_k}=-y_{D_k}\mathrm{\θ}{h}_{D_k}=-\frac{(a^2+b^2)}{\mathrm{\π}{a}^3{b}^{3}G}{M}_T{y}_{D_k}{h}_{D_k}$

$=-\frac{[{\left(0.04\right)}^2+{\left(0.03\right)}^2]}{3.14\×{\left(0.04\right)}^3\×{\left(0.03\right)}^3\×2.93\×{10}^6}\×5\×(-0.0275)\×0.1$

$=0.00216$

$\mathrm{\Δ}{y}_{D_k}=x_{\mathrm{Dk}}\mathrm{\θ}{h}_{D_k}=\frac{(a^2+b^2)}{\mathrm{\π}{a}^3{b}^{3}G}{M}_T{x}_{D_k}{h}_{D_k}$

$=\frac{[{\left(0.04\right)}^2+{\left(0.03\right)}^2]}{3.14\×{\left(0.04\right)}^3\×{\left(0.03\right)}^3\×2.93\×{10}^6}\×5\×\left(0.0159\right)\×0.1$

$=-0.00125$

椭圆弹性柱体中任意截面的翘曲变形与其离固定端面的垂直距离h无关，这里的翘曲指在扭矩作用下，截面上各点沿Z轴方向产生的位移分量，取任意横截面，对于区域K上的四个关键点的翘曲函数Δz是该截平面上自变量x，y坐标的函数，即：

$\mathrm{\Δ}{z}_{A_k}=f(x_{A_k},y_{A_k})=-\frac{a^2-b^2}{a^2+b^2}\mathrm{\θ}{x}_{A_k}{y}_{A_k}=-\frac{(a^2-b^2)}{\mathrm{\π}{a}^3{b}^{3}G}{M}_T{x}_{A_k}{y}_{A_k}$

$=-\frac{[{\left(0.04\right)}^2-{\left(0.03\right)}^2]}{3.14\×{\left(0.04\right)}^3\×{\left(0.03\right)}^3\×2.93\×{10}^6}\×5\×0.0159\×(-0.0275)$

$=0.000963$

$\mathrm{\Δ}{z}_{B_k}=f(x_{B_k},y_{B_k})=-\frac{a^2-b^2}{a^2+b^2}\mathrm{\θ}{x}_{B_k}{y}_{B_k}=-\frac{(a^2-b^2)}{\mathrm{\π}{a}^3{b}^{3}G}{M}_T{x}_{B_k}{y}_{B_k}$

$=-\frac{[{\left(0.04\right)}^2-{\left(0.03\right)}^2]}{3.14\×{\left(0.04\right)}^3\×{\left(0.03\right)}^3\×2.93\×{10}^6}\×5\×0.0159\×(-0.0275)$

$=0.000963$

$\mathrm{\Δ}{z}_{C_k}=f(x_{C_k},y_{C_k})=-\frac{a^2-b^2}{a^2+b^2}\mathrm{\θ}{x}_{C_k}{y}_{C_k}=-\frac{(a^2-b^2)}{\mathrm{\π}{a}^3{b}^{3}G}{M}_T{x}_{C_k}{y}_{C_k}$

$=-\frac{[{\left(0.04\right)}^2-{\left(0.03\right)}^2]}{3.14\×{\left(0.04\right)}^3\×{\left(0.03\right)}^3\×2.93\×{10}^6}\×5\×(-0.0159)\×(-0.0275)$

$=-0.000963$

$\mathrm{\Δ}{z}_{D_k}=f(x_{D_k},y_{D_k})=-\frac{a^2-b^2}{a^2+b^2}\mathrm{\θ}{x}_{D_k}{y}_{D_k}=-\frac{(a^2-b^2)}{\mathrm{\π}{a}^3{b}^{3}G}{M}_T{x}_{D_k}{y}_{D_k}$

$=-\frac{[{\left(0.04\right)}^2-{\left(0.03\right)}^2]}{3.14\×{\left(0.04\right)}^3\×{\left(0.03\right)}^3\×2.93\×{10}^6}\×5\×(-0.0159)\×(-0.0275)$

$=-0.000963$

因此，在k区域内四个关键点的位移变形量分别为：ΔA_k(0.00324，0.00188，0.000963)、ΔB_k(0.00216，0.00125，0.000963)、ΔC_k(0.00324，-0.00188，-0.000963)、ΔD_k(0.00216，-0.00125，-0.000963)。

则点Q_k在三坐标轴方向上的变形量

可由块状区域上的A_k、B_k、C_k、D_k关键点的位移线性插值求解。

$\left[\begin{array}{c}\mathrm{\Δ}{x}_{Q_k}\\ \mathrm{\Δ}{y}_{Q_k}\\ \mathrm{\Δ}{z}_{Q_k}\end{array}\right]={\mathrm{\ρ}}_{A_k}^{Q_k}\left[\begin{array}{c}\mathrm{\Δ}{x}_{A_k}\\ \mathrm{\Δ}{y}_{A_k}\\ \mathrm{\Δ}{z}_{A_k}\end{array}\right]+{\mathrm{\ρ}}_{B_k}^{Q_k}\left[\begin{array}{c}\mathrm{\Δ}{x}_{B_k}\\ \mathrm{\Δ}{y}_{B_k}\\ \mathrm{\Δ}{z}_{B_k}\end{array}\right]+{\mathrm{\ρ}}_{C_k}^{Q_k}\left[\begin{array}{c}\mathrm{\Δ}{x}_{C_k}\\ \mathrm{\Δ}{y}_{C_k}\\ \mathrm{\Δ}{z}_{C_k}\end{array}\right]+{\mathrm{\ρ}}_{D_k}^{Q_k}\left[\begin{array}{c}\mathrm{\Δ}{x}_{D_k}\\ \mathrm{\Δ}{y}_{D_k}\\ \mathrm{\Δ}{z}_{D_k}\end{array}\right]$

$=0.248\×\left[\begin{array}{c}0.00324\\ 0.00188\\ 0.000963\end{array}\right]+0.252\×\left[\begin{array}{c}0.00216\\ 0.00125\\ 0.000963\end{array}\right]+0.248\×\left[\begin{array}{c}0.00324\\ -0.00188\\ -0.000963\end{array}\right]+0.252\×\left[\begin{array}{c}0.00216\\ -0.00125\\ -0.000963\end{array}\right]$

$=\left[\begin{array}{c}0.00270\\ 0\\ 0\end{array}\right]$

4、根据上面方法，获得各个关键点及各个块状区域内任意一点在三维空间里的坐标，更新绘制扭转后的弹性柱体模型。

为验证本发明的实施效果，操作者通过6-DOF Delta手控器(工作空间范围为平移Φ360mm×L300mm和旋转±20deg/axis，能够实现沿Z轴的作用力矩触觉反馈)使得椭圆截面柱体产生扭转的变形，并将交互过程中产生的力触觉信息实时反馈给操作者。实验结果表明：该模拟方法不仅计算简单，而且能够保证触觉接触力和变形计算具有较高精度，模拟效果逼真。

Claims (1)

1、一种椭圆截面弹性柱体自由扭转/变形的模拟方法，其特征在于具体步骤如下：

步骤1将椭圆截面弹性柱体的一端固定，并以固定端面所在的平面为XY平面，以椭圆截面弹性柱体的中轴线为Z轴，建立XYZ空间坐标系，XYZ空间坐标系的原点为固定端面的形心O，以椭圆截面弹性柱体的另一端为自由端；在椭圆截面弹性柱体的侧面上选取出i条棱线且相邻棱线之间的距离相等，其中i的取值范围为2～360，取j个与XY平面平行的截面且相邻截面之间的距离相等，其中j的取值范围为2～360，且截面中包含固定端面和自由端的端面；以各个截面与各条棱线的交点作为关键点，以截面与棱线的交线为截面线，为棱线为经线，以截面线为纬线，相邻的经线、相邻的纬线将弹性柱体侧面分割成一系列块状区域，块状区域个数为：

s＝i·(j-1) i＝2，3，......，360 j＝2，3，......，360    (1)

其中s为椭圆截面弹性柱体侧面被分割成的块状区域个数，i为侧面上选取的棱线条数，j为选取的截面个数；

步骤2假定任一块状区域K由A_k、B_k、C_k、D_k四个关键点围成，且对应各个关键点在扭矩作用下的位移变形量分别为

块状区域内存在任一点Q_k，且Q_k离A_k、B_k、C_k、D_k四个关键点的直线距离分别为

这四个关键点对Q_k位移变形量影响系数分别记为

且其影响系数之和为1；

$\left[\begin{array}{c}{\mathrm{\ρ}}_{A_k}^{Q_k}\\ {\mathrm{\ρ}}_{B_k}^{Q_k}\\ {\mathrm{\ρ}}_{C_k}^{Q_k}\\ {\mathrm{\ρ}}_{D_k}^{Q_k}\end{array}\right]=\left[\begin{array}{c}{d}_{B_k}^{Q_k}{d}_{C_k}^{Q_k}{d}_{D_k}^{Q_k}\\ d_{A_k}^{Q_k}{d}_{C_k}^{Q_k}{d}_{D_k}^{Q_k}\\ d_{A_k}^{Q_k}{d}_{B_k}^{Q_k}{d}_{D_k}^{Q_k}\\ d_{A_k}^{Q_k}{d}_{B_k}^{Q_k}{d}_{C_k}^{Q_k}\end{array}\right]/\left(d_{B_k}^{Q_k}{d}_{C_k}^{Q_k}{d}_{D_k}^{Q_k}+d_{A_k}^{Q_k}{d}_{C_k}^{Q_k}{d}_{D_k}^{Q_k}+d_{A_k}^{Q_k}{d}_{B_k}^{Q_k}{d}_{D_k}^{Q_k}+d_{A_k}^{Q_k}{d}_{B_k}^{Q_k}{d}_{C_k}^{Q_k}\right)---\left(2\right)$

则Q_k点的位移变形量

可由A_k、B_k、C_k、D_k四个关键点的位移变形量及影响系数的线性插值求出：

${\mathrm{\ΔQ}}_k={\mathrm{\ρ}}_{A_k}^{Q_k}\·{\mathrm{\ΔA}}_k+{\mathrm{\ρ}}_{B_k}^{Q_k}\·{\mathrm{\ΔB}}_k+{\mathrm{\ρ}}_{C_k}^{Q_k}\·{\mathrm{\ΔC}}_k+{\mathrm{\ρ}}_{D_k}^{Q_k}\·{\mathrm{\ΔD}}_k---\left(3\right)$

上述任一关键点的位移变形量的计算方法相同，可由步骤3求得；

步骤3在椭圆截面弹性柱体的自由端施加一扭矩M_T，使椭圆截面弹性柱体产生自由扭转变形，在选取的任意水平截面XY方向上存在一个扭转应力函数所述扭转应力函数与沿XY方向的变量x、y满足：

其中，a、b为椭圆的长半轴、短半轴，G为剪切模量，其与弹性材料的性质有关，θ为单位长度的相对扭转角，这里单位长度指自由端面与固定端面之间设定的单位垂直距离；

扭转应力函数

与扭转力矩M_T存在以下函数关系：

其中

为扭转应力函数

在截面 $\frac{x^2}{a^2}+\frac{y^2}{b^2}\≤1$ 上沿XY方向的二重积分；

椭圆截面弹性柱体扭转后，任意截面的扭转角与该截面到固定端面的垂直距离成正比，则离固定端面垂直距离为h处截面的扭转角为h·θ；假定初始位置时，在离截面形心O′距离为r的任一点P_o(x_o，y_o，z_o)，其中r的取值范围介于a、b之间。在扭矩M_T的作用下，扭转后到达一个新的位置P_e(x_e，y_e，z_e)，则该点沿x、y方向的位移分量变形量Δx、Δy可分别表示为：

Δx＝x_e-x_o＝-r·sina·h·θ＝-y·h·θ   (6)

Δy＝y_e-y_o＝r·cosa·h·θ＝x·h·θ   (7)

椭圆截面弹性柱体中任意截面的翘曲变形与其离固定端面的垂直距离h无关，这里的翘曲指在扭矩作用下，截面上各点沿Z轴方向产生的位移分量，取任意截面，对于椭圆截面弹性柱体其对应的翘曲函数Δz是该截面上自变量x，y坐标的函数，即：

$\mathrm{\Δz}=f(x,y)=-\frac{a^2-b^2}{\mathrm{\π}{a}^3{b}^{3}G}{M}_T\mathrm{xy}---\left(8\right)$

步骤4根据步骤2、3，获得各个关键点及各个块状区域内任意一点在三维空间里的坐标，更新绘制扭转后的弹性柱体模型。

Images