撓性航天器動力學模型的非約束模態分析1)

宋新宇 戈新生

(北京信息科技大學機電工程學院，北京 100192)

引言

在撓性航天器動力學[1]的研究中，剛體模態表示撓性航天器系統整體的平移和轉動；彈性模態表示撓性附件的結構振動.彈性模態為微小振動擾動，疊加在系統整體運動之上.在撓性航天器系統的運動中，由于撓性附件結構的柔性變形所引起的慣性作用將可能導致撓性航天器運動失穩.因此，探討和研究撓性航天器系統的剛柔耦合作用問題是撓性航天器動力學研究的重點[2].在建立撓性航天器動力學模型時，可以用有限段法[3-4]、有限元法[5-6]、線性彈性理論及假設模態法[7]等方法，對太陽帆板等撓性結構的變形進行描述；撓性多體系統的模型可以用混合坐標法、運動--彈性動力學法和絕對節點坐標法[8]來描述；系統動力學方程可以用牛頓--歐拉法、拉格朗日法、凱恩方法、虛功原理等方法來建立.周志成等針對撓性組合航天器，采用假設模態法和有限元法描述柔性臂桿和太陽翼,通過動量守恒和拉格朗日方法推導了其動力學方程[9].繆炳琪等[10]采用模態綜合法建立了撓性航天器的動力學方程.劉倫[11]采用哈密頓原理建立了典型撓性航天器的動力學方程，并提出了一種直接獲取撓性航天器剛柔耦合模態的解析方法.袁秋凡等[12]針對小中心剛體--單側大撓性結構的航天器提出了一種全局模態，并基于哈密頓原理推導了其動力學方程.朱孟萍等[13-14]利用凱恩方法和假設模態法建立了不同類型的撓性航天器的動力學模型.方柳等[15]通過哈密頓變分原理和假設模態法建立了考慮動力剛化效應的撓性航天器的姿態運動和結構振動的偏微分方程.呂旺等[16]利用牛頓--歐拉向量法和模態分析方法建立了撓性衛星的動力學方程.王首喆等[17]采用虛功率原理建立了剛柔液耦合的航天器姿態軌道動力學方程.孟德山[18]針對撓性空間機器人系統，利用遞推組級法和模態綜合法建立了其動力學方程.魏進等[19]提出復合柔性結構全局模態的解析提取方法，通過全局模態離散得到系統非線性動力學模型，得到系統的固有頻率和解析函數表征的全局模態.孫家亮等[20]對近年來多柔體系統動力學建模、動力學優化及結構優化的研究成果和進展進行了綜述和總結，并提出了在多柔體系統動力學研究中值得關注的諸多問題.曹登慶等[21]從數個方面對大型撓性航天器動力學與振動控制問題，及在研究中迫切需要解決的諸多基礎科學問題作了較為全面的介紹.

利用約束模態和非約束模態，并結合混合坐標法來建立撓性航天器動力學方程是兩種主要方法.約束模態定義為中心剛體固定或固連在慣性坐標系上，撓性附件無阻尼自由振動，體現了撓性附件在固定端約束狀態下的固有特性.利用約束模態法建立的撓性系統動力學模型，接近撓性航天器在地面進行試驗測試的狀態.非約束模態定義為中心剛體不固定，整個系統(剛體、撓性附件)的無阻尼自由振動，體現了整星系統的固有特性.利用非約束模態法建立的撓性系統動力學模型，更接近衛星在軌實際飛行狀況[22-24].通常情況下，階數越高的動力學模型所反映的系統動力學特性越精確.但是，模型階數較高又不利于航天器系統控制設計和應用，因此存在動力學模型降階問題[25].慣性完備性準則[26-27]常在工程應用中被用于對系統低頻模態不太密集的動力學模型進行模態截斷和模型降階[28].所謂慣性完備性，是指對撓性結構質量或慣量的逼近程度.慣性完備性準則實質上就是忽略對撓性耦合系數(即平動撓性耦合系數和轉動撓性耦合系數)影響較小的模態，其是一個關于撓性耦合系數與撓性附件質量和慣量等特征相關的恒等式，在工程應用中有著重要的價值.Hughes 等[26,29]深入研究了模態恒等式，給出了慣性完備性準則.徐小勝等[30]曾對約束模態慣性完備性準則給出了詳細證明.

