基于Kirchhoff动力学比拟理论的柔性管件维修可视化

焦玉民 王 强 徐 婷 谢庆华

1.94679部队，南京，210038 2.解放军理工大学，南京，210007

0 引言

军用工程机械防护性能最薄弱的部位是暴露在机械外部的液压管路，液压管路一旦损坏，极易造成机械无法工作。虚拟维修训练是军用工程机械维修保障训练的重要手段之一，柔性零件（如软管、线缆、皮带等）的维修模拟是一个不容回避的问题。目前的虚拟维修仿真过程中，绝大多数维修对象均为刚性体（如壳体、齿轮、轴承），仿真过程仅考虑其在虚拟环境中的位置变化，不考虑其受力后的变形过程。而变形体的虚拟维修过程涉及单元类型、材料特性和外力载荷等问题，建模技术复杂。目前的虚拟现实环境中主要通过B样条描述柔性体的静态形状，采用颜色变换、声音提示以及动画等方式模拟其故障状态，难以模拟柔性体维修变形过程［1］。国内外学者在此方面开展了一系列研究。文献［2］采用了离散点方法来构建曲线模型，该模型在几何层面可以较好地模拟出曲线形态，但由于缺少曲线的物理特征描述，难以对不同类型的柔性体进行区分。文献［3］提出了基于弹簧－质点法的线缆弹性模型，使用弹簧的弹性系数描述线缆的物理特性，为柔性软管的运动求解提供了一种有效的解决办法，但弹簧－质点法依然存在局限性，该方法仅能计算软管的弯曲变形过程，无法对软管的扭转变形进行表达。刘检华等［4］、万毕乐等［5］对虚拟现实环境下线缆的建模与装配进行了研究，提出了离散控制点建模方法，对线缆的几何约束关系和装配方法进行了详细的分析，但并没有提及在外力影响下线缆的变形规律。刘延柱［6］和薛纭等［7］以DNA力学模型为参考，提出了基于Kirchhoff动力学比拟理论的弹性杆动力学模型，为虚拟现实环境下柔性体建模与分析提供了一种有效的解决方法。文献［8－9］对Kirchhoff动力学比拟理论进行了深入研究，并应用于不同的研究领域。

在上述研究的基础上，本文提出了虚拟现实环境下柔性管件空间形态的基本假设，并运用Kirchhoff动力学比拟理论对柔性体受力后的平衡状态进行建模分析，构建了柔性体动力学平衡方程，使用欧拉四元参数法对平衡方程进行求解，解决了虚拟现实环境下柔性体受力后的动态变形可视化问题，为在虚拟维修环境中实现柔性管件的布局优化以及维修过程的动态模拟等问题提供了切实有效的解决方案。

1 虚拟现实环境下柔性管件的基本假设

1.1 柔性管件描述

柔性体在外力作用下经受拉伸变形和扭转变形，其连接状态、材料性质和外力载荷都影响着变形过程，本文将软管受力后的变形简化为软管中心线的几何形状改变，为便于对柔性管件变形状态进行模拟计算，需对管件模型作出如下假设。几何特性方面：①将软管均分为n段，以每段软管端部横截面作为研究对象；②端部横截面均为刚性平面；③相邻截面可绕中心线做相对转动。物理特性方面：①忽略软管中心线的拉伸变形，任意两截面沿中心线的弧长不变；②横截面始终与中心线正交。

在以上假设条件下，软管形态如图1所示。

1.2 柔性体维修的特点

图1 软管空间形态假设

虚拟现实环境下，柔性体维修与刚性体维修有很大的区别，主要表现在以下两方面：①在描述虚拟维修对象的动态变化时，刚体建模只考虑其在空间内的位移，不考虑其形状改变；而柔性体维修过程中，为满足维修过程逼真度要求，需要实时更新柔性体的形态参数，并进行动态绘制。②虚拟现实环境中，刚性体模型形状固定，碰撞检测易于实现，柔性体在维修过程中存在弯曲、扭转等情况，碰撞干涉难以预测［10］。

2 软管模型的建立与平衡状态分析

2.1 软管参数化模型

基于软管几何特性和物理特性假设，将软管的变形过程简化为其中心线上节点受力前后的位置改变过程。具体描述如下：软管中心线上任一运动节点P的位置可由一矢径r表示，P点弧坐标为s，在t＝0时刻，柔性体的初始状态处于P点位置，如图2所示。

图2 截面中心变形状态描述

主矢r为弧坐标s的单值可微函数，曲线在P点处的单位切向量T可描述如下：

