沈 劍,武彥偉,高懷亮

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

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

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

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

1 空間飛網動力學模型

1.1 ANCF索單元

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

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

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

(1)

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

單元坐標可以表示為：

e=[ei,ej]T

(2)

索單元被認為各向同性，因此其與縱向拉伸變形相關的應變能可以表示為[10]：

(3)

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

縱向拉伸應變εl及變形后長度ls表達式為：

(4)

與縱向拉伸變形相關的彈性力表達式為：

(5)

相應剛度矩陣表達式：

(6)

類似地，與彎曲變形相關的應變能可以表示為：

(7)

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

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

(8)

式(8)中，雙下標xx代表單元上任一點的位置矢量r對物質坐標x的二階導數。

與彎曲變形相關的彈性力表達式為：

(9)

相應剛度矩陣表達式：

(10)

單元總剛度矩陣：

K=Kl+Kt=

(11)

1.2 動力學方程及數值算法

根據第一類拉格朗日方程，建立空間飛網的運動方程為：

(12)

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

將運動方程表示為矩陣形式，有：

(13)

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

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

2 軸向應變模型和捕獲碰撞模型

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

2.1 基于弧長積分的軸向應變模型

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

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

(14)

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

(15)

式(15)中，下標x代表單元上任一點的位置矢量r對物質坐標x的導數，即單元位置梯度矢量：

(16)

(17)

(18)

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

(19)

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

圖2 文獻[7]中的L1模型

2.2 捕獲碰撞模型

在空間飛網捕獲目標時，需要進行飛網與目標的碰撞檢測。采用文獻[18]的繩網捕獲碰撞模型，對飛網上的節(jié)點與剛體進行實時碰撞檢測，可以模擬出飛網在接觸、碰撞目標后的空間構型變化過程。仿真中實時判斷繩網節(jié)點與目標中心的距離，接觸后需考慮繩索與目標的接觸力。

2.2.1碰撞檢測

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

(20)

當Δr=R時節(jié)點與球面開始接觸。

碰撞的法線方向為：

(21)

碰撞侵入量為：

(22)

2.2.2碰撞力計算

假設碰撞接觸力為F，轉化為廣義形式為：

Q=STF

(23)

將繩網與目標之間的碰撞過程可以等效為非線性彈簧阻尼模型。其中，碰撞處的恢復力由Herz接觸理論中彈簧接觸力描述，能量損失由阻尼器模擬。碰撞體受到的彈簧接觸力的方向與它們在碰撞處相互穿透的方向相反，并且承受的都是壓力，受到的阻尼力的方向為它們相對運動速度的反方向[18]。根據Hertz理論，法向接觸力的大小為：

(24)

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

3 仿真分析

3.1 材料參數

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

表1 索單元材料參數

3.2 空間飛網自由懸掛下落仿真

以網體一端固定一端自由下落為算例，通過仿真分析，繩網在不同時刻的空間構型差異如圖3所示。

圖3 不同時刻下空間飛網構型示意圖

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

圖4 繩網末端單元軸向應變曲線

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

3.3 空間飛網捕獲動力學仿真

為了驗證本文建立的空間飛網模型、軸向應變模型和捕獲碰撞模型，對空間飛網捕獲太空目標的過程進行仿真分析。仿真過程中不考慮重力作用。飛網材料參數如表1所示，其余模型參數如表2所示。

表2 空間飛網捕獲模型參數

牽引力施加于網體四個角點，沿Z向正向為正。空間飛網捕獲太空目標過程如圖6所示。

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

圖6 空間飛網捕獲太空目標過程示意圖

4 結論

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

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

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