在撓性航天器動力學的研究中，一般使用約束模態法對航天器系統進行動力學建模及分析，并利用約束模態下慣性完備性準則對所建立的動力學模型進行模態截斷和模型降階.但是，隨著現代空間技術的不斷發展，航天器所承擔的任務越來越復雜，為保證任務的順利完成，航天器的規模也越來越大，結構也越來越復雜，中心剛體上通常帶有多個大型撓性附件，致使中心剛體在航天器系統中質量或慣量所占的比重越來越小，又因為要保證航天器滿足輕質要求，這樣勢必會使撓性航天器的撓性問題更突出[31].在撓性附件相對于系統質心的轉動慣量與中心剛體相對于系統質心的轉動慣量比值較大的情況下，采用傳統約束模態方法建立系統動力學模型并對其動力學問題進行分析，將會產生偏差[32].

在撓性航天器規模不斷變大，結構日益復雜的情況下，本文使用非約束模態法對航天器系統動力學問題進行分析.對于中心剛體附帶多個撓性附件的撓性航天器結構，首先利用約束模態和非約束模態方法建立了相應的動力學模型，隨后探究兩種動力學模型之間的關系，推導了非約束模態慣性完備性準則恒等式.最后，通過數值仿真計算了撓性航天器模型的非約束模態固有頻率和振型，并嘗試用非約束模態慣性完備性準則的質量關系式對撓性航天器模型進行了驗算.

1 撓性航天器系統動力學建模

圖1 所示為一個中心剛體B 帶有N個附件的撓性航天器系統結構示意圖，暫不考慮撓性附件相對中心剛體轉動的情況.建立慣性坐標系和本體坐標系分別為Oxyz和Obxbybzb，本體坐標系原點Ob為撓性附件未變形時系統的質心.Ro(t)為航天器質心相對慣性坐標系原點O的矢徑；RB(rB,t)和REn(rEn,t)分別為中心剛體上任一質點dmB與撓性附件En上任一質點dmEn到O的矢徑；rB和rEn分別為中心剛體上任一質點dmB與第n個撓性附件En上任一質點dmEn到航天器質心Ob的矢徑；撓性附件En上任一質點dmEn在變形后的位移為un(rEn,t)；航天器系統的絕對角速度為ω.假定撓性體為小變形，變形后產生的位移為一階小量，航天器的線位移、線速度也進行線性化處理.航天器的姿態用3-2-1 歐拉角描述，當歐拉角Θ(t)在小角度時可近似為=ω.

圖1 撓性航天器動力學模型Fig.1 Dynamic model of flexible spacecraft

圖1 航天器中心剛體上任一質點dmB與撓性附件上任一質點dmEn到O的矢徑可分別寫為

對式(1)和式(2)相對時間求導，可得

1.1 約束模態

為描述撓性附件的變形，引入模態坐標

式中，φjn(rEn)為第n個撓性附件En的第j階模態振型，ηjn(t)為相應的模態坐標，其滿足模態正交性條件

式中，ρ(rEn)為撓性附件質量密度，drEn為撓性附件質量微元體積，k為剛度算子，δij為符號函數，?jn為第n個撓性附件En的第j階約束模態固有頻率.

根據動量定理、角動量定理推導系統平動與轉動方程；根據第二類拉格朗日方程推導撓性附件振動方程.忽略與相關的二階小量，并注意到系統相對其質心的靜矩c=0，可得約束模態航天器系統動力學方程為

式中，mEn為撓性附件En的質量，I為單位陣，cEn和IEn分別為撓性附件En相對于Ob點的靜距和轉動慣量.令

則撓性航天器系統動力學方程(7)可寫為

1.2 非約束模態

為了描述非約束模態下系統振動情況，設撓性航天器系統中任一質點的振動位移為w(r)，則

式中，ro(t)和θ(t)為系統振動的剛性平移和轉動；r為系統中任一質點到航天器質心Ob的矢徑.引入模態坐標，可將系統中任一質點的振動位移為w(r,t)寫為

式中，Φα(r)為系統任一質點的非約束模態振型，ηα(t)為相應的非約束模態坐標

式中，roα和θα為剛體的平動和轉動模態，φαn(r)為第n個撓性附件En的第α 階非約束模態振型.其滿足非約束模態正交性條件[23]

式中，ωα為撓性航天器系統的非約束模態固有頻率，積分號中的下角標V 表示積分區域為整個航天器系統.

在非約束模態下，撓性航天器系統的平移和轉動可分別表示為

根據動量定理、角動量定理以及第二類拉格朗日方程，可導出撓性航天器系統基于非約束模態下的動力學方程為

因坐標選取在航天器系統質心，平動與轉動耦合系數分別滿足[23]

