旋轉艙體拋罩對其姿態影響的分析*

李維洲，宋曉東，王帥

(北京理工大學 宇航學院 空間運動模擬與控制實驗室，北京 100081)

0 引言

外罩分離是航天器運行過程中經常出現的問題之一，在分離過程中外罩與艙體之間有相互作用力，該作用力會導致分離后航天器本體運動狀態改變，因此需要分析外罩與艙體分離后艙體的運動狀態。

在航天器的分離問題中，由于自身結構、航天器制造工藝等因素的影響，其質心與設計位置存在一定的偏移，因此航天器分離時會受到由于質心偏移所產生力矩的干擾，從而影響其轉動狀態[1-3]。此外，物體的柔性變形也會對航天器分離過程產生一定影響[4-5]。針對航天器分離問題，盧麗穎等[6]針對衛星空間二次分離過程，用ADAMS對衛星分離后的姿態進行了定性和定量分析，研究了彈簧剛度等參數對航天器分離速度的影響；張兵[7]等利用ADAMS對三星分離的近場情況進行了計算，并詳細分析了各偏差因素對分離過程的影響；李中郢[8]等在以彈簧推桿為分離元件的前提下，建立了船箭分離過程的數學模型，通過仿真計算分析了推桿參數及分散性對飛船入軌精度、入軌安全性及平穩性的影響;舒燕[9]等基于ADAMS建立了包括航天平臺、載荷和兩類分離機構在內的仿真模型，分析了3種在軌釋放、分離載荷方案的分離過程動力學；張華[10]等建立了航天器分離動力學仿真分析模型，進行了航天器分離過程的參數敏感性分析、分離安全包絡分析，以此評估分離機構多個設計參數的偏差對航天器分離后姿態的影響；Jeyakumar D[11]建立了衛星的分離動力學模型，在此基礎上通過統計方法來計算了指定動態參數的變化。當航天器由大量部件組成時，其釋放與分離是一個具有多自由度、含時變約束的多體動力學問題，當多個航天器組成的復雜系統進行釋放時需要遵循一定次序，且要保證釋放后的航天器彼此之間不碰撞發生危險，基于這一考慮，羅操群[12]等在建立了單個剛體的動力學模型的基礎上，利用拉格朗日乘子法獲得含連接約束的非線性動力學模型，并以此為基礎設計了2種分離釋放方案，確保了航天器釋放的安全性。

為分析艙體在拋罩后的運動狀態以及造成該狀態的重要影響因素。本文基于多體動力學理論建立了旋轉艙體拋罩過程的含約束動力學模型，其中剛體的旋轉姿態由歐拉參數方法描述，時間積分算法采用可精確控制系統能量耗散特性的廣義α算法。最后，考慮彈射角、彈射速度、拋罩間隔時間等關鍵參數，進行了多種工況的計算，得到了相關的計算結果用于指導拋罩分離方案的設計和驗證。

1 旋轉艙體拋罩分離模型

1.1 旋轉艙體模型描述

在旋轉艙體拋罩分離模型中，外罩在拋罩前處于收攏狀態，該狀態下旋轉艙體的模型如圖1所示。

圖1 旋轉艙體模型Fig.1 Model of rotating cabin

圖中坐標系Oxyz為艙體坐標系，坐標系原點為艙體底面中心，Ox軸為艙體回轉軸線，Oy位于兩罩體貼合面內且與Ox垂直，Oz滿足右手坐標系。

模型由左右罩體和艙體3部分組成，在進行動力學建模僅考慮剛體間的相互作用力。艙體為一回轉體，罩體分離前與艙體固定連接，此狀態下3部分視作整體。其中，兩罩體為均質剛體，即兩罩體質心位于Oxz平面內。

如圖2所示的旋轉艙體模型中，整體模型的形心位于Ox軸上，而由于艙體質量分布不均導致整體質心偏離Ox軸一定距離。

圖2 旋轉艙體模型質心位置Fig.2 Centroid position of rotating cabin model

1.2 罩體分離方式

在該旋轉艙體拋罩分離問題中，由于氣囊充氣分離時艙體與罩體二者間的作用力十分復雜且時刻變化無法精確表達，但推著罩體逐漸充氣膨脹該力逐漸減小，因此在建模時用彈簧力近似代替合力，由于該合力不過質心，因此，分離時該合力對艙體與罩體均有力矩作用。初始時彈簧處于被壓縮狀態，艙體與兩罩體間均固連，罩體與艙體的固定副失效的同時彈簧力開始作用。

