多剛體系統(tǒng)分離策略及釋放動力學(xué)研究1)

羅操群 孫加亮 文 浩 胡海巖,? 金棟平,2)

*(南京航空航天大學(xué)機械結(jié)構(gòu)力學(xué)及控制國家重點實驗室,南京 210016)

?(北京理工大學(xué)宇航學(xué)院飛行器動力學(xué)與控制教育部重點實驗室,北京 100081)

引言

多剛體系統(tǒng)是指由多個獨立剛體通過連接約束組合而成的多體系統(tǒng),具有自由度高、約束條件復(fù)雜等特征[1-6].隨著航天技術(shù)發(fā)展及航天任務(wù)的多樣化,多剛體系統(tǒng)在到達目標(biāo)軌道后,需再次分離釋放成多個子系統(tǒng)以滿足任務(wù)需求,這些任務(wù)涉及通訊、監(jiān)控、氣象和導(dǎo)航等應(yīng)用方向[7],如衛(wèi)星互聯(lián)網(wǎng)星座任務(wù)[8]、雷達衛(wèi)星星座任務(wù)[9]等.此外,多剛體系統(tǒng)還廣泛用于航天器在軌組裝、衛(wèi)星編隊及空間碎片捕獲等領(lǐng)域[10-14].多剛體系統(tǒng)成功分離和釋放關(guān)鍵在于保證分離過程中剛體之間無碰撞發(fā)生,以使剛體按預(yù)定軌跡安全進入預(yù)定軌道.然而,在軌多剛體系統(tǒng)受科氏力及重力梯度影響,給分離過程中非線性運動及約束變化為無碰撞分離帶來困難,影響釋放安全性[15].

近年來,多剛體系統(tǒng)分離釋放動力學(xué)受到普遍關(guān)注.K?l?c? 等[16]研究了50 顆立方星的分離釋放問題,分析了不同分離頻次、分離方向和分離次序下衛(wèi)星動力學(xué)行為,發(fā)現(xiàn)所有衛(wèi)星同時釋放具有很高的碰撞風(fēng)險.Handschuh 等[17]考慮星座衛(wèi)星分離的短期碰撞風(fēng)險,對同時釋放時的分離速度和方向進行了優(yōu)化設(shè)計,通過星體間距離最大化來避免碰撞.Jeyakumar 等[18]對衛(wèi)星的分離動力學(xué)進行了分析,采用統(tǒng)計方法研究了動力學(xué)參數(shù)的變化,并且基于衛(wèi)星分離相對速度及距離等設(shè)計了分離系統(tǒng).Liu 等[19]針對多種納衛(wèi)星釋放時的安全分離問題,提出了一種基于相對軌道離心率及傾角向量的分離方法.Zhang等[20]研究了小衛(wèi)星彈射分離過程中含接觸、碰撞和摩擦的非光滑動力學(xué),提出了無干擾彈射分離的參數(shù)優(yōu)化方法.Bridges 等[21]研究了分布在火箭上級的星群間無碰撞分離釋放問題,基于釋放速度的方向,提出了兩種分離方案.Wermuth 等[22]設(shè)計了鈹衛(wèi)星從小衛(wèi)星中分離的一種分段式釋放策略,通過蒙特卡羅方法驗證了該策略的有效性.商顯揚等[23]針對一箭多星分離問題,提出了一種并聯(lián)斜置發(fā)射結(jié)構(gòu)布局及三星同時分離的方案,提高了分離可靠性.

上述研究大多關(guān)注多個剛體逐個或同時從運載航天器中分離釋放的動力學(xué)控制,僅考慮剛體與航天器之間的運動學(xué),未關(guān)注從航天器中脫離入軌后,進一步自主分離釋放的多剛體系統(tǒng)非線性動力學(xué),以及無碰撞分離問題.本文提出了一種由多個剛體緊密連接而成的在軌多剛體系統(tǒng),其分離釋放可通過剛體間的彈射機構(gòu)自主完成,以有效提高運載航天器的空間利用率,簡化分離釋放操作、降低碰撞風(fēng)險.為研究多剛體系統(tǒng)分離釋放動力學(xué),本文采用自然坐標(biāo)方法建立單剛體動力學(xué)模型,繼而連接而成為多剛體約束系統(tǒng),分析適用的連接約束條件.同時計入軌道運動、分離時剛體間相互作用,利用拉格朗日乘子法獲得多剛體系統(tǒng)分離-釋放動力學(xué)方程.基于分離方向和分離順序設(shè)計兩種分離釋放方案,規(guī)劃分離時序,實現(xiàn)了多剛體系統(tǒng)的無碰撞分離.

