沈 劍,武彥偉,高懷亮

(1.中北大學(xué) 機(jī)電工程學(xué)院， 太原 030051； 2.晉西集團(tuán)山西江陽(yáng)化工有限公司， 太原 030041；3.重慶長(zhǎng)安工業(yè)集團(tuán)有限責(zé)任公司， 重慶 401120)

根據(jù)美國(guó)宇航局(NASA)的數(shù)據(jù)，超過(guò)50萬(wàn)的碎片圍繞地球運(yùn)行，但技術(shù)上的限制使得只有大約2萬(wàn)個(gè)碎片可以被跟蹤。它們以每秒約7 800 m的速度飛行，并對(duì)人類(lèi)太空活動(dòng)構(gòu)成嚴(yán)重威脅[1]。有關(guān)于太空垃圾清除的分析研究和工程新方法有很多類(lèi)，例如美國(guó)航天局的TSS-1R任務(wù)[2]。國(guó)外關(guān)于空間飛網(wǎng)的研究主要集中在空間飛網(wǎng)捕捉以及繩索動(dòng)力學(xué)兩個(gè)方面[3]。空間飛網(wǎng)是由多個(gè)物理節(jié)點(diǎn)連接橫縱雙向多條繩索(帶)組成的，結(jié)構(gòu)輕盈、易折疊、收縮比例大[4]。在繩網(wǎng)動(dòng)力學(xué)研究方面，國(guó)內(nèi)外的研究主要還是采用集中質(zhì)量法等簡(jiǎn)化模型，如果要精確地模擬空間飛網(wǎng)的運(yùn)動(dòng)過(guò)程和構(gòu)型變化，則需要采用多體系統(tǒng)動(dòng)力學(xué)方法來(lái)進(jìn)行繩網(wǎng)柔性體建模，例如絕對(duì)節(jié)點(diǎn)坐標(biāo)法(ANCF)[5]。

不同于集中質(zhì)量模型，ANCF離散的是柔性體，并且單元的應(yīng)變模型是影響仿真精度的主要物理因素[6]。M.Berzeri和A.A.Shabana提出的縱向和橫向雙向彈性力計(jì)算模型[7]，根據(jù)是否假設(shè)縱向變形為小變形來(lái)選擇不同的計(jì)算模型。但對(duì)于單元尺寸或迭代步長(zhǎng)較大時(shí)，采用這種模型計(jì)算就會(huì)發(fā)散。Wago T等人對(duì)文獻(xiàn)[7]中的單元應(yīng)變模型進(jìn)行了改進(jìn)，建立了柔性索單元準(zhǔn)確、合適的剛度矩陣，可更真實(shí)地描述柔性體的變形過(guò)程[8]。張?jiān)降柔槍?duì)傳統(tǒng)絕對(duì)節(jié)點(diǎn)坐標(biāo)方法(ANCF)繩索模型，考慮纖維繩索初始松弛余量,推導(dǎo)了松弛繩索動(dòng)力學(xué)模型，通過(guò)仿真驗(yàn)證了松弛模型能夠更好地反映繩索在松弛狀態(tài)下的動(dòng)力學(xué)特性[9]。

本文基于ANCF法建立柔性索單元模型，針對(duì)軸向應(yīng)變模型影響計(jì)算精度的問(wèn)題，采用基于弧長(zhǎng)積分的軸向應(yīng)變模型來(lái)進(jìn)行計(jì)算。通過(guò)仿真分析研究空間飛網(wǎng)在捕獲太空目標(biāo)過(guò)程中的動(dòng)力學(xué)特性。

1 空間飛網(wǎng)動(dòng)力學(xué)模型

1.1 ANCF索單元

Gerstmayr J和Shabana A A于2006年在文獻(xiàn)[5]中提出了一種低階的索單元，該單元是基于連續(xù)介質(zhì)理論，從Euler-Bernoulli梁?jiǎn)卧葑兌鴣?lái)，通過(guò)Green應(yīng)變描述單元位移、轉(zhuǎn)動(dòng)、應(yīng)變等物理特性。對(duì)于柔性繩索，可考慮繩索軸向和彎曲應(yīng)變，忽略剪切應(yīng)變，適合于空間柔性繩網(wǎng)這種大變形、非線性的動(dòng)力學(xué)模型建模需要。單元的假設(shè)為各向同性，只受拉力作用不受壓，并且不考慮剪切、扭轉(zhuǎn)變形。