经历时间t1后，质点P移动至P′，该状态为变形状态，P′点位置由矢径r1表示，r1＝r0＋Δr，因此，由矢径r可以完全确定软管中心线上各节点的空间位置，根据式（1），可由下式求解矢径r：

式中，r0为起始点矢径。

2.2 约束关系模型

2.2.1软管的Kirchhoff方程

本文应用Kirchhoff动力学比拟理论［6］求解轴线节点的动平衡问题。首先将节点所受到的外力转换成适当的力和力矩，在定参考坐标系Oξηζ中，考虑微元弧PP′的内力平衡，设P点所在截面受临近截面作用的主矢和主矩分别为－F和－M，P′点所在截面受临近截面作用的主矢和主矩为F＋ΔF和M＋ΔM，如图3所示。

图3 软管微元弧的动力平衡

在受力平衡状态下，上述作用力对P点进行简化后的主矢和主矩均为0，仅保留主矢F和主矩M 增量的一阶小量时，可以得到［6］

其中，0为零向量，将式（3）各项除以Δs，令Δs→0，导出微元弧力平衡与力矩平衡微分方程：

建立P点所在截面的主坐标系Pxyz，根据柔性体几何特性假设，z轴与P点切向量T重合，建立截面主坐标系基向量e＝ex＋ey＋ez，则P点单位切向量T＝ez。为便于求解，将求解过程由相对虚拟环境中定坐标系Oξηζ转换成相对截面的主坐标系Pxyz中进行，将式（4）改写为相对于截面主坐标系Pxyz的投影，即

上式中的“～”表示主矢F与主矩M 相对截面主坐标系Pxyz的局部导数，ω为主坐标系Pxyz相对于定坐标系Oξηζ的角位移变化率，它是与软管弯曲变形和扭转变形相关的矢量。角位移变化率ω、主矢F和主矩M相对于截面主坐标系Pxyz的投影为

根据式（6），可将式（5）改写成如下形式：

为便于计算，将定坐标系中Oζ轴的方向选取为与主矢F方向同向，若将F相对坐标系Pxyz各轴的方向余弦记为α、β、γ，则有：F1＝Fα，F2＝F β，F3＝Fγ，其中，F＝｜F｜。设软管的初始曲率和扭矩均为0，则截面主矩M 可表示为：M1＝Aω1，M2＝Bω2，M3＝Cω3。其中，参数A、B为截面对于x轴和y轴的抗弯刚度，参数C为截面对于z轴的抗扭刚度，参数A、B、C由变形软管的弹性系数和截面的几何形状确定。其中，A＝EIx，B＝EIy，C＝GIz。软管为均匀各向同性时，有G＝E／［2（1＋ν）］，其中，E、G、ν分别为软管的弹性模量、剪切模量和泊松比。Ix、Iy分别为截面相对x轴和y轴的惯性矩，Iz为截面相对z轴的极惯性矩，对于半径为a的圆截面，Ix＝Iy＝ πa4／4，Iy＝πa4／2。

结合上述分析，在忽略了管件的原始曲率和扭率时，可将软管动力平衡方程（式（7）、式（8））转化为如下形式：

常微分方程组（式（9）、式（10））为求解软管动力学平衡的Kirchhoff方程。Kirchhoff方程中的拉 伸 变 量 α（s）、β（s）、γ（s）及 扭 转 变 量 ω1（s）、ω2（s）、ω3（s）决定刚性截面在空间坐标系Oξηζ 中的姿态。

2.2.2欧拉参数形式的Kirchhoff方程

常微分方程组（式（9）、式（10））中，需要求解6个未知数：α、β、γ、ω1、ω2、ω3，本文采用欧拉四元参数法（λ1，λ2，λ3，λ4）将未知数进行统一并对方程进行求解［11］。刚体有限转动的欧拉定理表明：截面绕定点O的任意有限转动可由绕O点的某个轴p的一次有限转动实现，如图4所示。

图4 刚体的有限转动

其中，θ为截面中P点绕Op轴旋转至P′点的旋转角，用以下参数定义欧拉参数：

式中，p1、p2、p3为转动轴Op 相对定坐标系Oξηζ 的方向余弦。

根据欧拉四元参数定义，参数的平方和恒等于1，即

设Qk为λk对弧坐标s的导数，得到

运用刚体无限小转动理论，将横截面角位移变化率ω与欧拉四元参数及其导数的关系表达如下：

令图3中的微元弧PP′长度趋近于零，则P与P′点接近重合。在定坐标系Oξηζ中，忽略软管自身重力，将主矢F向O点简化，得到作用于O点等效力FO，如图5所示。

