焦玉民 王 強(qiáng) 徐 婷 謝慶華

1.94679部隊(duì)，南京，210038 2.解放軍理工大學(xué)，南京，210007

0 引言

軍用工程機(jī)械防護(hù)性能最薄弱的部位是暴露在機(jī)械外部的液壓管路，液壓管路一旦損壞，極易造成機(jī)械無法工作。虛擬維修訓(xùn)練是軍用工程機(jī)械維修保障訓(xùn)練的重要手段之一，柔性零件（如軟管、線纜、皮帶等）的維修模擬是一個不容回避的問題。目前的虛擬維修仿真過程中，絕大多數(shù)維修對象均為剛性體（如殼體、齒輪、軸承），仿真過程僅考慮其在虛擬環(huán)境中的位置變化，不考慮其受力后的變形過程。而變形體的虛擬維修過程涉及單元類型、材料特性和外力載荷等問題，建模技術(shù)復(fù)雜。目前的虛擬現(xiàn)實(shí)環(huán)境中主要通過B樣條描述柔性體的靜態(tài)形狀，采用顏色變換、聲音提示以及動畫等方式模擬其故障狀態(tài)，難以模擬柔性體維修變形過程［1］。國內(nèi)外學(xué)者在此方面開展了一系列研究。文獻(xiàn)［2］采用了離散點(diǎn)方法來構(gòu)建曲線模型，該模型在幾何層面可以較好地模擬出曲線形態(tài)，但由于缺少曲線的物理特征描述，難以對不同類型的柔性體進(jìn)行區(qū)分。文獻(xiàn)［3］提出了基于彈簧－質(zhì)點(diǎn)法的線纜彈性模型，使用彈簧的彈性系數(shù)描述線纜的物理特性，為柔性軟管的運(yùn)動求解提供了一種有效的解決辦法，但彈簧－質(zhì)點(diǎn)法依然存在局限性，該方法僅能計算軟管的彎曲變形過程，無法對軟管的扭轉(zhuǎn)變形進(jìn)行表達(dá)。劉檢華等［4］、萬畢樂等［5］對虛擬現(xiàn)實(shí)環(huán)境下線纜的建模與裝配進(jìn)行了研究，提出了離散控制點(diǎn)建模方法，對線纜的幾何約束關(guān)系和裝配方法進(jìn)行了詳細(xì)的分析，但并沒有提及在外力影響下線纜的變形規(guī)律。劉延柱［6］和薛紜等［7］以DNA力學(xué)模型為參考，提出了基于Kirchhoff動力學(xué)比擬理論的彈性桿動力學(xué)模型，為虛擬現(xiàn)實(shí)環(huán)境下柔性體建模與分析提供了一種有效的解決方法。文獻(xiàn)［8－9］對Kirchhoff動力學(xué)比擬理論進(jìn)行了深入研究，并應(yīng)用于不同的研究領(lǐng)域。

在上述研究的基礎(chǔ)上，本文提出了虛擬現(xiàn)實(shí)環(huán)境下柔性管件空間形態(tài)的基本假設(shè)，并運(yùn)用Kirchhoff動力學(xué)比擬理論對柔性體受力后的平衡狀態(tài)進(jìn)行建模分析，構(gòu)建了柔性體動力學(xué)平衡方程，使用歐拉四元參數(shù)法對平衡方程進(jìn)行求解，解決了虛擬現(xiàn)實(shí)環(huán)境下柔性體受力后的動態(tài)變形可視化問題，為在虛擬維修環(huán)境中實(shí)現(xiàn)柔性管件的布局優(yōu)化以及維修過程的動態(tài)模擬等問題提供了切實(shí)有效的解決方案。

1 虛擬現(xiàn)實(shí)環(huán)境下柔性管件的基本假設(shè)

1.1 柔性管件描述

柔性體在外力作用下經(jīng)受拉伸變形和扭轉(zhuǎn)變形，其連接狀態(tài)、材料性質(zhì)和外力載荷都影響著變形過程，本文將軟管受力后的變形簡化為軟管中心線的幾何形狀改變，為便于對柔性管件變形狀態(tài)進(jìn)行模擬計算，需對管件模型作出如下假設(shè)。幾何特性方面：①將軟管均分為n段，以每段軟管端部橫截面作為研究對象；②端部橫截面均為剛性平面；③相鄰截面可繞中心線做相對轉(zhuǎn)動。物理特性方面：①忽略軟管中心線的拉伸變形，任意兩截面沿中心線的弧長不變；②橫截面始終與中心線正交。

