基于Cosserat彈性桿理論的柔性線纜物理建模方法

林海立， 劉檢華， 唐承統(tǒng)， 劉佳順

(北京理工大學(xué)機(jī)械與車輛學(xué)院，北京 100081)

基于Cosserat彈性桿理論的柔性線纜物理建模方法

林海立， 劉檢華， 唐承統(tǒng)， 劉佳順

(北京理工大學(xué)機(jī)械與車輛學(xué)院，北京 100081)

針對(duì)虛擬環(huán)境中線纜敷設(shè)過(guò)程仿真的線纜建模問(wèn)題，提出一種基于Cosserat彈性桿理論的柔性線纜物理建模方法。該方法在線纜離散模型的基礎(chǔ)上根據(jù)Cosserat彈性桿理論導(dǎo)出線纜的總體勢(shì)能，并通過(guò)優(yōu)化算法得到線纜總體勢(shì)能最小時(shí)的線纜形態(tài)，即線纜的平衡狀態(tài)。模型考慮到線纜的柔性和連續(xù)性，能夠模擬線纜的彎曲和扭轉(zhuǎn)變形；通過(guò)對(duì)線纜進(jìn)行離散表達(dá)，實(shí)現(xiàn)了對(duì)線纜總體勢(shì)能的數(shù)值求解；引入罰函數(shù)將求解線纜總體勢(shì)能最小值問(wèn)題轉(zhuǎn)化為非線性無(wú)約束優(yōu)化問(wèn)題；非線性無(wú)約束優(yōu)化問(wèn)題的求解采用信賴域方法，提高了求解的穩(wěn)定性。設(shè)計(jì)并開發(fā)了虛擬環(huán)境中的線纜敷設(shè)過(guò)程仿真原型系統(tǒng)，并進(jìn)行實(shí)例驗(yàn)證，證明該方法的可行性。

線纜；物理建模；敷設(shè)仿真；Cosserat彈性桿理論；能量最小

作為能量和信號(hào)的傳輸通道，線纜在機(jī)電產(chǎn)品中應(yīng)用廣泛，而線纜的布局和敷設(shè)安裝質(zhì)量直接影響機(jī)電產(chǎn)品的質(zhì)量及可靠性。近年來(lái)，隨著CAD技術(shù)的發(fā)展，在虛擬環(huán)境中進(jìn)行線纜布局設(shè)計(jì)與敷設(shè)過(guò)程仿真逐漸成為國(guó)內(nèi)外的研究熱點(diǎn)。

在工程應(yīng)用中，線纜具有復(fù)雜的拓?fù)浣Y(jié)構(gòu)，涵蓋了大量的幾何拓?fù)湫畔?、工程語(yǔ)義信息和分散的布線信息[1]。柔性線纜的模型表達(dá)是進(jìn)行虛擬環(huán)境下線纜敷設(shè)過(guò)程仿真的基礎(chǔ)，該模型不僅要支持電纜的幾何和拓?fù)涮匦?，還要支持線纜的物理特性表達(dá)，并滿足虛擬現(xiàn)實(shí)系統(tǒng)特有的實(shí)時(shí)性要求。與布料和人體組織等柔性體不同，線纜最大的特點(diǎn)是其長(zhǎng)度通常遠(yuǎn)大于直徑，工程中通常以線纜的中心線為掃描路徑，以一定的截面信息掃描來(lái)形成線纜的三維模型，因此中心線的確定是線纜建模的關(guān)鍵，先前表示線纜的中心線的方法包括折線段，或自由曲線如貝塞爾曲線、B樣條曲線[2-3]等，但這些方法屬于線纜的幾何建模，沒(méi)有考慮線纜的物理屬性，進(jìn)行線纜敷設(shè)過(guò)程仿真時(shí)缺乏真實(shí)性。

柔性物體建模與仿真研究一直受到國(guó)內(nèi)外學(xué)者的關(guān)注，IEEE(computer graphics and applications) 與SIGGRAPH曾數(shù)次召開專門的小組討論會(huì)對(duì)該方向進(jìn)行探討[4]。由于柔性物體的形狀不像剛體那樣固定不變，而是隨著外界條件的不同發(fā)生改變，因此建模難度很大。線纜物理建模是典型的柔性體的物理建模，由于線纜敷設(shè)過(guò)程中，其整體的運(yùn)動(dòng)與機(jī)器人機(jī)械臂的運(yùn)動(dòng)具有相似性，基于該思想，魏發(fā)遠(yuǎn)等[5]提出了一種基于蛇形機(jī)器人的線纜建模與敷設(shè)仿真方法，將來(lái)源于機(jī)器人理論的逆運(yùn)動(dòng)學(xué)方法應(yīng)用于線纜建模中。德國(guó)學(xué)者Hergenr?ther 和D?hne[6]開發(fā)的虛擬裝配系統(tǒng)中采用了類似的模型，通過(guò)一系列相連接的剛性桿表示線纜，每個(gè)連接處稱為關(guān)節(jié)，在關(guān)節(jié)處加入扭簧和質(zhì)點(diǎn)以考慮抗彎剛度和重力這兩個(gè)物理屬性，當(dāng)線纜端部移動(dòng)到某一位置時(shí)，采用能量?jī)?yōu)化法和逆運(yùn)動(dòng)學(xué)方法，求解出每一關(guān)節(jié)所應(yīng)有的變量并使系統(tǒng)具有最小能量。以上建模方法能夠體現(xiàn)線纜在敷設(shè)仿真過(guò)程中“長(zhǎng)度不變”的性質(zhì)，但考慮的線纜物理參數(shù)較少，真實(shí)性不足。