图5 主矢的等效力表示

设定主坐标系ζ轴与主矢F方向共线，则可确定F1、F2、F3与欧拉参数的关系，表示如下：

桌面式虚拟现实环境中，可设置作用于软管活动端截面上的主矢F和主矩M，因此，可以采用有限差分法将常微分方程组（式（14）、式（15））转化为代数公式进行计算，求解每一截面上的8个未知参数Q1、Q2、Q3、Q4、λ1、λ2、λ3、λ4。可将起始位置（即活动端截面）作为迭代第一步进行求解，将每次迭代所得到的解向量作为下一步迭代的初始向量，固定端位置（参考坐标系中心）作为迭代的最终精确解，依据式（11）～式（13）对每一步迭代所得的解进行修正，可求得各截面微元中心P点的矢径r，以确定活动软管的空间姿态。

设软管长度为L，图5所示的截面中心P点矢径为r＝ （ξ，η，ζ），依据Kirchhoff方程的欧拉角表示方法［6］，矢径r求解方程为

将长度为L的软管均分为n段，每段弧长为s，s＝L／n，第i个节点的矢径为ri＝ （ξ（Pi），η（Pi），ζ（Pi）），设r0＝ （0，0，0），则软管长度约束方程为

2.3 碰撞检测

在装配体中，软管、线缆等柔性体往往不会单独存在，特别是在大型机械设备中，大多数管路走向错综复杂，因此，柔性体之间的碰撞干涉情况是必须考虑的问题。

根据软管模型假设，软管横截面始终为刚性平面，且变形过程中忽略中心线的拉伸变形，所以可以通过检验柔性体之间的最小距离来作为碰撞检测的判断依据，以此来检查柔性体是否满足装配规范。

（1）对于共面或异面的两条管线，可直接检测轴线上节点间的距离D，距离D应满足如下条件：

式中，min｜ri－rj｜为管路节点之间的最小距离；ai＋aj为两管路的半径和。

找出节点间的最短距离Dmin，进行碰撞检测分析。

（2）对于柔性管件的自碰撞问题，可采用球形包围盒方法进行求解，节点间距离d应满足如下条件：

其中，ri＋rj为与节点i和节点j对应的球形包围盒的半径和，对于dmin＜ri＋rj的两个节点，需进行碰撞检测分析。

3 实例分析

以天然橡胶软管（其几何和物理参数如表1所示）为例，将橡胶软管进行n等份划分，构成节点集P＝ ｛P0，P1，…，Pn｝，需对集合内Pi点的位置参数ξ（Pi）、η（Pi）、ζ（Pi）进行求解，形成软管空间形态。

表1 天然橡胶软管几何和物理参数

初始条件设定如下：忽略柔性软管的重力及初始弯扭矩，软管运动端受恒力F作用（｜F｜＝100N）。a→b→c→d即为运动端节点的受力方向，设定ab段、bc段、cd段直线方程如下：

依据关键帧记录的橡胶软管运动端节点的空间位置，已知软管长度L＝2m，均分为n＝40等份。初始状态下，活动端r0＝ （0，0，0），自由端r40＝（0，0，2），自由端在主矢F作用下的移动过程中，对软管40个运动节点的平衡状态方程进行差分运算，确定整根软管的空间姿态。

假设软管受拉伸时受力均匀，终止时刻软管固定端矢径r0＝（ξ0，η0，ζ0）＝（0，0，0），自由端节点矢径r40＝ （ξ40，η40，ζ40）＝ （0.8，1.6，0），欧拉四 元 参 数 （λ1，0，λ2，0，λ3，0，λ4，0）＝ （1，0，0，0），（λ1，40，λ2，40，λ3，40，λ4，40）＝ （0.8507，0.2351，0.4702，0），使用有限差分方法求解柔性软管各个节点的解向量，将节点依次连接，采用均匀B样条曲线拟合法显示在输出设备上，得到软管轴线的形态仿真结果，如图6所示。

图6 软管中轴线运动仿真

研究中发现，采用均匀B样条对各节点间曲线进行拟合，达到了曲线形状拟实性要求。但如果节点数目过少，会导致曲线长度发生变化，远远小于曲线的初始长度；而节点数过多，会加重计算负担，影响计算机的实时绘制效果，因此设计过程中必须对节点数目进行控制。选用具有表1所示特性的天然橡胶软管进行对比分析（实测图如图7所示），为避免重力因素影响，采用水平方向受力的方式进行实验，将实测数据与仿真结果进行对比后发现，天然橡胶软管的运动趋势与软管轴线的形态仿真数据基本吻合，但由于受摩擦因素和弯扭度非线性变化影响，误差依然存在。