索單元屬于一維二節(jié)點(diǎn)單元，如圖1所示。在三維建模中，每個(gè)節(jié)點(diǎn)坐標(biāo)由位置矢量和該點(diǎn)的位置梯度矢量組成，共6個(gè)坐標(biāo)，每個(gè)單元含有12個(gè)坐標(biāo)。節(jié)點(diǎn)i坐標(biāo)為：

ei=[e1,e2,e3,e4,e5,e6]=

(1)

圖1 一維二節(jié)點(diǎn)索單元示意圖

單元坐標(biāo)可以表示為：

e=[ei,ej]T

(2)

索單元被認(rèn)為各向同性，因此其與縱向拉伸變形相關(guān)的應(yīng)變能可以表示為[10]：

(3)

式(3)中，εl代表縱向拉伸應(yīng)變，下標(biāo)l代表縱向Longitudinal的英文首字母；E代表單元彈性模量；A代表單元橫截面積。

縱向拉伸應(yīng)變?chǔ)舕及變形后長(zhǎng)度ls表達(dá)式為：

(4)

與縱向拉伸變形相關(guān)的彈性力表達(dá)式為：

(5)

相應(yīng)剛度矩陣表達(dá)式：

(6)

類(lèi)似地，與彎曲變形相關(guān)的應(yīng)變能可以表示為：

(7)

式(7)中，κ代表當(dāng)前構(gòu)型單元圓弧的曲率，下標(biāo)t代表橫向Transverse的英文首字母；E代表單元彈性模量；I代表單元截面慣性矩。

在建模中，假設(shè)索單元縱向變形為小變形，則曲率的表達(dá)式為：

(8)

式(8)中，雙下標(biāo)xx代表單元上任一點(diǎn)的位置矢量r對(duì)物質(zhì)坐標(biāo)x的二階導(dǎo)數(shù)。

與彎曲變形相關(guān)的彈性力表達(dá)式為：

(9)

相應(yīng)剛度矩陣表達(dá)式：

(10)

單元總剛度矩陣：

K=Kl+Kt=

(11)

1.2 動(dòng)力學(xué)方程及數(shù)值算法

根據(jù)第一類(lèi)拉格朗日方程，建立空間飛網(wǎng)的運(yùn)動(dòng)方程為：

(12)

式(12)中，C(e,t)表示約束方程。

將運(yùn)動(dòng)方程表示為矩陣形式，有：

(13)

式(13)中，λ為拉格朗日乘子；Ce為約束方程雅克比矩陣；R為約束方程對(duì)自變量的二階偏導(dǎo)數(shù)。

方程式(13)屬于微分-代數(shù)方程(DAE)。對(duì)于求解長(zhǎng)時(shí)間歷程的DAE方程，傳統(tǒng)算法不穩(wěn)定，容易產(chǎn)生數(shù)值耗散[11]。對(duì)此，Bathe[12-14]提出了一種適用于長(zhǎng)時(shí)間歷程非線性問(wèn)題的求解方法，后來(lái)Tian Q等[15]對(duì)Bathe法應(yīng)用于DAE方程求解的策略。方程式(13)系數(shù)矩陣是非對(duì)稱(chēng)的稀疏矩陣，屬于大型非對(duì)稱(chēng)稀疏矩陣問(wèn)題，可采用 GMRES 方法求解[16]。為減小內(nèi)存占用，提高調(diào)用效率，本文采用Bathe積分策略和Matlab內(nèi)置函數(shù)gmres聯(lián)合的方法進(jìn)行方程式(13)的求解。

2 軸向應(yīng)變模型和捕獲碰撞模型

在空間飛網(wǎng)動(dòng)力學(xué)仿真的過(guò)程中，計(jì)算的收斂能力與迭代步長(zhǎng)和單元結(jié)構(gòu)尺寸之間有著密切的關(guān)系。當(dāng)這兩方面因素滿足要求且計(jì)算能夠收斂后，單元的應(yīng)變曲線和節(jié)點(diǎn)的速度曲線等仍然存在高頻振動(dòng)，計(jì)算精度偏低。帶來(lái)這種計(jì)算低精度的原因主要與單元的軸向應(yīng)變模型有關(guān)[17]?？臻g飛網(wǎng)的實(shí)時(shí)構(gòu)型求解對(duì)于目標(biāo)捕獲任務(wù)的順利完成有很大影響。

2.1 基于弧長(zhǎng)積分的軸向應(yīng)變模型