則式(16)可寫為

則式(18)可寫為

根據式(15)和式(17)，并利用式(19)中的第1 個式子，可將式(19)改寫為

2 非約束模態慣性完備性準則

對約束模態下航天器動力學方程(9)進行拉普拉斯變換，可得

對非約束模態下航天器動力學方程(20)進行拉普拉斯變換，得到

對式(28)求極限，并利用式(17)可得

因系統慣性只與系統本身有關，則式(30)與式(31)應相等，即有

對式(32)兩邊移項化簡后可得

由式(33)矩陣對應元素相等得到非約束模態下慣性完備性準則

若航天器中心剛體質心與系統質心重合，則恒等式(34)可簡化為

3 數值仿真分析

考慮兩種撓性航天器模型，即中心剛體帶雙側太陽帆板和帶單側太陽帆板組成的撓性航天器模型.中心剛體為立方體，其材料參數為[33]：楊氏模量E=2.0×1011N/m2，泊松比g=0.3，密度為7800 kg/m3，厚度為0.005 m；太陽帆板的材料參數為[33]：楊氏模量E=2.62×1011N/m2，泊松比v=0.3，密度為920 kg/m3，厚度為0.0015 m，每側太陽帆板質量為5.209 5 kg.模型整體結構及太陽帆板尺寸圖2 所示，其中太陽帆板短邊長度為0.1 m.

圖2 撓性航天器模型及太陽帆板結構尺寸Fig.2 Flexible spacecraft model and structure size of solar array

3.1 算例1

算例1 模型的中心剛體邊長尺寸為1.5 m，其質量為523 kg，中心剛體相對系統質心的轉動慣量為，撓性附件相對系統質心的轉動慣量為IE=kg· m2.不考慮太陽帆板相對中心剛體的轉動，利用Ansys 軟件對圖2 所示的模型進行非約束模態仿真分析，可計算出前6 階非約束模態固有頻率和振型，如圖3 所示.

圖3 前6 階非約束模態固有頻率和振型Fig.3 Natural frequency and mode shapes of unconstrained mode of the 1st to 6th

圖3 前6 階非約束模態固有頻率和振型(續)Fig.3 Natural frequency and mode shapes of unconstrained mode of the 1st to 6th(continued)

在中心剛體相對系統質心的轉動慣量大于撓性附件相對系統質心的轉動慣量的情況下，計算兩種模態的慣性完備性準則關于質量特性恒等式.分別求解太陽帆板前10 階約束模態平動撓性耦合系數Pjn和前10 階非約束模態平動撓性耦合系數pα，如表1 和表2 所示.

表1 前10 階約束模態平動撓性耦合系數PjnTable 1 Translational flexible coupling coefficient Pjn of constrained mode of the 1st to 10th

表2 前10 階非約束模態平動撓性耦合系數pαTable 2 Translational flexible coupling coefficient pα of unconstrained mode of the 1st to 10th

根據式(8)和式(35)以及表1 和表2 的計算數據，分別對兩種模態情況下慣性完備性準則中質量特性恒等式的理論值與前10 階數值求解，得到

從計算結果可以看出，矩陣中表征太陽帆板在x和y方向上平動的影響系數與理論值不符，其原因在于前10 階模態對于太陽帆板在x和y方向上的平移影響甚微，可近似視為無影響.

在撓性附件相對系統質心的轉動慣量小于中心剛體相對系統質心的轉動慣量時，由上述計算結果可以看出，計算結果與模型仿真實驗情況相一致，兩種模態下前10 階模態均在z方向上近似滿足關于質量特性的慣性完備性準則.

3.2 算例2

在算例1 的基礎上，改變中心剛體邊長尺寸，其他參數均不變.將中心剛體邊長尺寸改為0.5 m，其質量為57.338 kg，中心剛體相對系統質心的轉動慣量為IB=dig[3.903 2 3.903 2 3.903 2]kg·m2，撓性附件相對系統質心的轉動慣量為IE=dig[70.494 0.837 2 71.331]kg·m2.不考慮太陽帆板相對中心剛體的轉動，利用Ansys 軟件對模型進行非約束模態仿真計算，可以得到前6 階非約束模態固有頻率和振型，如圖4 所示.

圖4 前6 階非約束模態固有頻率和振型Fig.4 Natural frequency and mode shapes of unconstrained mode of the 1st to 6th