1 動力學(xué)建模

圖1 系統(tǒng)運動的參考坐標(biāo)系Fig.1 The reference frames of system motions

多剛體系統(tǒng)M處于近地赤道圓軌道,建立慣性坐標(biāo)系和軌道坐標(biāo)系O1XYZ和O2xyz,它們均為右手正交坐標(biāo)系,如圖1 所示.慣性坐標(biāo)系以地球中心O1為原點,O1X軸在赤道平面內(nèi)并指向春分點,O1Z軸與地球自轉(zhuǎn)軸重合;軌道坐標(biāo)系以系統(tǒng)質(zhì)心O2為原點,O2x軸和O2y軸分別指向天頂和系統(tǒng)飛行方向.

將單個剛體及系統(tǒng)均視為長方體構(gòu)型,多剛體系統(tǒng)由這N個相同的獨立剛體緊密連接而成.如圖2 所示,設(shè)多剛體系統(tǒng)共有n層,每層剛體個數(shù)為N/n.單個剛體質(zhì)量為m,質(zhì)心慣量矩陣為J,尺寸為l×b×h,其中l(wèi),b及h分別表示剛體長、寬和高.

圖2 多剛體系統(tǒng)構(gòu)型Fig.2 The configuration of MRB system

1.1 空間剛體力學(xué)模型

考慮到多剛體系統(tǒng)約束復(fù)雜,這里采用自然坐標(biāo)方法進行剛體系統(tǒng)動力學(xué)建模[24-25].自然坐標(biāo)方法較常規(guī)方法的優(yōu)勢：一是可以避免諸如歐拉角描述時的奇異性;二是大大降低施加約束的難度.在常規(guī)方法中,對于距離等簡單約束需要建立高次的約束方程且約束方程求解雅可比時需要進行多次變換.在自然坐標(biāo)方法中,僅需線性方程即可描述連接約束,并且連接約束的雅可比均為常數(shù)矩陣,利于約束分析及力學(xué)建模.需要注意的是,為便于后續(xù)分析,本節(jié)建模過程中未考慮剛體軌道運動,僅考慮剛體相對軌道坐標(biāo)系的運動.

自然坐標(biāo)方法采用剛體上兩個固定點的位矢及兩個不共面的單位矢量作為廣義坐標(biāo),以描述剛體的空間運動.兩個固定點分別為剛體質(zhì)心C和剛體上表面一點j,如圖3 所示.此外,剛體上有一固連坐標(biāo)系Cξηζ.初始時刻,坐標(biāo)軸Cξ,Cη 和Cζ 分別與軌道系的3 個主軸平行.

圖3 自然坐標(biāo)描述剛體運動Fig.3 A rigid body described by NCF method

剛體廣義坐標(biāo)在軌道坐標(biāo)系中給出,表示為

式中,rC和rj分別為固定點C和j的坐標(biāo)向量,u和v為兩個不共面單位向量.這里注意,本文將u和v選為坐標(biāo)軸Cη 和Cζ 的單位向量且固定點j位于Cξ軸上.

這樣,剛體上任意一點P在坐標(biāo)系O2xyz中的位置坐標(biāo)可以表示為

式中,I3為3×3 的單位矩陣.a1,a2和a3表示為局部坐標(biāo)系下的C和j點位置及單位向量u和v的函數(shù),即

式中,上橫線指矢量在局部坐標(biāo)系Cξηζ 下的坐標(biāo)表示.對于剛體上固定的一點P而言,A矩陣為常數(shù),故剛體上任一點的速度及加速度為