能量曲線方法由加拿大學(xué)者 Terzopoulos和Qin[7]提出，其認(rèn)為曲線具有質(zhì)量和抗彎特性，在一定的約束下，當(dāng)曲線具有最小物理變形能時(shí)曲線處于穩(wěn)定狀態(tài)，因此通過(guò)數(shù)學(xué)規(guī)劃或優(yōu)化方法得到具有最小物理變形能的曲線，將該曲線作為線纜的中心線，即可得到處于穩(wěn)定狀態(tài)的線纜。王志斌等[8]在線纜布局設(shè)計(jì)中利用該模型對(duì)線纜的長(zhǎng)度進(jìn)行取樣，解決了線纜取樣長(zhǎng)度不精確的問(wèn)題。

彈簧質(zhì)點(diǎn)模型[9-11]在線纜建模與敷設(shè)仿真中有較多應(yīng)用，該模型具有易于實(shí)施、計(jì)算效率高的優(yōu)點(diǎn)。在該方法中線纜由一系列質(zhì)點(diǎn)和連接質(zhì)點(diǎn)的直線彈簧表示，在質(zhì)點(diǎn)處加入扭簧或者在兩個(gè)質(zhì)點(diǎn)間隔加入直線彈簧，可以考慮線纜的抗彎剛度。在彈簧質(zhì)點(diǎn)模型解算方面，Loock等[12]在Baraff基礎(chǔ)上采用隱式計(jì)算，該算法允許較長(zhǎng)的時(shí)間步，具有較強(qiáng)的實(shí)時(shí)性，但真實(shí)性無(wú)法保證。

由于逆運(yùn)動(dòng)學(xué)、能量曲線和彈簧質(zhì)點(diǎn)模型在處理線纜受約束形變時(shí)部分考慮了線纜的物理屬性，具有一定真實(shí)性且具有實(shí)時(shí)性強(qiáng)的特點(diǎn)，是目前線纜敷設(shè)仿真中主要采用的物理模型。但以上建模方法考慮線纜物理屬性并不全面，其中最重要的是不能對(duì)線纜的扭轉(zhuǎn)變形進(jìn)行模擬，使得在線纜敷設(shè)過(guò)程仿真中其真實(shí)性并不十分理想。

線纜在敷設(shè)過(guò)程中的局部應(yīng)變很小，基本保持在材料的彈性范圍內(nèi)，因此有學(xué)者將其本構(gòu)關(guān)系抽象為彈性關(guān)系，并用彈性桿模型進(jìn)行線纜物理建模。彈性細(xì)桿靜力學(xué)理論是Kirchhoff于1859年建立的，工程中繩索、鉆桿、纖維等都曾將彈性細(xì)桿作為其力學(xué)模型[13]。劉檢華等[14]提出的以Kirchhoff彈性細(xì)桿非線性力學(xué)理論為基礎(chǔ)的虛擬環(huán)境下活動(dòng)線纜建模與運(yùn)動(dòng)仿真方法，對(duì)活動(dòng)線纜的仿真有較好的效果，但該方法需求解微分方程組的邊界值，并利用“打靶法”在已知兩端位置和朝向的基礎(chǔ)上對(duì)線纜的整體外形進(jìn)行求解，其問(wèn)題是求解所需要的初值難以給定，而且求解時(shí)間較長(zhǎng)，無(wú)法滿足線纜敷設(shè)過(guò)程仿真的實(shí)時(shí)性需求。在物理模型求解方面，Bertails等[15]利用離散方法求解Kirchhoff方程，并提出一種離散和連續(xù)相結(jié)合的模型，將發(fā)絲模擬為一系列相連的螺旋曲線段，由通用坐標(biāo)描述線纜，并用積分的方法重構(gòu)線纜，從而消除了加入長(zhǎng)度不變約束產(chǎn)生的硬度，提高了數(shù)值求解的穩(wěn)定性，該方法為目前較為先進(jìn)的彈性桿建模方法。而 Bergou等[16]采用離散微分幾何的方法，顯示以最少自由度的簡(jiǎn)化坐標(biāo)方程表示線纜的中心線，并以標(biāo)量表達(dá)線纜材料框架，以避免使用聯(lián)系中心線和材料框架的硬化約束。