在以上假設(shè)條件下，軟管形態(tài)如圖1所示。

1.2 柔性體維修的特點(diǎn)

圖1 軟管空間形態(tài)假設(shè)

虛擬現(xiàn)實(shí)環(huán)境下，柔性體維修與剛性體維修有很大的區(qū)別，主要表現(xiàn)在以下兩方面：①在描述虛擬維修對象的動態(tài)變化時，剛體建模只考慮其在空間內(nèi)的位移，不考慮其形狀改變；而柔性體維修過程中，為滿足維修過程逼真度要求，需要實(shí)時更新柔性體的形態(tài)參數(shù)，并進(jìn)行動態(tài)繪制。②虛擬現(xiàn)實(shí)環(huán)境中，剛性體模型形狀固定，碰撞檢測易于實(shí)現(xiàn)，柔性體在維修過程中存在彎曲、扭轉(zhuǎn)等情況，碰撞干涉難以預(yù)測［10］。

2 軟管模型的建立與平衡狀態(tài)分析

2.1 軟管參數(shù)化模型

基于軟管幾何特性和物理特性假設(shè)，將軟管的變形過程簡化為其中心線上節(jié)點(diǎn)受力前后的位置改變過程。具體描述如下：軟管中心線上任一運(yùn)動節(jié)點(diǎn)P的位置可由一矢徑r表示，P點(diǎn)弧坐標(biāo)為s，在t＝0時刻，柔性體的初始狀態(tài)處于P點(diǎn)位置，如圖2所示。

圖2 截面中心變形狀態(tài)描述

主矢r為弧坐標(biāo)s的單值可微函數(shù)，曲線在P點(diǎn)處的單位切向量T可描述如下：

經(jīng)歷時間t1后，質(zhì)點(diǎn)P移動至P′，該狀態(tài)為變形狀態(tài)，P′點(diǎn)位置由矢徑r1表示，r1＝r0＋Δr，因此，由矢徑r可以完全確定軟管中心線上各節(jié)點(diǎn)的空間位置，根據(jù)式（1），可由下式求解矢徑r：

式中，r0為起始點(diǎn)矢徑。

2.2 約束關(guān)系模型

2.2.1軟管的Kirchhoff方程

本文應(yīng)用Kirchhoff動力學(xué)比擬理論［6］求解軸線節(jié)點(diǎn)的動平衡問題。首先將節(jié)點(diǎn)所受到的外力轉(zhuǎn)換成適當(dāng)?shù)牧土?，在定參考坐?biāo)系Oξηζ中，考慮微元弧PP′的內(nèi)力平衡，設(shè)P點(diǎn)所在截面受臨近截面作用的主矢和主矩分別為－F和－M，P′點(diǎn)所在截面受臨近截面作用的主矢和主矩為F＋ΔF和M＋ΔM，如圖3所示。

圖3 軟管微元弧的動力平衡

在受力平衡狀態(tài)下，上述作用力對P點(diǎn)進(jìn)行簡化后的主矢和主矩均為0，僅保留主矢F和主矩M 增量的一階小量時，可以得到［6］

其中，0為零向量，將式（3）各項(xiàng)除以Δs，令Δs→0，導(dǎo)出微元弧力平衡與力矩平衡微分方程：

建立P點(diǎn)所在截面的主坐標(biāo)系Pxyz，根據(jù)柔性體幾何特性假設(shè)，z軸與P點(diǎn)切向量T重合，建立截面主坐標(biāo)系基向量e＝ex＋ey＋ez，則P點(diǎn)單位切向量T＝ez。為便于求解，將求解過程由相對虛擬環(huán)境中定坐標(biāo)系Oξηζ轉(zhuǎn)換成相對截面的主坐標(biāo)系Pxyz中進(jìn)行，將式（4）改寫為相對于截面主坐標(biāo)系Pxyz的投影，即

上式中的“～”表示主矢F與主矩M 相對截面主坐標(biāo)系Pxyz的局部導(dǎo)數(shù)，ω為主坐標(biāo)系Pxyz相對于定坐標(biāo)系Oξηζ的角位移變化率，它是與軟管彎曲變形和扭轉(zhuǎn)變形相關(guān)的矢量。角位移變化率ω、主矢F和主矩M相對于截面主坐標(biāo)系Pxyz的投影為