當撓性附件相對于質心的轉動慣量大于中心剛體相對于質心的轉動慣量時，因約束模態方法求解的頻率和振型將產生誤差，僅考慮非約束模態慣性完備性準則的檢測驗算.前10 階系統非約束模態平動撓性耦合系數pα如表3 所示.

表3 前10 階非約束模態平動撓性耦合系數pαTable 3 Translational flexible coupling coefficient pα of unconstrained mode of the 1st to 10th

根據式(35)和表3 的計算數據，分別可以求出非約束模態慣性完備性準則中質量恒等式的理論值與前10 階檢測數值，即

上述計算結果中表征太陽帆板在x和y方向上平動影響系數與理論值不符，其原因在于前10 階模態對于太陽帆板在x和y方向上的影響甚微，可近似視為無影響.在z方向上平動影響系數逐漸向理論值靠近，但誤差偏大.若考慮系統模態階數分別擴大到前20 階、前30 階、前40 階，分別計算可得

對比上述計算結果可以看出，當模態階數取前30 階時，表征太陽帆板在x方向上平動的影響系數變為6.0033，這是由于在第29 階模態時，太陽帆板出現了在x方向上的對稱平動.隨著所取的模態階數的增加，實驗計算結果逐漸接近理論值，可以近似滿足關于質量特性的非約束模態慣性完備性準則.

3.3 算例3

在算例2 的基礎上，只保留y軸正方向一側的帆板，中心剛體質心相對系統質心的坐標為rc=[0 ?0.1965 0]Tm，中心剛體相對其質心的轉動慣量為JB=dig[3.9032 3.9032 3.9032] kg·m2，其它參數保持不變.此時撓性附件相對系統質心的轉動慣量大于中心剛體相對系統質心的轉動慣量.不考慮太陽能帆板相對中心剛體的轉動，利用Ansys 軟件對該單側帆板模型進行非約束模態仿真分析，計算出前6 階非約束模態固有頻率和振型，如圖5 所示.

圖5 前6 階非約束模態固有頻率和振型Fig.5 Natural frequency and mode shapes of unconstrained mode of the 1st to 6th

求得前10 階系統非約束模態平動撓性耦合系數pα，如表4 所示.

根據式(34)和表4 的計算數據，分別計算非約束模態慣性完備性準則中質量恒等式的理論值與前10 階數值，得到

表4 前10 階非約束模態平動撓性耦合系數pαTable 4 Translational flexible coupling coefficient pα of unconstrained mode of the 1st to 10th

對比計算結果可以看出，前10 階計算值與理論值不符，其原因在于，當撓性附件相對系統質心的轉動慣量大于中心剛體相對系統質心的轉動慣量時，單側帆板模型的各階非約束模態對帆板在3 個方向上平動影響甚微，帆板的低階模態集中表現為彎曲和扭轉，因此采用關于質量特性的非約束模態慣性完備性準則對此類模型進行降階無實際意義，應考慮利用關于慣量特性恒等式的慣性完備性準則對其進行檢測.

4 結論

本文利用了混合坐標法建立了撓性航天器動力學模型，采用約束模態和非約束模態對航天器的振動進行展開，從而建立了撓性航天器動力學方程，探討了兩種模態情況的系統動力學方程關系.雖然約束模態和非約束模態的系統動力學方程形式有差異，但是其與系統本身相關的特征(即系統慣性)是不變的，推導出了非約束模態慣性完備性準則；隨后，利用Ansys 軟件對算例模型進行非約束模態數值仿真分析，檢測驗算了關于質量特性的非約束模態慣性完備性準則.仿真實驗結果表明：

(1)當撓性附件相對于系統質心的轉動慣量小于中心剛體時，約束模態和非約束模態的計算結果均在z方向上近似滿足關于質量關系的慣性完備性準則.

(2)當撓性附件相對于系統質心的轉動慣量大于中心剛體時，隨著所選取模態階數的增加，實驗計算結果越來越接近理論值，可近似滿足關于質量特性的非約束模態慣性完備性準則.

(3)對于撓性附件相對于系統質心的轉動慣量大于中心剛體的單側帆板撓性航天器模型，非約束模態對太陽帆板在3 個方向上的平動影響甚微，計算結果不滿足關于質量特性的非約束模態慣性完備性準則，應考慮利用關于慣量特性的非約束模態慣性完備性準則對其進行檢測驗算.

綜上所述，受限于模型各階模態對帆板在3 個方向上平動的影響，應選用慣量恒等式的慣性完備性準則進行模型檢測以及降階研究，這也是作者下一步開展研究工作的重點問題.