如圖3所示，彈簧兩端分別與艙體和罩體兩物體上的固定點相連，兩作用力等大反向。

圖3 彈簧作用示意Fig.3 Illustration of spring

2 旋轉艙體分離動力學模型及算法

2.1 艙體動力學模型建立

如圖4建立描述艙體分離過程的相關坐標系。

圖4 描述艙體分離過程的坐標系Fig.4 Coordinate system for describing the separation process of cabin

其中坐標系Oxyz為慣性坐標系，坐標系OAxAyAzA，OBxByBzB和O0x0y0z0分別為兩罩體和艙體的連體坐標系，其中x軸指向物體縱向慣量主軸方向，z軸與x軸垂直且位于物體的對稱平面內，y軸滿足右手坐標系。

首先定義用歐拉參數描述的物體旋轉參數為

(1)

式中：θ為轉動的軸-角表示法中的轉動角；n=[n1,n2,n3]為軸-角表示法中的單位旋轉向量。

描述剛體姿態的歐拉參數僅有3個獨立變量，根據定義有:

(2)

定義

λ=[λ1λ2λ3]T.

(3)

將歐拉參數p寫成列矢量形式為

p=[λ0λ1λ2λ3]T=[λ0λT]T.

(4)

由歐拉參數定義的方向余弦矩陣R為

(5)

據此將該式分解得到:

R=EGT,

(6)

式中：

(7)

(8)

據此，在連體基下，剛體角速度矢量可由歐拉參數速度表示為如下的運動學關系:

(9)

在慣性基下，剛體角速度矢量則表示為

(10)

連體基和慣性基下的角速度轉換關系為

(11)

艙體間依靠彈簧力模擬氣囊釋放產生的推力，以實現艙體與罩體的分離，以艙體和罩體A為例，彈簧力作用位置如圖5所示。

圖5 偏心彈簧力Fig.5 Spring force deviating from the center of mass

圖中r0，rA為慣性基下兩物體質心的位置矢量，a0,aA為質心到彈簧力作用點的位移矢量。設彈簧原長為l0，彈簧剛度為k，彈簧變形后長度矢量l定義如下：

l=(r0+a0)-(rA+aA).

(12)

則據此彈簧力表示為

(13)

式中：l0為彈簧原長；k為彈簧的剛度系數。

由于彈簧力并不過兩物體的質心，在計算作用于剛體的廣義外力和外力矩時，需要將物體所受彈簧力平移至質心處，平移過后彈簧力會對物體產生額外的力矩。單個分離體的力矢量平移示意圖如圖6所示。

圖6 力矢量平移Fig.6 Translation of force vector

圖6中，C為物體質心，a為物體質心指向力作用點的矢量，L為力矢量平移后的附加力矩，根據以上描述有:

(14)

式中：L1,L2均為3×1的列矢量。

由式(9)的變分形式可得

δθ=2Gδp.

(15)

根據虛功原理

δW=L·δθ=(2GTL)Tδp.

(16)

因此，L對應的系統所受廣義外力矩為

QL=2GTL.

(17)

經過上述變換后將力矩L由3×1列矢量轉換為一個4×1的廣義外力矩列矢量。而第i個物體平動所對應的廣義外力則無需變換，表示為

QFi=(F1i，F2i，F3i)T.

(18)

因此，第i個物體所受的廣義外載荷可寫為

Qi=[QFiQLi]T.

(19)

每個剛體的廣義外載荷矩陣均可表示為一個7×1列矢量。

該多剛體系統共含有3個剛體，對應的廣義外載荷矩陣則可寫為

Q=[Q1，Q2，Q3]T.

(20)

在多剛體系統動力學模型中，采用歐拉參數描述轉動的第i個剛體的質量矩陣可寫為

(21)

式中：

(22)

其中：mi為第i個物體的質量;mipp對應質量矩陣的轉動部分，為第i個物體4×4矩陣形式的轉動慣量；Ji為第i個物體矩陣形式的轉動慣量，

(23)

綜上所述，可得到單個無約束剛體的動力學方程為

(24)

以艙體與罩體A為例，在分離模型中，艙體分離前物體間以固定副相連接，如圖7所示。

圖7 固定連接副在不同坐標下的描述Fig.7 Description of fixed connection pair in different coordinates

在慣性基下，定義u為從慣性基原點到兩剛體固接點的矢量，即

(25)

又因為

(26)

因此有

(27)