图7 作用力F下软管变形实测结果

4 结论

（1）针对虚拟维修过程中柔性管件维修难的问题，构建了基于刚性截面假设的柔性体空间形态模型，该模型可以对柔性体截面形状进行修改（如圆形截面的软管、方形截面的皮带等），满足虚拟维修环境下柔性体形态多样性的需要。

（2）引入Kirchhoff动力平衡理论，在定参考坐标系Oξηζ和动坐标系Pxyz中对横截面的运动状态进行模拟，并构建与之匹配的平衡状态方程，使用欧拉四元参数法对方程进行求解，该方法能够实时求解出柔性体形态参数，实现虚拟维修环境下柔性管件维修可视化。

（3）基于刚性截面理论假设，提出了基于该假设的碰撞检测方法，该方法只需要考虑节点划分及截面形态等因素，可以较为准确地判断出软管间的碰撞干涉情况。下一步将对空间网格模型假设的不规则柔性体模型进行研究，解决刚性截面假设只能描述细长管件模型的问题。

［1］刘佳，刘毅.虚拟维修技术发展综述［J］.计算机辅助设计与图形学学报，2009，21（11）：1520－1534.Liu Jia，Liu Yi.A Survey of Virtual Maintenance Technology［J］.Journal of Computer－aided Design and Computer Graphics，2009，21（11）：1520－1534.

［2］Ehud K，Jan W.Toward Assembly Sequence Planning with Flexible Parts［C］／／Proceeding of IEEE International Conference on Robotics and Automation.Minneapolis，USA：IEEE Press，1996：517－524.

［3］Loock A，Schömer E.A Virtual Environment for Interactive Assembly Simulation：From Rigid Bodies to Deformable Cable［C］／／5th World Multiconference on Systemics，Cybernetics and Informatics.Florida，USA：Springer Press，2001：325－332.

［4］刘检华，万毕乐，宁汝新.虚拟环境下基于离散控制点的线缆装配规划技术［J］.机械工程学报，2006，42（8）：125－130.Liu Jianhua，Wan Bile，Ning Ruxin.Realization Technology of Cable Harness Process Planning in Virtual Environment Based on Discrete Control Point Modeling Method［J］.Journal of Mechanical Engineering，2006，42（8）：125－130.

［5］万毕乐，宁汝新，刘检华.面向交叉装配的线缆虚拟装配信息规划信息管理技术［J］.中国机械工程，2008，19（7）：793－797.Wan Bile，Ning Ruxin，Liu Jianhua.Information Management Technology of Cable Harness Virtual Assembly Planning for Cross Assembly［J］.China Mechanical Engineering，2008，19（7）：793－797.

［6］刘延柱.弹性细杆的非线性力学［M］.北京：清华大学出版社，2003.

［7］薛纭，刘延柱.Kirchhoff弹性杆分析动力学的准坐标表达［J］.力学季刊，2006，27（4）：550－556.Xue Yun，Liu Yanzhu.Analytical Dynamics of a Kirchhoff Elastic Rod Expressed in Quasi－coordinates［J］.Chinese Quarterly of Mechanics，2006，27（4）：550－556.

［8］翁玉权.超细长弹性圆截面杆模型的Noether对称性、Lie对称性及其守恒量［D］.青岛：青岛大学，2008.

［9］刘检华，赵涛，王春生，等.虚拟环境下的活动线缆物理特性建模与运动仿真技术［J］.机械工程学报，2011，47（9）：117－124.Liu Jianhua，Zhao Tao，Wang Chunsheng，et al.Motion Cable Harness Physical Characteristic Oriented Modeling and Kinetic Simulation Technology in Virtual Environment［J］.Journal of Mechanical Engineering，2011，47（9）：117－124.

［10］刘鹏远，李瑞华，胡昌林.虚拟环境下柔性线缆的建模方法［J］.军械工程学院学报，2008，20（3）：63－65.Liu Pengyuan，Li Ruihua，Hu Changlin.Flexible Cable Modeling Method in Virtual Environment［J］.Journal of Ordnance Engineering College，2008，20（3）：63－65.

［11］刘延柱.关于刚体姿态的数学表达［J］.力学与实践，2008，30（1）：98－101.Liu Yanzhu.Math Expression of Pose of Rigid Body［J］.Mechanics in Engineering，2008，30（1）：98－101.