Cosserat理論作為Kirechhoff理論的改進(jìn)，考慮了彈性桿的軸向線應(yīng)變和彎曲剪應(yīng)變等物理參數(shù)，建立更精確的平衡方程。Cosserat理論由Cosserat兄弟于19世紀(jì)提出，目前該方法已經(jīng)成為非線性彈性力學(xué)的基本方法。Cosserat理論首先由 Pai[17]用于彈性桿這類細(xì)長(zhǎng)的可變形結(jié)構(gòu)的建模，該方法同樣需通過(guò)“打靶法”數(shù)值求解微分方程得到桿的靜態(tài)外形；由于數(shù)值解法在初始條件上難以確定，同時(shí)在處理自接觸和與其他物體接觸的問(wèn)題上也存在困難，之后的方法多關(guān)注于模型的離散表達(dá)。Spillmann和Teschner[18]提出用于柔性桿動(dòng)力學(xué)模型，該模型以Cosserat理論為基礎(chǔ)導(dǎo)出桿件的連續(xù)能量表達(dá)式，通過(guò)對(duì)桿的離散計(jì)算出桿的每個(gè)單元具有的能量，進(jìn)而由拉格朗日方程得到桿的動(dòng)力學(xué)微分方程，數(shù)值求解該微分方程得到桿的動(dòng)態(tài)變化。由于動(dòng)力學(xué)模型需要較長(zhǎng)時(shí)間進(jìn)入穩(wěn)定狀態(tài)，黃勁等[19]提出了一個(gè)穩(wěn)定、快速的優(yōu)化策略對(duì)線纜進(jìn)行形變模擬，并著重研究了大步長(zhǎng)、準(zhǔn)靜態(tài)模擬過(guò)程中的接觸處理這一關(guān)鍵問(wèn)題。Grégoire 和 Sch?mer[20]提出了通用彈簧質(zhì)點(diǎn)模型表達(dá)線纜這類柔性零件的外形，并以Cosserat理論對(duì)其彎曲和扭轉(zhuǎn)進(jìn)行建模，通過(guò)能量最小化過(guò)程求解線纜的平衡狀態(tài)，此方法避免了動(dòng)力學(xué)模型中的振動(dòng)現(xiàn)象，具有較好的穩(wěn)定性。但該方法采用比例積分控制實(shí)現(xiàn)模型的定長(zhǎng)約束，引入了約束力變量，并通過(guò)迭代的方法確定約束力，使得模型變得復(fù)雜，降低了求解效率；四元數(shù)標(biāo)準(zhǔn)化約束中，為了避免線纜外形出現(xiàn)“V”形的異常形變，需改變了約束表達(dá)式，將四元數(shù)的模變?yōu)槠涞箶?shù)，這使得雅克比矩陣變得復(fù)雜，降低了求解效率；在優(yōu)化問(wèn)題求解中級(jí)聯(lián)使用了多種優(yōu)化方法以保證算法收斂，但初值的選擇仍對(duì)求解效率和算法收斂性影響較大。

本文基于以上研究，提出采用Cosserat彈性桿理論對(duì)線纜進(jìn)行物理建模，考慮敷設(shè)仿真中線纜做低速運(yùn)動(dòng)，提出了線纜的運(yùn)動(dòng)過(guò)程仿真可通過(guò)每一時(shí)刻線纜的靜力學(xué)平衡狀態(tài)進(jìn)行表達(dá)，線纜的靜力學(xué)平衡狀態(tài)即為線纜具有最小彈性勢(shì)能狀態(tài)。可由Cosserat理論推導(dǎo)出線纜的彈性勢(shì)能，通過(guò)選取適當(dāng)?shù)牧P因子系數(shù)將約束轉(zhuǎn)化為罰函數(shù)，并通過(guò)信賴域方法求解非線性優(yōu)化問(wèn)題得到線纜具有最小能量的形態(tài)。由于Cosserat彈性桿模型考慮到了彈性桿的彎曲變形，并通過(guò)四元數(shù)表示線纜截面的扭轉(zhuǎn)變形(相對(duì)于基于 Kirchhoff彈性桿理論的線纜建模方法，考慮了更多物理特性)，因而本文提出的方法更為精確，能夠在線纜敷設(shè)過(guò)程的仿真中更加真實(shí)地表現(xiàn)出線纜的形態(tài)。

1　Cosserat彈性桿模型

線纜是細(xì)長(zhǎng)的柔性體，實(shí)際敷設(shè)過(guò)程中其應(yīng)變很小，可以認(rèn)為其本構(gòu)關(guān)系為彈性關(guān)系。因此采用Cosserat彈性桿理論進(jìn)行線纜物理建模，并對(duì)線纜做兩方面的假設(shè)：

(1) 幾何方面。在剛性截面假設(shè)基礎(chǔ)上，幾何形態(tài)由中心線的移動(dòng)和截面沿中心線的轉(zhuǎn)動(dòng)所體現(xiàn)，即將線纜作為一個(gè)軸向長(zhǎng)度遠(yuǎn)遠(yuǎn)大于半徑的空間圓柱體進(jìn)行研究，且需滿足：①中心線在變形前后均為二階以上的光滑曲線；②活動(dòng)線纜的長(zhǎng)度和曲率半徑遠(yuǎn)大于橫截面的尺度；③相鄰橫截面可繞中心線作相對(duì)扭轉(zhuǎn)，扭角為連續(xù)函數(shù)；④線纜在松弛狀態(tài)下無(wú)原始曲率和扭率。

(2) 物理特性方面。將線纜等效為勻質(zhì)、各向同性的圓截面彈性細(xì)桿，滿足以下條件：①桿為均勻各向同性，應(yīng)力和應(yīng)變滿足線性本構(gòu)關(guān)系；②橫截面為剛性平面，忽略彎曲引起的剪切變形，橫截面與中心線正交。

1.1彈性桿的空間位姿描述