根據(jù)式（6），可將式（5）改寫成如下形式：

為便于計算，將定坐標(biāo)系中Oζ軸的方向選取為與主矢F方向同向，若將F相對坐標(biāo)系Pxyz各軸的方向余弦記為α、β、γ，則有：F1＝Fα，F(xiàn)2＝F β，F(xiàn)3＝Fγ，其中，F(xiàn)＝｜F｜。設(shè)軟管的初始曲率和扭矩均為0，則截面主矩M 可表示為：M1＝Aω1，M2＝Bω2，M3＝Cω3。其中，參數(shù)A、B為截面對于x軸和y軸的抗彎剛度，參數(shù)C為截面對于z軸的抗扭剛度，參數(shù)A、B、C由變形軟管的彈性系數(shù)和截面的幾何形狀確定。其中，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。

結(jié)合上述分析，在忽略了管件的原始曲率和扭率時，可將軟管動力平衡方程（式（7）、式（8））轉(zhuǎn)化為如下形式：

常微分方程組（式（9）、式（10））為求解軟管動力學(xué)平衡的Kirchhoff方程。Kirchhoff方程中的拉 伸 變 量 α（s）、β（s）、γ（s）及 扭 轉(zhuǎn) 變 量 ω1（s）、ω2（s）、ω3（s）決定剛性截面在空間坐標(biāo)系Oξηζ 中的姿態(tài)。

2.2.2歐拉參數(shù)形式的Kirchhoff方程

常微分方程組（式（9）、式（10））中，需要求解6個未知數(shù)：α、β、γ、ω1、ω2、ω3，本文采用歐拉四元參數(shù)法（λ1，λ2，λ3，λ4）將未知數(shù)進(jìn)行統(tǒng)一并對方程進(jìn)行求解［11］。剛體有限轉(zhuǎn)動的歐拉定理表明：截面繞定點(diǎn)O的任意有限轉(zhuǎn)動可由繞O點(diǎn)的某個軸p的一次有限轉(zhuǎn)動實(shí)現(xiàn)，如圖4所示。

圖4 剛體的有限轉(zhuǎn)動

其中，θ為截面中P點(diǎn)繞Op軸旋轉(zhuǎn)至P′點(diǎn)的旋轉(zhuǎn)角，用以下參數(shù)定義歐拉參數(shù)：

式中，p1、p2、p3為轉(zhuǎn)動軸Op 相對定坐標(biāo)系Oξηζ 的方向余弦。

根據(jù)歐拉四元參數(shù)定義，參數(shù)的平方和恒等于1，即

設(shè)Qk為λk對弧坐標(biāo)s的導(dǎo)數(shù)，得到

運(yùn)用剛體無限小轉(zhuǎn)動理論，將橫截面角位移變化率ω與歐拉四元參數(shù)及其導(dǎo)數(shù)的關(guān)系表達(dá)如下：

令圖3中的微元弧PP′長度趨近于零，則P與P′點(diǎn)接近重合。在定坐標(biāo)系Oξηζ中，忽略軟管自身重力，將主矢F向O點(diǎn)簡化，得到作用于O點(diǎn)等效力FO，如圖5所示。

圖5 主矢的等效力表示

設(shè)定主坐標(biāo)系ζ軸與主矢F方向共線，則可確定F1、F2、F3與歐拉參數(shù)的關(guān)系，表示如下：

桌面式虛擬現(xiàn)實(shí)環(huán)境中，可設(shè)置作用于軟管活動端截面上的主矢F和主矩M，因此，可以采用有限差分法將常微分方程組（式（14）、式（15））轉(zhuǎn)化為代數(shù)公式進(jìn)行計算，求解每一截面上的8個未知參數(shù)Q1、Q2、Q3、Q4、λ1、λ2、λ3、λ4?？蓪⑵鹗嘉恢茫椿顒佣私孛妫┳鳛榈谝徊竭M(jìn)行求解，將每次迭代所得到的解向量作為下一步迭代的初始向量，固定端位置（參考坐標(biāo)系中心）作為迭代的最終精確解，依據(jù)式（11）～式（13）對每一步迭代所得的解進(jìn)行修正，可求得各截面微元中心P點(diǎn)的矢徑r，以確定活動軟管的空間姿態(tài)。