根據(jù)虛功原理[23],慣性力和外力在剛體虛位移上的虛功為零,即

式中,δWF和δWI分別為作用在剛體上的外力和慣性力所產(chǎn)生的虛功,它們是

式中,F為作用在剛體上的廣義外力,Ms為自然坐標(biāo)方法描述下剛體的質(zhì)量陣.

通過式(6),可給出剛體質(zhì)量陣為

單個空間剛體有6 個自由度,自然坐標(biāo)方法采用12 個廣義坐標(biāo),意味著這些廣義坐標(biāo)并非相互獨立,故存在6 個約束,稱之為剛體自然坐標(biāo)固有約束方程.以單個剛體為例,設(shè)廣義坐標(biāo)為q1-q12,則固有約束可以表示為

式中,L表示點C和j之間的距離.上述6 個約束方程中,前3 個代表剛體局部坐標(biāo)系主軸單位向量的數(shù)值約束,后3 個代表局部坐標(biāo)系3 個主軸互相垂直的方向約束.

將式(8)記為Φs,利用拉格朗日乘子法,未計入軌道運動的單個剛體動力學(xué)方程為

式中,Φs,q表示約束方程相對于q的雅可比矩陣,λs為拉氏乘子向量.

1.2 系統(tǒng)約束分析

根據(jù)式(8)可知,單個剛體固有約束方程為6 個,故N個剛體組裝成的多剛體系統(tǒng),其固有約束方程為6×N個.

圖4 中剛體1 和2 通過連接約束組裝成一個整體系統(tǒng),剛體灰色面代表剛體正面.連接時,兩個剛體背面相互貼合、無間隙.兩個剛體質(zhì)心分別為C1和C2,連接約束點A為兩個剛體所共有.A點坐標(biāo)可以用剛體1 和2 的質(zhì)心坐標(biāo)表示為

式中,rC1和rC2表示剛體1 和2 質(zhì)心位矢,rC1A和rC2A分別表示剛體1 和2 質(zhì)心指向約束點A的位矢.另外,剛體1 和2 的還存在方向約束

將式(10)和式(11)組成的剛體間連接約束投影到各個坐標(biāo)軸上,即得9 個約束方程.通過施加這些連接約束,剛體1 和2 組裝成一個多剛體系統(tǒng).若該剛體系統(tǒng)需要與新增加的剛體連接,則需再確定一個連接約束點,并建立與相連接剛體的局部坐標(biāo)系間的關(guān)系,即每多一個剛體連接,則多出9 個約束方程.一層中N/n個剛體間的連接約束為9×(N/n-1)個,n層共9×(N-n)個約束方程.

圖4 兩個剛體之間的連接約束Fig.4 The connecting constraints between two rigid bodies

對于層與層之間的連接約束,每層內(nèi)的剛體已經(jīng)通過施加約束連接成一個整體,故視為一個剛體.因此,層與層之間的連接只需要建立9 個約束方程即可.對于n層多剛體系統(tǒng),層間連接約束共有9×(n-1)個.

基于上述分析,由N個獨立剛體組合連接而成的n層多剛體系統(tǒng),總約束方程Φt數(shù)目為

注意的是,進行分離釋放時,需要根據(jù)分離構(gòu)型減少連接約束,從而更新約束方程.

1.3 系統(tǒng)動力學(xué)及解算方法

O

1

XYZ

式中,θ 為O2x與O1X軸之間的夾角,,ω 表示軌道角速度,如圖5 所示.

圖5 慣性系和軌道系下的廣義坐標(biāo)Fig.5 The generalized coordinates between inertia and orbital frames

根據(jù)式(13),在慣性系中,系統(tǒng)廣義速度及廣義加速度為

未計入約束時,系統(tǒng)在慣性系中的動力學(xué)方程表示為

式中,M為系統(tǒng)廣義質(zhì)量陣,FT為系統(tǒng)受到的分離彈射力,Fg表示系統(tǒng)萬有弓力[27].

