空間飛網(wǎng)捕獲過程動力學(xué)研究

沈 劍,武彥偉,高懷亮

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

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

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

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

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

1.1 ANCF索單元

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

索單元屬于一維二節(jié)點(diǎn)單元，如圖1所示。在三維建模中，每個節(jié)點(diǎn)坐標(biāo)由位置矢量和該點(diǎn)的位置梯度矢量組成，共6個坐標(biāo)，每個單元含有12個坐標(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)變εl及變形后長度ls表達(dá)式為：

(4)

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

(5)

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

(6)

類似地，與彎曲變形相關(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對物質(zhì)坐標(biāo)x的二階導(dǎo)數(shù)。

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

(9)

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

(10)

單元總剛度矩陣：

K=Kl+Kt=

(11)

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

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

(12)

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

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

(13)

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

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

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

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

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

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

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

(14)

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

(15)

式(15)中，下標(biāo)x代表單元上任一點(diǎn)的位置矢量r對物質(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)時，需要進(jìn)行飛網(wǎng)與目標(biāo)的碰撞檢測。采用文獻(xiàn)[18]的繩網(wǎng)捕獲碰撞模型，對飛網(wǎng)上的節(jié)點(diǎn)與剛體進(jìn)行實(shí)時碰撞檢測，可以模擬出飛網(wǎng)在接觸、碰撞目標(biāo)后的空間構(gòu)型變化過程。仿真中實(shí)時判斷繩網(wǎng)節(jié)點(diǎn)與目標(biāo)中心的距離，接觸后需考慮繩索與目標(biāo)的接觸力。

2.2.1碰撞檢測

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

(20)

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

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

(21)

碰撞侵入量為：

(22)

2.2.2碰撞力計(jì)算

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

Q=STF

(23)

將繩網(wǎng)與目標(biāo)之間的碰撞過程可以等效為非線性彈簧阻尼模型。其中，碰撞處的恢復(fù)力由Herz接觸理論中彈簧接觸力描述，能量損失由阻尼器模擬。碰撞體受到的彈簧接觸力的方向與它們在碰撞處相互穿透的方向相反，并且承受的都是壓力，受到的阻尼力的方向?yàn)樗鼈兿鄬\(yùn)動速度的反方向[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)體一端固定一端自由下落為算例，通過仿真分析，繩網(wǎng)在不同時刻的空間構(gòu)型差異如圖3所示。

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

圖4表示繩網(wǎng)末端單元的軸向應(yīng)變變化情況。在重力作用下懸垂下落，末端單元軸向變形最大，在t=3.852 s時刻軸向應(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)捕獲動力學(xué)仿真

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

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

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

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

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

4 結(jié)論

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

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

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