彈性桿為細(xì)長(zhǎng)的可變形體，其可以呈現(xiàn)出彎曲、扭轉(zhuǎn)等復(fù)雜的形態(tài)。彈性桿截面的幾何中心連成空間曲線 C稱為彈性桿的中心線，在剛性截面假定基礎(chǔ)上，彈性桿的幾何形態(tài)由截面沿中心線的移動(dòng)和轉(zhuǎn)動(dòng)所體現(xiàn)，因此可以通過(guò)彈性桿中心線構(gòu)建出彈性桿的外形。

為了方便描述彈性桿的空間姿態(tài)，需要建立多個(gè)坐標(biāo)系。如圖1所示，對(duì)于彈性桿的中心線C，以曲線上的一固定點(diǎn) P0為原點(diǎn)建立弧坐標(biāo)系 s，以空間中一固定點(diǎn)O建立參考坐標(biāo)系(O?ξηζ) ，O點(diǎn)到空間中一點(diǎn)的向量為r。在曲線上任意一點(diǎn)P可定義一個(gè)依附于曲線的右手坐標(biāo)系(P?NBT)稱為Frenet坐標(biāo)系，如圖2所示。其中T為該點(diǎn)處切線方向上的單位矢量，N為該點(diǎn)處的法線方向上的單位矢量，矢量B由B=T×N得到。

圖1　描述彈性桿的坐標(biāo)系

圖2　Frenet坐標(biāo)系

僅通過(guò) Frenet坐標(biāo)系無(wú)法表示線纜的扭轉(zhuǎn)變形，由于彈性桿的幾何形態(tài)由截面沿中心線的移動(dòng)和轉(zhuǎn)動(dòng)所體現(xiàn)，因此在截面中心點(diǎn) P建立與剛性截面相固連的主軸坐標(biāo)系(P?xyz)，各坐標(biāo)軸的單位矢量稱為基矢量為d1、d2、d3，其中z軸與P點(diǎn)處的切線軸 T重合，即3= dT。主軸坐標(biāo)系x軸與Frenet坐標(biāo)系的N軸的夾角為θ，夾角θ為截面扭轉(zhuǎn)變形的體現(xiàn)。

1.2主軸坐標(biāo)系的空間位姿描述

主軸坐標(biāo)系的姿態(tài)通過(guò)參考坐標(biāo)系的旋轉(zhuǎn)量表示，本文采用四元數(shù)表示旋轉(zhuǎn)。

主軸坐標(biāo)系基矢量為d1、d2、d3，參考坐標(biāo)系基矢量為即：

1.3彈性桿的彎扭度

彈性桿的彎曲、扭轉(zhuǎn)變形程度用彎扭度矢量u表示，中心線上任意一點(diǎn) P處的彎扭度是該點(diǎn)處主軸坐標(biāo)系的角位移對(duì)弧坐標(biāo)s的變化率。與角速度的概念類似，且不同的是角位移是對(duì)時(shí)間 t的變化率。彎扭度等價(jià)于動(dòng)點(diǎn)P沿著中心線沿弧坐標(biāo)正向以單位速度做勻速運(yùn)動(dòng)時(shí)主軸坐標(biāo)系的角速度。

在計(jì)算中使用四元數(shù)對(duì)彎扭度矢量進(jìn)行表達(dá)，Schwab和Meijaard[21]推導(dǎo)出了四元數(shù)與旋轉(zhuǎn)的角速度矢量的關(guān)系為：

計(jì)算得到在主軸坐標(biāo)系中角速度矢量的各個(gè)分量為：

由于彎扭度矢量與角速度都是角位移的變化率，將式(6)中對(duì)時(shí)間變量t的微分改寫為對(duì)弧坐標(biāo)s的微分，可以得到彎扭度矢量u在主軸坐標(biāo)系中的表達(dá)：

1.4彈性桿的能量表達(dá)式

處于平衡狀態(tài)的曲桿只具有彈性勢(shì)能，所以無(wú)需考慮運(yùn)動(dòng)過(guò)程中的動(dòng)能以及耗散能。彈性勢(shì)能包括拉伸形變具有的彈性勢(shì)能 Vs和彎曲、扭轉(zhuǎn)形變具有的彈性勢(shì)能 Vb。文獻(xiàn)[18]給出彈性桿的Cossserat理論中彈性勢(shì)能表達(dá)，得到拉伸形變具有的彈性勢(shì)能Vs的表達(dá)式為：

其中剪切變形被忽略。ks為抗拉剛度，其與拉伸楊氏模量Es有關(guān)，可通過(guò)計(jì)算得到。

彎曲和扭轉(zhuǎn)變形所具有的彈性勢(shì)能Vb與彎扭度矢量有關(guān)，見式(9)：

其中，E為彎曲變形中的楊氏模量，G為扭轉(zhuǎn)變形時(shí)的切變模量，r為截面半徑。

對(duì)于理想桿，拉伸變形中的楊氏模量 Es等于彎曲變形時(shí)的楊氏模量E，由于實(shí)際線纜中存在多條線纜組成的線束，其截面并不均勻而且存在空隙，在實(shí)際測(cè)量中線纜拉伸時(shí)的楊氏模量 Es相比于彎曲時(shí)的楊氏模量E大。