當剛體以固定副連接時，兩剛體的固接位置矢量在全局坐標系中始終重合且描述姿態的歐拉參數始終相等。

因此，兩罩體與艙體間的固定副約束方程可寫作：

(28)

除上述方程外，每個剛體在運動時其歐拉參數在任意時刻均滿足歐拉參數的基本性質，即式(2)。

該多剛體系統由3個剛體組成，系統的總質量矩陣可表示為如下對角形式：

(29)

經過上述推導得出了該多剛體系統的動力學模型，采用拉格朗日乘子法引入相關約束方程后得到描述該多剛體系統動力學狀態的微分-代數方程組為

(30)

式中：M為統質量陣；q為系統廣義坐標；Φ為系統約束方程；λ為拉格朗日乘子；Q為系統廣義外載荷。

2.2 約束多體系統動力學方程求解算法

根據文獻[13-15]，廣義α算法在Newmark算法基礎上引入了新的算法矢量a，現采用廣義α算法對上述動力學方程進行求解。將上述描述多剛體系統動力學狀態的動力學方程在tn+1時刻進行時間離散得到：

(31)

式中：

(32)

其中，h為迭代步長；矢量a參數有如下關系：

(33)

各個參數的選取方法如下：

(34)

式中：ρ∈[0,1]為算法的譜半徑，該參數影響算法能量耗散分布的頻率范圍。針對下一個時間步tn+1，對廣義坐標、廣義速度、廣義加速度和拉格朗日乘子可按如下方式更新：

(35)

(36)

上述公式中的修正項Δq和Δλ由牛頓迭代解非線性方程組得到，可按照如下方法計算：

[ΔqΔλ]T=-J-1bn+1，

(37)

式中：b為殘差項，且其計算表達式為

(38)

J為方程關于q和λ的雅各比矩陣：

(39)

3 計算結果與分析

3.1 罩體分離前后對比

罩體分離前后，建立在艙體上的主軸坐標系原點和三主軸方向均有所變化，由于罩體質量遠遠小于艙體質量，該變化極其微小，但仍然使得原先在坐標系Oxyz的Ox方向角速度在分離后的艙體主軸坐標系Ox0y0z0中的Oy0和Oz0方向有投影，如圖8所示。

圖8 罩體分離前后的艙體連體坐標系Fig.8 Body-fixed coordinate system before and after cover’s separation

3.2 計算結果

在該計算中，彈簧在開始時即受到壓縮，固定副與彈簧同時被釋放，彈簧彈力變化趨勢如圖9所示。

圖9 彈簧力曲線Fig.9 Deflection of spring force

罩體分離過程持續時間短暫，約為1 s左右，且彈簧不會被拉伸，即彈簧從壓縮狀態伸長到原長時彈簧立刻失效，整個過程中彈簧始終保持受壓狀態。

下文中的波動幅值所表示的是波動的波峰值與中心值之差，轉速中心值即轉速平均值。

現針對不同初始條件的計算結果如下。

(1) 不同分離速度下的±z方向分離

在圖1所描述的坐標系中，±z向分離示意如圖10所示。

圖10 ±z方向分離Fig.10 Separation in ±z direction

在該分離模式下,彈簧彈力矢量位于Oxz平面內且與z軸平行，通過調節彈簧剛度使得罩體以不同的速度飛出。

圖11 罩體以不同分離速度在±z向同步 分離時艙體自旋角速度Fig.11 Spin angular velocity of the cabin when the cover is separated synchronously in ±z direction at different separation velocities

從圖11中可以看出，增大彈簧剛度時使得罩體以較大速度拋出時，由于彈簧剛度的增加，彈簧力對艙體的力矩也增大，分離后艙體在x方向角速度也越大，而y和z方向分離后變化趨勢也基本相同。

此外，對于y和z2個方向而言，當分離時刻該方向的角速度足夠小，則后續角速度會在零值附近波動，而x向角速度由于初值過大導致在圖中無法看到該波動。

(2) 罩體不同分離時間間隔下在±z方向以2 m/s速度分離

在該分離條件下，第1個罩體分離后，第2個罩體間隔2，4，6 s后再分離。在該分離模式中，當第1個罩體分離后剩余部件質心發生較大移動，慣量主軸方向變化劇烈，而當第2個罩體分離后，質心又移動到距離分離前整體質心位置的不遠處。整個分離過程中質心位置移動2次，質心移動趨勢如圖12所示。