設(shè)軟管長度為L，圖5所示的截面中心P點(diǎn)矢徑為r＝ （ξ，η，ζ），依據(jù)Kirchhoff方程的歐拉角表示方法［6］，矢徑r求解方程為

將長度為L的軟管均分為n段，每段弧長為s，s＝L／n，第i個節(jié)點(diǎn)的矢徑為ri＝ （ξ（Pi），η（Pi），ζ（Pi）），設(shè)r0＝ （0，0，0），則軟管長度約束方程為

2.3 碰撞檢測

在裝配體中，軟管、線纜等柔性體往往不會單獨(dú)存在，特別是在大型機(jī)械設(shè)備中，大多數(shù)管路走向錯綜復(fù)雜，因此，柔性體之間的碰撞干涉情況是必須考慮的問題。

根據(jù)軟管模型假設(shè)，軟管橫截面始終為剛性平面，且變形過程中忽略中心線的拉伸變形，所以可以通過檢驗(yàn)柔性體之間的最小距離來作為碰撞檢測的判斷依據(jù)，以此來檢查柔性體是否滿足裝配規(guī)范。

（1）對于共面或異面的兩條管線，可直接檢測軸線上節(jié)點(diǎn)間的距離D，距離D應(yīng)滿足如下條件：

式中，min｜ri－rj｜為管路節(jié)點(diǎn)之間的最小距離；ai＋aj為兩管路的半徑和。

找出節(jié)點(diǎn)間的最短距離Dmin，進(jìn)行碰撞檢測分析。

（2）對于柔性管件的自碰撞問題，可采用球形包圍盒方法進(jìn)行求解，節(jié)點(diǎn)間距離d應(yīng)滿足如下條件：

其中，ri＋rj為與節(jié)點(diǎn)i和節(jié)點(diǎn)j對應(yīng)的球形包圍盒的半徑和，對于dmin＜ri＋rj的兩個節(jié)點(diǎn)，需進(jìn)行碰撞檢測分析。

3 實(shí)例分析

以天然橡膠軟管（其幾何和物理參數(shù)如表1所示）為例，將橡膠軟管進(jìn)行n等份劃分，構(gòu)成節(jié)點(diǎn)集P＝ ｛P0，P1，…，Pn｝，需對集合內(nèi)Pi點(diǎn)的位置參數(shù)ξ（Pi）、η（Pi）、ζ（Pi）進(jìn)行求解，形成軟管空間形態(tài)。

表1 天然橡膠軟管幾何和物理參數(shù)

初始條件設(shè)定如下：忽略柔性軟管的重力及初始彎扭矩，軟管運(yùn)動端受恒力F作用（｜F｜＝100N）。a→b→c→d即為運(yùn)動端節(jié)點(diǎn)的受力方向，設(shè)定ab段、bc段、cd段直線方程如下：

依據(jù)關(guān)鍵幀記錄的橡膠軟管運(yùn)動端節(jié)點(diǎn)的空間位置，已知軟管長度L＝2m，均分為n＝40等份。初始狀態(tài)下，活動端r0＝ （0，0，0），自由端r40＝（0，0，2），自由端在主矢F作用下的移動過程中，對軟管40個運(yùn)動節(jié)點(diǎn)的平衡狀態(tài)方程進(jìn)行差分運(yùn)算，確定整根軟管的空間姿態(tài)。

假設(shè)軟管受拉伸時受力均勻，終止時刻軟管固定端矢徑r0＝（ξ0，η0，ζ0）＝（0，0，0），自由端節(jié)點(diǎn)矢徑r40＝ （ξ40，η40，ζ40）＝ （0.8，1.6，0），歐拉四 元 參 數(shù) （λ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），使用有限差分方法求解柔性軟管各個節(jié)點(diǎn)的解向量，將節(jié)點(diǎn)依次連接，采用均勻B樣條曲線擬合法顯示在輸出設(shè)備上，得到軟管軸線的形態(tài)仿真結(jié)果，如圖6所示。

圖6 軟管中軸線運(yùn)動仿真