1.5模型的約束

Cosserat桿模型是具有內(nèi)在約束的系統(tǒng)，由于中心線上任意一點(diǎn)處的切線方向與該點(diǎn)處主軸坐標(biāo)系的基矢量 d3的方向相同，所以存在方向一致約束：

對(duì)于四元數(shù)λ，只有當(dāng)其模等于1時(shí)，才能表示為純轉(zhuǎn)動(dòng)，因此存在四元數(shù)的標(biāo)準(zhǔn)化約束：

2　線纜的離散模型

在給出的Cosserat彈性桿的連續(xù)模型，包含兩個(gè)關(guān)于彈性勢(shì)能的積分式(8)和(9)以及兩個(gè)對(duì)坐標(biāo)的約束式(11)和(12)。在Cosserat彈性桿的應(yīng)用中通常采用離散模型，文獻(xiàn)[18]將彈性桿中心線離散為節(jié)點(diǎn)和線段，文獻(xiàn)[19]采用彈簧質(zhì)點(diǎn)模型模擬彈性桿。本文，采用Cosserat彈性桿理論建立線纜物理模型，并對(duì)線纜進(jìn)行離散表達(dá)，根據(jù)線纜截面的剛性假設(shè)，線纜的幾何形態(tài)由中心線決定。離散之后的線纜中心線為折線，已知各折線段中點(diǎn)的四元數(shù)和折線頂點(diǎn)的坐標(biāo)，用代數(shù)式表達(dá)折線上各點(diǎn)的四元數(shù)和坐標(biāo)，通過(guò)數(shù)值計(jì)算得到線纜的總能量。

2.1中心線的離散

如圖3所示，對(duì)于沒(méi)有分支的線纜取其中心線上的N個(gè)節(jié)點(diǎn)，順序連接這些節(jié)點(diǎn)得到N?1條中心線段。第i個(gè)節(jié)點(diǎn)的位置由參考坐標(biāo)系原點(diǎn)到節(jié)點(diǎn)的向量［1,N］表示，第j條中心線段對(duì)應(yīng)的四元數(shù)由1］表示。

圖3　線纜的離散模型

在對(duì)線纜進(jìn)行離散后，式(8)和(9)中的微分項(xiàng)可以近似求出。第i與i+1個(gè)節(jié)點(diǎn)間的空間向量ir對(duì)弧坐標(biāo)s的微分為：

其中，li是線纜處于自然狀態(tài)時(shí)，第i條中心線段的長(zhǎng)度。

第j條中心線段上四元數(shù)對(duì)弧坐標(biāo)s的微分為：

其中，lj是曲桿處于自然狀態(tài)時(shí)，第j條和第j+1條中心線段和的一半。

計(jì)算中近似給出第j條中心線段上任意一點(diǎn)處的四元數(shù)為：

2.2離散模型的能量表達(dá)式

將式(13)帶入式(8)中，對(duì)第i條中心線段在長(zhǎng)度li上進(jìn)行積分，得到拉伸變形具有的彈性勢(shì)能為：

將彎扭度用四元數(shù)表示，對(duì)第j條中心線段在長(zhǎng)度lj上積分，得到彎曲、扭轉(zhuǎn)變形具有的彈性勢(shì)能為：

將上式代入式(14)和(15)可以進(jìn)一步計(jì)算出到第j條中心線段彎曲、扭轉(zhuǎn)變形具有的彈性勢(shì)能，假設(shè)線纜在松弛狀態(tài)下無(wú)原始彎扭度，那么uk(1≤k≤3)都為0。

3　求解離散模型的平衡狀態(tài)

線纜敷設(shè)仿真中的核心是解決在線纜定長(zhǎng)的情況下，其某一點(diǎn)運(yùn)動(dòng)時(shí)對(duì)線纜其他點(diǎn)的“聯(lián)動(dòng)”位置求解(即線纜的運(yùn)動(dòng)過(guò)程求解)[22]。Spillmann 和 Teschner[18]采用動(dòng)力學(xué)方法求解彈性桿的運(yùn)動(dòng)過(guò)程，其方法存在振動(dòng)現(xiàn)象，與現(xiàn)實(shí)中的線纜運(yùn)動(dòng)有一定差別，另外動(dòng)力學(xué)方法中的時(shí)間步選擇不當(dāng)還會(huì)導(dǎo)致計(jì)算效率過(guò)低或結(jié)果不收斂。在線纜敷設(shè)過(guò)程中，對(duì)線纜的操作大多是在低速下完成的，沒(méi)有明顯的振動(dòng)現(xiàn)象發(fā)生，忽略速度造成的影響認(rèn)為其處于靜力學(xué)平衡狀態(tài)，所以線纜敷設(shè)仿真的關(guān)鍵是求解線纜處于平衡狀態(tài)時(shí)的形態(tài)。

本文通過(guò)求解線纜彈性勢(shì)能的優(yōu)化問(wèn)題，獲得其有最小彈性勢(shì)能的姿態(tài)，即線纜的平衡狀態(tài)。其中線纜的彈性勢(shì)能由Cosserat理論推導(dǎo)得出，并由罰函數(shù)將優(yōu)化問(wèn)題轉(zhuǎn)化為無(wú)約束優(yōu)化問(wèn)題，最后采用信賴域方法對(duì)該優(yōu)化問(wèn)題進(jìn)行求解。