將式(16)代入式(17)后,并左乘CT得到軌道系下的系統(tǒng)動力學(xué)方程

式中,左邊第一項為廣義慣性力,第二項為廣義科氏力,第三和第四項為廣義離心力,且有

基于上述分析,計入系統(tǒng)約束并利用拉格朗日乘子法,即可得到多剛體系統(tǒng)軌道動力學(xué)方程

式中,Φr表示系統(tǒng)總約束方程,λt為拉氏乘子向量,為系統(tǒng)廣義力,即

下面采用基于Newmark 方法[28-29]發(fā)展而來的廣義α 方法[30-31]來求解約束微分方程(20)的高效求解.廣義α 方法將指標(biāo)為3 的方程(20)經(jīng)過差分直接離散成代數(shù)方程進行求解,其迭代過程如下

式中,qr,n和qr,n+1為離散后的第n和第n+1 迭代步的系統(tǒng)廣義坐標(biāo)向量,并滿足如下關(guān)系

式中,h為迭代步長,矢量參數(shù)a為新弓入的算法輔助參數(shù)列陣,由下式確定

式(24)中各參數(shù)的選取方法如下[32]

對于時間步tn+1,為了迭代求解式(23),可根據(jù)下式對廣義坐標(biāo)qr,n+1、廣義速度、廣義加速度和拉氏乘子λr,n+1進行更新將算法參數(shù)

這樣,式(26)中的修正項Δqr和Δλr可由下式計算得到

式中,g表示式(20)的殘差項

而D表示g關(guān)于qr和λr的雅可比矩陣

2 分離釋放方案設(shè)計

為實現(xiàn)多剛體系統(tǒng)無碰撞分離釋放,需要基于具體構(gòu)型設(shè)計分離釋放方案.如圖6 所示,本文考慮一種2×2×4 的多剛體系統(tǒng),即每層2 行2 列4 個剛體連接,n=4 層共N=16 個剛體組裝而成的系統(tǒng).為便于分析,按如圖6 所示方式對剛體進行編號.規(guī)定：層數(shù)從上往下分別為第1～4 層.每層衛(wèi)星編號,以第一層為例,前側(cè)從左至右分別為第1～2,后側(cè)從左至右分別為第3～4 號衛(wèi)星,以下每層按照第一層順序依次疊加進行編號.

圖6 16 個剛體組成的多剛體系統(tǒng)Fig.6 The MRB system of 16 rigid bodies

基于實際工程應(yīng)用,依靠較成熟的彈射分離方式進行多剛體的分離釋放[33-34].彈射裝置對稱安裝在剛體每個表面的4 個角上,每次分離時,剛體接觸面4 個角上的彈射裝置同時開啟,產(chǎn)生大小相同且方向垂直于剛體接觸面的彈射力FT以完成分離機動.需要注意的是,每次分離時的接觸面為分離成兩部分的物體之間的接觸面.比如,第1 層與第2 層分離時,接觸面為第1 層下表面和第2 層上表面,并且此時開啟的彈射機構(gòu)位于第1 層下表面4 個角與第2 層上表面4 個角上.

采用不同彈射方向及彈射分離順序,設(shè)計兩種無碰撞分離方案,如圖7 和圖8 所示.第一種即首先利用彈射機構(gòu)同時分離第1 層剛體與第4 層剛體;隨后第1 層剛體與第4 層剛體利用相同方式逐漸分離,即按照圖中所示方向先左右分離,再前后分離的順序;中間第2 層與第3 層按照先前后分離,再左右分離,最后上下分離的順序進行機動.圖中紅色箭頭表示剛體分離釋放方向,序號表示分離次序.因此,全部步驟分為4 步進行,需進行4 次彈射分離釋放.需要注意的是,圖中所示分離方向為示意圖,各剛體姿態(tài)與初始構(gòu)型相同.隨著分離進行,分離方向會隨著剛體姿態(tài)變化,但方向始終垂直于分離面.

圖7 第一種分離釋放方案Fig.7 The first scheme of separation deployment

圖8 第二種分離釋放方案Fig.8 The second scheme of separation deployment