圖12 分離過程中的質心位置移動Fig.12 Movement of the center of mass in the separation process

在該分離條件下兩罩體分離后艙體的3個方向自旋角速度如圖13所示。

圖13 罩體以相同速度在±z向異步 分離時艙體自旋角速度Fig.13 Spin angular velocity of the cabin when the cover is asynchronously separated in ±z direction at the same speed

從圖13中可看出異步分離對艙體姿態影響較大，艙體在異步分離時根據當前該方向艙體自旋角速度的不同，在第2個分離且自旋穩定后，艙體3個方向的自旋角速度會呈現較為劇烈的周期性波動。

(3) 罩體以相同速度且同步分離但彈簧力在Oxz面與Oz軸成不同角度

該分離條件下，彈簧力在圖1的坐標系中的Oxz面與Oz軸分別成15°,30°,45°,60°,75°和90°夾角，示意圖如圖14所示。

圖14 Oxz面不同角度下的分離Fig.14 Separation at different angles on Oxz plane

不同條件下分離后艙體3個方向角速度如圖15 所示。

圖15 罩體以相同速度且在Oxz面與Oz軸成 不同角度同步分離時艙體自旋角速度Fig.15 Spin angular velocity of the cabin when cover separates synchronously at the same speed at different angles to Oz axis on Oxz plane

圖15a)反應了與Oz軸夾角的大小與x方向角速度的關系，當夾角越大時，分離后對艙體的影響越小，當夾角為90°時，分離后艙體在x方向自旋角速度幾乎不變化，但對于另外2個方向卻完全相反，針對這2個方向而言，夾角越大分離后該向速度波動越大。

(4) 彈簧力在Oxz面與Oz軸成30°,罩體以2 m/s速度但在不同時間間隔分離

在該分離條件下分離時，艙體3個方向自旋角速度如圖16所示。

圖16 罩體以相同速度且在Oxz面與Oz軸成 相同角度異步分離時艙體自旋角速度Fig.16 Spin angular velocity of the cabin when cover separates asynchronously at the same speed at the same angles to Oz axis on Oxz plane

從圖16中可以發現，在這種有時間間隔的異步分離情況下，分離后x向自旋角速度波動幅值介于0.5 (°)/s(間隔4 s)至1.5 (°)/s(間隔6 s)，對比罩體沿±z軸異步分離的情況，兩者的波動情況具有很大的相似性，罩體的異步分離會造成飛行器在yz軸上較大的速度波動，其中，間隔2 s分離，yz軸的轉速波動幅值約為5 (°)/s，間隔4 s分離，yz軸的轉速波動幅值約為8 (°)/s。對比罩體沿±z軸異步分離的情況，兩者的波動幅值也基本相當。

(5) 罩體以相同速度且同步分離但彈簧力在Oyz面與Oy軸對稱成不同角度

該分離條件下，分離時彈簧力在圖1的坐標系中的Oyz面與Oy軸分別對稱成15°,30°,45°,60°,75°和90°夾角，示意圖如圖17所示。

圖17 Oyz面對稱夾角不同角度下的分離Fig.17 Separation on Oyz plane at different symmetrical angles

在該分離條件下艙體分離后的角速度如圖18所示。

圖18 罩體以相同速度且在Oyz平面與Oy對稱成 不同角度同步分離時艙體自旋角速度Fig.18 Spin angular velocity of the cabin when cover separates synchronously at the same speed at the different symmetrical angles to Oy axis on Oyz plane

通過分析圖18可以發現，罩體在Oyz平面內沿+y軸對稱打開后同步分離會造成艙體轉速下降，且打開角度越大，轉速下降越嚴重，在打開90°情況下，會使x軸平均轉速下降到68 (°)/s左右。

罩體在Oyz平面內沿+y軸對稱打開后同步分離會造成艙體轉速波動，且波動幅值隨打開角度的增加略有增加，在打開90°情況下，y軸方向最大角速度能達到-25 (°)/s，而z軸方向則是在打開角度為30°時達到最大角速度25 (°)/s。

罩體在Oyz平面內沿+y軸對稱打開后同步分離會造成艙體yz軸整體呈現類周期性章動特點，且隨著分離時張角的增大，對yz軸姿態影響整體程增大趨勢。15°張角對應章動角速度幅值最小，波動中心值和幅值均約為1 (°)/s；90°張角對應章動角速度y軸幅值約為8 (°)/s、中心值約為-17 (°)/s，z軸幅值約為8 (°)/s、中心值約為6 (°)/s。