3.1線纜彈性勢(shì)能優(yōu)化模型

通過(guò)式(16)和(17)計(jì)算出彈性桿具有的彈性勢(shì)能為：

在一定約束下，當(dāng)線纜具有最小彈性勢(shì)能時(shí)，其處于平衡狀態(tài)。本文通過(guò)求彈性勢(shì)能 V的最小值，得到彈性桿的平衡狀態(tài)。屬于約束非線性最優(yōu)化問(wèn)題，優(yōu)化模型為：

3.2優(yōu)化模型中的約束施加

模型中存在多個(gè)等式約束，在處理計(jì)算機(jī)圖形學(xué)約束的模型時(shí)，常用到罰函數(shù)方法。對(duì)于等式約束C（x）=0，將其轉(zhuǎn)化為二次能量項(xiàng)加入目標(biāo)函數(shù)中。當(dāng)解不滿足約束條件時(shí)，就增加目標(biāo)函數(shù)的值。

其中方向一致約束變?yōu)椋?/p>

四元數(shù)標(biāo)準(zhǔn)化約束變?yōu)椋?/p>

在離散模型中對(duì)式(22)和(23)進(jìn)行計(jì)算，得到每一段的罰函數(shù)能量表達(dá)式，其中方向一致約束為：

四元數(shù)標(biāo)準(zhǔn)化約束為：

引入罰函數(shù)把約束改寫能量表達(dá)式(22)和(23)，將約束最優(yōu)化式(19)轉(zhuǎn)化為無(wú)約束最優(yōu)化問(wèn)題：

文獻(xiàn)[20]指出當(dāng)抗彎剛度 kb1、kb2和抗扭剛度kb3較大時(shí)，彈性桿會(huì)發(fā)生“V”形的變形，其原因?yàn)閺澢团まD(zhuǎn)變形集中于發(fā)生急劇變形處，其他部分不發(fā)生形變，導(dǎo)致“V”形部分的四元數(shù)的模接近于零。解決方法為改變四元數(shù)標(biāo)準(zhǔn)化約束能量式為：

當(dāng)四元數(shù)的模為零時(shí)，能量將會(huì)變?yōu)闊o(wú)窮大，從而避免了以上問(wèn)題的發(fā)生。但此方法會(huì)增加求解的困難，降低求解的效率。實(shí)際上產(chǎn)生上述問(wèn)題的原因在于抗彎剛度和抗扭剛度過(guò)大使得罰函數(shù)的約束作用消失，通過(guò)增大罰函數(shù)罰因子系數(shù)同樣可以防止這種情況發(fā)生，從而避免了改寫四元數(shù)標(biāo)準(zhǔn)化約束形式引起的求解困難。

3.3能量?jī)?yōu)化模型的求解

文獻(xiàn)[20]在一次求解步驟中順序地采用了5種優(yōu)化算法，即牛頓-高斯法、非線性共軛梯度法、線性共軛梯度法、最速下降法、一維搜索法，以保證得到能量?jī)?yōu)化問(wèn)題的最優(yōu)解。因?yàn)槠湓蚴桥ｎD-高斯法只有初始值和最優(yōu)解足夠接近時(shí)，才能收斂。線纜的空間形態(tài)復(fù)雜，用以表示線纜形態(tài)的狀態(tài)量較多，如式(24)中存在較多的變量，其初值往往難以給定，如果初始點(diǎn)與最優(yōu)值差別較大，采用以上方法無(wú)法求得最優(yōu)解。為了避免初始點(diǎn)選取的困難，本文采用信賴域方法[23]，該算法為全局收斂方法，對(duì)初始點(diǎn)要求不高，使用該方法求解具有較高的穩(wěn)定性。

3.4算法的計(jì)算效率

確保線纜模型求解的實(shí)時(shí)性是實(shí)現(xiàn)虛擬環(huán)境下線纜敷設(shè)過(guò)程仿真必須滿足的要求，在確定罰函數(shù)罰因子系數(shù)及終止條件后，求解的效率主要由采用的優(yōu)化算法決定。信賴域方法是全局收斂算法，初值不影響最優(yōu)值的取得，但對(duì)求解時(shí)間有較大影響。當(dāng)初值在最優(yōu)值附近時(shí)，求解速度較快，由于操作過(guò)程相鄰位置線纜狀態(tài)變化不大，因此可將前一位置的線纜狀態(tài)作為其后線纜位姿求解的初值，能夠顯著提高求解效率。圖 4所示為17個(gè)位置求解時(shí)間，由于模型求解的第一位置，初值任意給定，求解時(shí)間較長(zhǎng)，其后各位置求解耗時(shí)降低(計(jì)算時(shí)間控制在60 ms)，滿足了實(shí)時(shí)性要求。

圖4　線纜各位置求解時(shí)間

4　線纜敷設(shè)仿真實(shí)例驗(yàn)證

采用本文提出的方法，開發(fā)了虛擬環(huán)境中的線纜敷設(shè)過(guò)程仿真軟件系統(tǒng)，系統(tǒng)采用2.80 GHz Intel Core i5-2300 CPU，顯卡為NVIDIA GeForce GT 240，內(nèi)存4 GB，操作系統(tǒng)為Window 7 SP1。