第二種方案即首先利用彈射機構(gòu)分離上三層剛體與第4 層剛體;隨后,第1 層與第3 層首先同時以相反方向和第2 層分離,隨后這三層再同時左右分離,最后分別前后分離成12 顆獨立剛體;第4 層先左右分離分成兩列,再前后分離為4 個獨立剛體.

上述兩種分離方案的分離方式均采取等間隔分離,彈射裝置可使彈射力FT的沖擊時間持續(xù)ΔtT.隨后進行ΔtV的自由運動,即每一次分離間隔為Δt=ΔtT+ΔtV.特別地,第一次分離開始于t1,隨后每次彈射開始于Δt·(k-1),k表示第k次彈射.

3 算例研究

為研究多剛體系統(tǒng)無碰撞分離釋放過程中的非線性動力學(xué)及驗證分離方式的有效性,本文基于兩種分離釋放方案進行動力學(xué)仿真.設(shè)多剛體系統(tǒng)位于400 km 赤道圓軌道.單個剛體質(zhì)量m=150 kg,尺寸為0.8 m×0.8 m×0.4 m.一般情況下,在軌飛行器由于硬件安裝及特殊的結(jié)構(gòu)設(shè)計,質(zhì)量分布并不均勻.因此,考慮一種質(zhì)量分布不均勻的剛體,質(zhì)心慣量矩陣為

每個彈射機構(gòu)可產(chǎn)生推力FT=1500 N,分離釋放階段的時間規(guī)劃為

為進行剛體分離碰撞分析,首先計算每個時刻每個剛體與其他15 個剛體質(zhì)心之間的間距,再將每個時刻的間距進行對比,分別取出每個時刻這些間距中的最小值與最大值

式中,dΔ,tk指的是tk時刻下第i和第j個剛體之間的距離的集合,i=1,2,···,N,j=1,2,···,N.dmin,tk和dmax,tk即tk時刻這些間距的最小值與最大值.

圖9 和圖10 分別給出了兩種分離方案中的剛體質(zhì)心之間最小間距和最大間距的時間歷程.最小間距時間歷程如圖9 所示,圖中Case 1 和Case 2 分別代表第一種和第二種分離釋放方案,下面稱為方案1 和方案2.可以看出,方案1 中最小間距在前15 s內(nèi)一直維持在0.4 m,這是因為剛體系統(tǒng)沒有完全分離,最小間距為剛體包絡(luò)最小尺寸,即上下連接的兩剛體質(zhì)心間距.而在方案2 中,第5 s 后最小間距變化為0.8 m,這是由于此時系統(tǒng)中沒有層間連接存在,最小間距變?yōu)橥粚又袆傮w連接的最小間距即剛體的長或?qū)?符合系統(tǒng)方案2 中的構(gòu)型變化.兩種方案結(jié)果相同的是在完全分離后,剛體最小間距持續(xù)增大且都大于衛(wèi)星本體包絡(luò)尺寸,表明兩種方案分離釋放過程中及分離之后剛體之間的間距均未小于剛體包絡(luò)最小尺寸,即沒有碰撞發(fā)生.此外,完全分離后,方案2 的最小間距上一直大于第一種分離方式.

圖10 給出了兩種分離方案下的剛體最大間距的時間歷程.從10 圖可見,兩種方案的最大間距在第一次分離后一直增大,并且方案2 的最大間距僅略小于方案1 的最大間距,兩種方案的最大分離范圍近乎相等.

圖9 兩種分離方案剛體最小間距時間歷程Fig.9 Time histories of the minimum distance between rigid bodies in two separation deployment schemes

圖10 兩種分離方案剛體最大間距時間歷程Fig.10 Time histories of the maximum distance between rigid bodies in two separation deployment schemes

上述結(jié)果表明,兩種分離方案的分離范圍幾乎相同,但完全分離后方案2 最小間距較大,分離效果更好.