(6) 罩體以相同速度但異步分離且彈簧力在Oyz面與Oy軸對稱成相同角度

在該分離條件下，罩體分離后艙體的3個方向自旋角速度如圖19所示。

圖19 罩體以相同速度且在Oyz平面與Oy對稱成 相同角度異步分離時艙體自旋角速度Fig.19 Spin angular velocity of the cabin when cover separates asynchronously at the same speed at the same symmetrical angles to Oy axis on Oyz plane

根據圖19分析可以發現罩體在Oyz平面內沿+y軸對稱打開30°后異步分離，對艙體的x，y，z3個軸的轉速均帶來了一定的干擾作用，且均明顯大于同步分離的干擾作用。

異步分離，對飛行器的干擾作用，不敏感于分離時間間隔。異步分離對飛行器x軸造成的速度波動幅值在0.5～2 (°)/s；對y軸造成的速度波動幅值在5.5～8.5 (°)/s；對z軸造成的速度波動幅值在5～9 (°)/s。另外異步分離對飛行器yz軸的平均轉速中心也造成了明顯的影響，其中y軸轉速平均值介于0～7 (°)/s，z軸的轉速平均值介于-1～17 (°)/s。

(7) 罩體以相同速度且同步分離但彈簧力在Oyz面與Oy軸反對稱成不同角度

該分離條件下，分離時彈簧力在圖1的坐標系中的Oyz面與Oy軸分別反對稱成15°,30°,45°,60°,75°和90°夾角，示意圖如圖20所示。

圖20 Oyz面反對稱夾角不同角度下的分離Fig.20 Separation on Oyz plane at different antisymmetrical angles

在圖20所示分離條件下，分離后艙體3個方向的自旋角速度如圖21所示。

圖21 罩體以相同速度且在Oyz平面與Oy反對稱 成不同角度同步分離時艙體自旋角速度Fig.21 Spin angular velocity of the cabin when cover separates synchronously at the same speed at the different antisymmetrical angles to Oy axis on Oyz plane

根據圖21分析可以發現，罩體的打開過程會帶來飛行器自旋速度的下降，且下降幅值與打開角度正相關，打開90°時，自旋速度下降幅值最大，為8 (°)/s。

罩體分離過程會造成艙體自旋速度的上升，且上升幅值與打開角度正相關，打開90°時，自旋速度上升幅值最大，為20 (°)/s，

罩體在Oyz平面內沿+y/-y對向打開一定角度后，以2 m/s的速度同步分離會造成飛行器yz軸產生一定的速度波動，波動幅值在0.5～1.4 (°)/s。

(8) 罩體以相同速度但異步分離且彈簧力在Oyz面與Oy軸反對稱成相同角度在該分離條件下，罩體分離后艙體3個方向自旋角速度如圖22所示。

圖22 罩體以相同速度且在Oyz平面與Oy反對稱 成相同角度異步分離時艙體自旋角速度Fig.22 Spin angular velocity of the cabin when cover separates asynchronously at the same speed at the same antisymmetrical angles to Oy axis on Oyz plane

根據圖22分析可以發現異步分離，對艙體的x,y,z3個軸的轉速均帶來了一定的干擾作用，且均明顯大于同步分離的干擾作用。

異步分離對艙體x軸造成的速度波動在0.6～1.2 (°)/s；對Y軸造成的速度波動在5～8 (°)/s；對z軸造成的速度波動在5～8 (°)/s。

異步分離對艙體yz軸的平均轉速中心也造成了明顯的影響，其中y軸轉速平均值介于-7～7 (°)/s，z軸的轉速平均值介于5～15 (°)/s。

4 結束語

通過建立旋轉艙體的動力學模型，描述了艙體分離過程分離中分離作用力方向、大小以及分離時間間隔對艙體分離這一過程的影響。

質心偏離且處于旋轉狀態的航天器在空間環境中進行分離時，根據其分離工況的不同，分離后的艙體也會有不同的運動狀態。無控狀態下的航天器分離后其3個方向角速度必然呈現等周期性的波動，波動幅值與分離時該方向的初始自旋速度有關。為了保證艙體分離后在軸線方向的旋轉穩定性，2個罩體需要被同時拋離。且為了使得分離過程二者的罩體和艙體的多次接觸碰撞，需要增大作用力使得二者分離過程加快。

通過對艙體分離過程過程的計算分析了相關因素對該過程的影響，為后續分析艙體的運動提供了參考。