文獻(xiàn)[7]中建立了ANCF的索單元軸向應(yīng)變模型L1。如圖2所示，假設(shè)單元原長(zhǎng)l，變形后縱向長(zhǎng)度為ls，模型L1計(jì)算得到縱向長(zhǎng)度為lL1。對(duì)于圖2中示例1的情況，發(fā)生彎曲變形的單元實(shí)際軸向長(zhǎng)度ls大于原長(zhǎng)l，而利用模型L1計(jì)算出來(lái)的兩節(jié)點(diǎn)間直線長(zhǎng)度lL1卻小于原長(zhǎng)，與事實(shí)不符；當(dāng)發(fā)生示例2的情況時(shí)，實(shí)際的軸向長(zhǎng)度大于原長(zhǎng)，而模型L1計(jì)算出來(lái)的數(shù)值卻等于原長(zhǎng)。當(dāng)單元變形后其軸向長(zhǎng)度大于原長(zhǎng)(即縱向應(yīng)變大于零)時(shí)，L1模型所計(jì)算出的軸向應(yīng)變不能反映真實(shí)情況，導(dǎo)致單元彈性力計(jì)算為零甚至為負(fù)情況出現(xiàn)[8]。文獻(xiàn)[8]中提出一種基于絕對(duì)節(jié)點(diǎn)坐標(biāo)公式的三維Bernoulli-Euler梁?jiǎn)卧S向彈性力的改進(jìn)形式。針對(duì)柔性梁彎曲變形較大的算例，將傳統(tǒng)的軸向彈性力計(jì)算公式(例如上述L1模型)與文獻(xiàn)[8]中的計(jì)算公式進(jìn)行了比較。對(duì)比結(jié)果表明，該公式可以用較少單元精確地表達(dá)梁?jiǎn)卧拇笞冃?。因此，針?duì)空間飛網(wǎng)運(yùn)動(dòng)過(guò)程中的大變形動(dòng)力學(xué)特性，可以采用該公式計(jì)算單元的軸向應(yīng)變。

假設(shè)變形后單元軸向長(zhǎng)度為ls，可以沿著索單元中心線對(duì)微元弧長(zhǎng)進(jìn)行積分得到：

(14)

微元的弧長(zhǎng)可以表示成單元位置梯度矢量的函數(shù),即：

(15)

式(15)中，下標(biāo)x代表單元上任一點(diǎn)的位置矢量r對(duì)物質(zhì)坐標(biāo)x的導(dǎo)數(shù)，即單元位置梯度矢量：

(16)

(17)

(18)

略去高階項(xiàng)，只取前兩項(xiàng)，并代入式(17)可得：

(19)

將式(19)代入式(4)就可以更精確地計(jì)算索單元縱向拉伸應(yīng)變，進(jìn)而計(jì)算相應(yīng)彈性力和剛度矩陣。

圖2 文獻(xiàn)[7]中的L1模型

2.2 捕獲碰撞模型

在空間飛網(wǎng)捕獲目標(biāo)時(shí)，需要進(jìn)行飛網(wǎng)與目標(biāo)的碰撞檢測(cè)。采用文獻(xiàn)[18]的繩網(wǎng)捕獲碰撞模型，對(duì)飛網(wǎng)上的節(jié)點(diǎn)與剛體進(jìn)行實(shí)時(shí)碰撞檢測(cè)，可以模擬出飛網(wǎng)在接觸、碰撞目標(biāo)后的空間構(gòu)型變化過(guò)程。仿真中實(shí)時(shí)判斷繩網(wǎng)節(jié)點(diǎn)與目標(biāo)中心的距離，接觸后需考慮繩索與目標(biāo)的接觸力。

2.2.1碰撞檢測(cè)

將目標(biāo)構(gòu)型簡(jiǎn)化為半徑為R的球體。首先通過(guò)判斷檢測(cè)點(diǎn)和剛體的空間位置對(duì)兩者進(jìn)行檢測(cè)，如果檢測(cè)到運(yùn)動(dòng)體相互之間存在侵入，則計(jì)算目標(biāo)剛體對(duì)繩索的侵入量，進(jìn)而根據(jù)繩索與目標(biāo)間的接觸模型計(jì)算接觸碰撞力，并更新運(yùn)動(dòng)體的狀態(tài)。 假設(shè)球體中心坐標(biāo)為(x0,y0,z0)，繩網(wǎng)節(jié)點(diǎn)(xi,yi,zi)是否與目標(biāo)球體接觸的判據(jù)為：

(20)

當(dāng)Δr=R時(shí)節(jié)點(diǎn)與球面開(kāi)始接觸。

碰撞的法線方向?yàn)椋?/p>