為研究剛體分離釋放過程中的非線性動力學(xué)行為,選取5 號剛體作為分析對象,給出其分離過程中位移及歐拉角的變化,如圖11～圖14 所示.圖11 和圖12 分別給出了在軌道系O2xyz中兩種分離釋放方案下5 號剛體的位移時間歷程,顯示了剛體位移受彈射力影響的變化.可以看出,彈射后剛體逐漸遠離初始軌道.此外,第二種分離方案下,5 號剛體位移僅發(fā)生3 次變化,這是由于進行第二步分離,即在第5 s開始的彈射分離時,5 號剛體受到大小相等方向相反的彈射力,其運動等價于未受彈射力影響.

圖11 第一種分離方案5 號剛體位移時間歷程Fig.11 Time histories of displacement of the No.5 rigid body in the first separation deployment scheme

圖12 第二種分離方案5 號剛體位移時間歷程Fig.12 Time histories of displacement of the No.5 rigid body in the second separation deployment scheme

為研究兩種分離釋放方案下剛體相對于軌道系的姿態(tài)變化,同樣選取5 號剛體作為分析對象.為便于直觀理解,采用歐拉四元數(shù)λ=[λ0λ1λ2λ3]T描述姿態(tài)運動,其中

對于5 號剛體,其廣義坐標(biāo)如式(1)所示,則其體軸系3 個主軸單位向量在軌道系中的表示為

顯然,5 號剛體體軸系3 個主軸方向單位向量組成的坐標(biāo)矩陣在體軸系下的表示是單位陣I,在軌道系下表示為坐標(biāo)矩陣B=[uξuηuζ],兩者存在如下關(guān)系

因此,四元數(shù)可以由B中元素確定為

式中,Bij為B中第i行第j列元素.

圖13 和圖14 給出了兩種分離釋放方案下5 號剛體姿態(tài)四元數(shù)的時間歷程.可以看出,由于衛(wèi)星質(zhì)量分布不均勻,剛體在彈射力作用下發(fā)生連續(xù)旋轉(zhuǎn)現(xiàn)象.

圖13 第一種分離方案5 號剛體姿態(tài)四元數(shù)時間歷程Fig.13 Time histories of Euler quaternions of the No.5 rigid body in the first separation deployment scheme

圖14 第二種分離方案5 號剛體姿態(tài)四元數(shù)時間歷程Fig.14 Time histories of Euler quaternions of the No.5 rigid body in the second separation deployment scheme

最后,給出兩種釋放分離方案下系統(tǒng)在軌道系中的最終構(gòu)型,以便更加直觀地說明分離情況.如圖15 和圖16 所示,方案1 的剛體分離后呈無規(guī)則分布,而方案2 中剛體有序分布,呈對稱構(gòu)型,有利于分離完成后對剛體進行控制.

圖15 第一種分離方案系統(tǒng)最終構(gòu)型Fig.15 The final configuration of system in the first separation deployment scheme

圖16 第二種分離方案系統(tǒng)最終構(gòu)型Fig.16 The final configuration of system in the second separation deployment scheme

4 結(jié)論

本文提出了一種在軌可自主分離的多剛體系統(tǒng),研究了其分離釋放的動力學(xué)行為.采用自然坐標(biāo)方法建立了剛體動力學(xué)模型,描述其相對軌道坐標(biāo)系的運動及姿態(tài)變化,并給出了系統(tǒng)連接組合的約束條件.基于彈射分離方式,設(shè)計了兩種多剛體系統(tǒng)分離釋放方案.結(jié)果表明,對于質(zhì)量分布不均勻的剛體,兩種分離釋放方案均可實現(xiàn)多剛體系統(tǒng)的無碰撞分離,分離后多剛體系統(tǒng)的分離范圍逐漸增大,有效避免了碰撞風(fēng)險.兩種分離方案的分離范圍幾乎相同,第二種分離方案最小間距較大且分離后剛體分布規(guī)則易于分離后施加控制,分離效果更好.由于剛體質(zhì)量分布不均勻且受到分離彈射作用力的影響,分離過程中可以發(fā)生連續(xù)偏轉(zhuǎn).若分離完成后對剛體姿態(tài)有所要求,需進一步進行姿態(tài)控制.