該系統(tǒng)利用三維造型引擎 ACIS和三維顯示交互工具包HOOPS建立三維環(huán)境，其中的結(jié)構(gòu)件模型是在Pro/E中建模，并通過(guò)數(shù)據(jù)轉(zhuǎn)換接口導(dǎo)入到本系統(tǒng)中。線纜的外形由物理模型計(jì)算得出，使用Matlab優(yōu)化工具箱(Optimation Toolbox)求解模型中的能量?jī)?yōu)化問(wèn)題，利用Matlab與C++的混合編程實(shí)現(xiàn)算法集成，以鍵盤和鼠標(biāo)為輸入設(shè)備，實(shí)現(xiàn)對(duì)線纜及其他虛擬物體的裝配操作。

彈簧質(zhì)點(diǎn)、能量曲線及逆運(yùn)動(dòng)學(xué)模型無(wú)法對(duì)線纜扭轉(zhuǎn)變形進(jìn)行模擬，但本文提出的物理模型能模擬這種形變，圖 5所示為線纜操作中接頭旋轉(zhuǎn)，可真實(shí)反映線纜的扭轉(zhuǎn)變形，并與真實(shí)情況相符。

圖5　線纜插裝操作

對(duì)于多根線纜在敷設(shè)過(guò)程中通常需要捆扎，操作扎帶可以帶動(dòng)線纜一起運(yùn)動(dòng)，捆扎效果如圖6所示。

線纜敷設(shè)過(guò)程仿真實(shí)例如圖 7所示，首先將所需的結(jié)構(gòu)件模型和線纜模型導(dǎo)入虛擬環(huán)境中；然后將線纜移動(dòng)到裝配位置附近，操作接頭可對(duì)接頭進(jìn)行插裝，之后通過(guò)卡箍將線纜固定到結(jié)構(gòu)件上，對(duì)于多根線纜用扎帶進(jìn)行捆扎，最終完成線纜的敷設(shè)。

圖6　線纜捆扎及扎帶調(diào)整

圖7　線纜敷設(shè)過(guò)程仿真實(shí)例

5　結(jié)　論

(1) 針對(duì)面向柔性線纜敷設(shè)過(guò)程仿真的線纜物理建模問(wèn)題，提出并建立了基于Cosserat彈性桿理論的柔性線纜物理模型，其考慮了線纜的彎曲和扭轉(zhuǎn)變形，能夠?qū)Ψ笤O(shè)過(guò)程中線纜的形變進(jìn)行真實(shí)的模擬。

(2) 通過(guò)罰函數(shù)方法對(duì)模型中的約束進(jìn)行轉(zhuǎn)化，選取適當(dāng)?shù)牧P因子系數(shù)，實(shí)現(xiàn)了模型中的四元數(shù)標(biāo)準(zhǔn)化約束和方向一致約束。

(3) 采用信賴域方法求解線纜模型中的非線性優(yōu)化問(wèn)題，該方法具有全局收斂特性，避免了初始點(diǎn)難以確定的問(wèn)題，并采用將上一位置的線纜狀態(tài)作為線纜下一位置模型計(jì)算的初始點(diǎn)的方法，提高了線纜物理模型求解的效率，滿足了虛擬環(huán)境下線纜敷設(shè)仿真的實(shí)時(shí)性要求。

(4) 開發(fā)了虛擬環(huán)境中的線纜敷設(shè)過(guò)程仿真軟件系統(tǒng)，通過(guò)實(shí)例驗(yàn)證了線纜物理模型的可行性，后續(xù)工作將考慮建立分支線纜的物理模型。

[1] 王發(fā)麟, 廖文和, 郭宇, 等. 復(fù)雜機(jī)電產(chǎn)品線纜信息本體表達(dá)與存儲(chǔ)分析[J]. 圖學(xué)學(xué)報(bào), 2015, 36(3): 376-383.

[2] 劉檢華, 萬(wàn)畢樂(lè), 寧汝新. 虛擬環(huán)境下基于離散控制點(diǎn)的線纜裝配規(guī)劃技術(shù)[J]. 機(jī)械工程學(xué)報(bào), 2006, 42(8): 125-130.

[3] 楊守勇. 三維線纜建模及布線設(shè)計(jì)的研究與實(shí)現(xiàn)[D].大連: 大連理工大學(xué), 2012.

[4] Eischen J W, Deng S, Clapp T G. Finite-element modeling and control of flexible fabric parts [J]. IEEE Computer Graphics and Applications, 1996, 16(5): 71-80.

[5] 魏發(fā)遠(yuǎn), 王峰軍, 陳新發(fā). 含有柔性電纜的復(fù)雜系統(tǒng)裝配仿真[J]. 工程設(shè)計(jì)學(xué)報(bào), 2007, 14(1): 25-31.

[6] Hergenr?ther E, D?hne P. Real-time virtual cables based on kinematic simulation [C]//WSCG ’2000, The 8th International Conference in Central Europe on Computer Graphics, Visualization and Interactive Digital Media 2000. Plzen, Czech Republic, 2000: 402-409.

[7] Terzopoulos D, Qin H. Dynamic NURBS with geometricconstraints for interactive sculpting [J]. ACM Transactions on Graphics (TOG), 1994, 13(2): 103-136.