研究中發(fā)現(xiàn)，采用均勻B樣條對各節(jié)點(diǎn)間曲線進(jìn)行擬合，達(dá)到了曲線形狀擬實(shí)性要求。但如果節(jié)點(diǎn)數(shù)目過少，會導(dǎo)致曲線長度發(fā)生變化，遠(yuǎn)遠(yuǎn)小于曲線的初始長度；而節(jié)點(diǎn)數(shù)過多，會加重計算負(fù)擔(dān)，影響計算機(jī)的實(shí)時繪制效果，因此設(shè)計過程中必須對節(jié)點(diǎn)數(shù)目進(jìn)行控制。選用具有表1所示特性的天然橡膠軟管進(jìn)行對比分析（實(shí)測圖如圖7所示），為避免重力因素影響，采用水平方向受力的方式進(jìn)行實(shí)驗(yàn)，將實(shí)測數(shù)據(jù)與仿真結(jié)果進(jìn)行對比后發(fā)現(xiàn)，天然橡膠軟管的運(yùn)動趨勢與軟管軸線的形態(tài)仿真數(shù)據(jù)基本吻合，但由于受摩擦因素和彎扭度非線性變化影響，誤差依然存在。

圖7 作用力F下軟管變形實(shí)測結(jié)果

4 結(jié)論

（1）針對虛擬維修過程中柔性管件維修難的問題，構(gòu)建了基于剛性截面假設(shè)的柔性體空間形態(tài)模型，該模型可以對柔性體截面形狀進(jìn)行修改（如圓形截面的軟管、方形截面的皮帶等），滿足虛擬維修環(huán)境下柔性體形態(tài)多樣性的需要。

（2）引入Kirchhoff動力平衡理論，在定參考坐標(biāo)系Oξηζ和動坐標(biāo)系Pxyz中對橫截面的運(yùn)動狀態(tài)進(jìn)行模擬，并構(gòu)建與之匹配的平衡狀態(tài)方程，使用歐拉四元參數(shù)法對方程進(jìn)行求解，該方法能夠?qū)崟r求解出柔性體形態(tài)參數(shù)，實(shí)現(xiàn)虛擬維修環(huán)境下柔性管件維修可視化。

（3）基于剛性截面理論假設(shè)，提出了基于該假設(shè)的碰撞檢測方法，該方法只需要考慮節(jié)點(diǎn)劃分及截面形態(tài)等因素，可以較為準(zhǔn)確地判斷出軟管間的碰撞干涉情況。下一步將對空間網(wǎng)格模型假設(shè)的不規(guī)則柔性體模型進(jìn)行研究，解決剛性截面假設(shè)只能描述細(xì)長管件模型的問題。

［1］劉佳，劉毅.虛擬維修技術(shù)發(fā)展綜述［J］.計算機(jī)輔助設(shè)計與圖形學(xué)學(xué)報，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］劉檢華，萬畢樂，寧汝新.虛擬環(huán)境下基于離散控制點(diǎn)的線纜裝配規(guī)劃技術(shù)［J］.機(jī)械工程學(xué)報，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］萬畢樂，寧汝新，劉檢華.面向交叉裝配的線纜虛擬裝配信息規(guī)劃信息管理技術(shù)［J］.中國機(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］劉延柱.彈性細(xì)桿的非線性力學(xué)［M］.北京：清華大學(xué)出版社，2003.

［7］薛紜，劉延柱.Kirchhoff彈性桿分析動力學(xué)的準(zhǔn)坐標(biāo)表達(dá)［J］.力學(xué)季刊，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］翁玉權(quán).超細(xì)長彈性圓截面桿模型的Noether對稱性、Lie對稱性及其守恒量［D］.青島：青島大學(xué)，2008.

［9］劉檢華，趙濤，王春生，等.虛擬環(huán)境下的活動線纜物理特性建模與運(yùn)動仿真技術(shù)［J］.機(jī)械工程學(xué)報，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］劉鵬遠(yuǎn)，李瑞華，胡昌林.虛擬環(huán)境下柔性線纜的建模方法［J］.軍械工程學(xué)院學(xué)報，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］劉延柱.關(guān)于剛體姿態(tài)的數(shù)學(xué)表達(dá)［J］.力學(xué)與實(shí)踐，2008，30（1）：98－101.Liu Yanzhu.Math Expression of Pose of Rigid Body［J］.Mechanics in Engineering，2008，30（1）：98－101.