(21)

碰撞侵入量為：

(22)

2.2.2碰撞力計(jì)算

假設(shè)碰撞接觸力為F，轉(zhuǎn)化為廣義形式為：

Q=STF

(23)

將繩網(wǎng)與目標(biāo)之間的碰撞過(guò)程可以等效為非線性彈簧阻尼模型。其中，碰撞處的恢復(fù)力由Herz接觸理論中彈簧接觸力描述，能量損失由阻尼器模擬。碰撞體受到的彈簧接觸力的方向與它們?cè)谂鲎蔡幭嗷ゴ┩傅姆较蛳喾?，并且承受的都是壓力，受到的阻尼力的方向?yàn)樗鼈兿鄬?duì)運(yùn)動(dòng)速度的反方向[18]。根據(jù)Hertz理論，法向接觸力的大小為：

(24)

式(24)中,K等效接觸剛度;C等效接觸阻尼。

3 仿真分析

3.1 材料參數(shù)

空間飛網(wǎng)大多采用柔性材料編織而成，例如高強(qiáng)度尼龍織帶等。仿真中所用到的材料參數(shù)如表1所示。

表1 索單元材料參數(shù)

3.2 空間飛網(wǎng)自由懸掛下落仿真

以網(wǎng)體一端固定一端自由下落為算例，通過(guò)仿真分析，繩網(wǎng)在不同時(shí)刻的空間構(gòu)型差異如圖3所示。

圖3 不同時(shí)刻下空間飛網(wǎng)構(gòu)型示意圖

圖4表示繩網(wǎng)末端單元的軸向應(yīng)變變化情況。在重力作用下懸垂下落，末端單元軸向變形最大，在t=3.852 s時(shí)刻軸向應(yīng)變達(dá)到了0.115 7。圖5所示為繩網(wǎng)末端單元的節(jié)點(diǎn)Y向位移曲線，可知呈周期性變化。隨著能量的衰減，最終趨于0附近。

圖4 繩網(wǎng)末端單元軸向應(yīng)變曲線

圖5 繩網(wǎng)末端單元節(jié)點(diǎn)Y向位移曲線

3.3 空間飛網(wǎng)捕獲動(dòng)力學(xué)仿真

為了驗(yàn)證本文建立的空間飛網(wǎng)模型、軸向應(yīng)變模型和捕獲碰撞模型，對(duì)空間飛網(wǎng)捕獲太空目標(biāo)的過(guò)程進(jìn)行仿真分析。仿真過(guò)程中不考慮重力作用。飛網(wǎng)材料參數(shù)如表1所示，其余模型參數(shù)如表2所示。

表2 空間飛網(wǎng)捕獲模型參數(shù)

牽引力施加于網(wǎng)體四個(gè)角點(diǎn)，沿Z向正向?yàn)檎??？臻g飛網(wǎng)捕獲太空目標(biāo)過(guò)程如圖6所示。

由圖6可以看出，t=0.2 s時(shí)網(wǎng)體接觸目標(biāo)后開(kāi)始變形，并開(kāi)始包覆目標(biāo)；t=0.3 s時(shí)，目標(biāo)體一半的表面被網(wǎng)體覆蓋；t=0.4 s時(shí)，網(wǎng)體覆蓋目標(biāo)面積超過(guò)4/5，可認(rèn)為捕獲成功。仿真中網(wǎng)體變形穩(wěn)定，流線型效果較好，在一定程度上能夠反映真實(shí)的空間飛網(wǎng)捕獲過(guò)程中的姿態(tài)和構(gòu)型變化。

圖6 空間飛網(wǎng)捕獲太空目標(biāo)過(guò)程示意圖

4 結(jié)論

1) 本文采用ANCF索單元建立空間柔性繩網(wǎng)模型，引入基于弧長(zhǎng)積分的軸向應(yīng)變模型精確計(jì)算單元軸向應(yīng)變，結(jié)合繩網(wǎng)捕獲碰撞模型，可以較好地模擬空間飛網(wǎng)捕獲太空目標(biāo)的過(guò)程。

2) 本文研究成果可用于工程中各類(lèi)大變形柔性繩網(wǎng)的建模和仿真工作。

3) 計(jì)算發(fā)現(xiàn)，算法的高效性對(duì)于求解大規(guī)模網(wǎng)體的動(dòng)力學(xué)方程組十分重要。未來(lái)需要對(duì)算法的收斂性、海量數(shù)據(jù)存儲(chǔ)等問(wèn)題進(jìn)行深入研究。