[8] 王志斌, 劉檢華, 劉佳順, 等. 電纜虛擬布線中的物理特性分析與布局設(shè)計(jì)技術(shù)[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 2014, 26(7): 1193-1202.

[9] 尚煒, 寧汝新, 劉檢華, 等. 復(fù)雜機(jī)電產(chǎn)品中的柔性線纜裝配過(guò)程仿真技術(shù)[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 2012, 24(6): 822-831.

[10] 王志斌. 虛擬環(huán)境下基于物理屬性的線纜布局與裝配仿真技術(shù)[D]. 北京: 北京理工大學(xué), 2014.

[11] 王志斌, 劉檢華, 劉佳順, 等. 面向電纜虛擬裝配仿真的多分支彈簧質(zhì)點(diǎn)模型[J]. 機(jī)械工程學(xué)報(bào), 2014, 50(3): 174-183.

[12] Loock A, Sch?mer E, Stadtwald I. A virtual environment for interactive assembly simulation: from rigid bodies to deformable cables [C]//5th World Multiconference on Systemics, Cybernetics and Informatics (SCI’01). Orlando, Florida, USA, 2001: 325-332.

[13] 劉延柱. 彈性細(xì)桿非線性力學(xué)[M]. 北京: 清華大學(xué)出版社, 2006: 15-19.

[14] 劉檢華, 趙濤, 王春生, 等. 虛擬環(huán)境下的活動(dòng)線纜物理特性建模與運(yùn)動(dòng)仿真技術(shù)[J]. 機(jī)械工程學(xué)報(bào), 2011, 47(9): 117-124.

[15] Bertails F, Audoly B, Cani M P, et al. Super-helices for predicting the dynamics of natural hair [J]. ACM Transactions on Graphics (TOG), 2006, 25(3): 1180-1187.

[16] Bergou M, Wardetzky M, Robinson S, et al. Discrete elastic rods [J]. ACM Transactions on Graphics (TOG), 2008, 27(3): 63:1-63:12.

[17] Pai D K. Strands: interactive simulation of thin solids using cosserat models [J]. Computer Graphics Forum, 2002, 21(3): 347-352.

[18] Spillmann J, Teschner M. C o R d E: Cosserat rod elements for the dynamic simulation of one-dimensional elastic objects [C]//Proceedings of the 2007 ACM SIGGRAPH, Eurographics Symposium on Computer Animation. Eurographics Association, 2007: 63-72.

[19] 黃勁, 沈中偉, 王青, 等. 基于優(yōu)化的大步長(zhǎng)準(zhǔn)靜態(tài)線纜模擬[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 2011, 23(1): 1-10.

[20] Grégoire M, Sch?mer E. Interactive simulation of one-dimensional flexible parts [J]. Computer-Aided Design, 2007, 39(8): 694-707.

[21] Schwab A L, Meijaard J P. How to draw Euler angles and utilize Euler parameters [C]//ASME 2006 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference. American Society of Mechanical Engineers, 2006: 259-265.

[22] 劉檢華, 萬(wàn)畢樂(lè), 孫剛, 等. 線纜虛擬布線與敷設(shè)過(guò)程仿真技術(shù)[J]. 計(jì)算機(jī)集成制造系統(tǒng), 2012, 18(4): 787-795.

[23] 張光澄, 王文娟, 韓會(huì)磊, 等. 非線性優(yōu)化計(jì)算方法[M]. 北京: 高等教育出版社, 2005: 146-159.

Physical Characteristic Oriented Modeling for Flexible Cable Harness Based on Cosserat Elastic Rod Theory

Lin Haili,Liu Jianhua,Tang Chengtong,Liu Jiashun

(School of Mechanical Engineering, Beijing Institute of Technology, Beijing 100081, China)

To simulate laying process of cable harness, a physical characteristic oriented modeling for flexible cable harness is proposed based on Cosserat elastic rod theory. The cable’s total energy is obtained by applying Cosserat elastic rod theory on discrete cable model. The cable’s configuration with the lowest energy can be obtained by using optimization algorithm, which is the equilibrium state. The model takes into consideration the flexibility and continuity of the cable, which could simulate bending and twisting. From the discrete expression of the cable, the cable’s energy can be calculated by numerical method. And penalty method is used to convert the constraint problem to an unconstrained problem. The trust region algorithm is implemented to get the configuration of the cable with minimum potential energy. This method is applied in the simulation of cable laying process in a virtual environment and a prototype system is developed, and the model is validated by this system.

cable harness; physical modeling; laying simulation; Cosserat elastic rod theory; energy minimizing

TP 391.9

10.11996/JG.j.2095-302X.2016010034

A

2095-302X(2015)01-0034-09

2015-09-24；定稿日期：2015-10-15

林海立(1991–)，男，四川西昌人，碩士研究生。主要研究方向?yàn)橛?jì)算機(jī)輔助設(shè)計(jì)。E-mail：linhailisc@163.com

劉檢華(1977–)，男，江西萍鄉(xiāng)人，教授，博士，博士生導(dǎo)師。主要研究方向?yàn)閺?fù)雜產(chǎn)品裝配與檢測(cè)。E-mail：jeffliu@bit.edu.cn