張 瑾，翟 坤，王天舒

（清華大學 航天航空學院，北京 100084）

0 引言

太陽帆航天器是一種新型的深空探測航天器，依靠太陽光壓力產(chǎn)生推力。太陽帆航天器自身無需攜帶大量的燃料用于任務，但為在微弱的光壓力下獲得最大推力，太陽帆有巨大的帆面及很輕的結(jié)構(gòu)質(zhì)量，這決定了太陽帆具大轉(zhuǎn)動慣量及大撓性結(jié)構(gòu)［1］。一般大撓性結(jié)構(gòu)有剛度低、阻尼弱、頻率低和模態(tài)密集等動力學特性，即具強非線性，這為其精確動力學模型的建立帶來了困難［2］。同時在大撓性結(jié)構(gòu)產(chǎn)生運動時，結(jié)構(gòu)的振動對整體結(jié)構(gòu)的動力學特性與控制的影響不可忽略［3］。

因太陽帆利用太陽光壓力產(chǎn)生推進力，為大尺度、大撓性結(jié)構(gòu)，故太陽帆是典型的姿態(tài)、軌道、振動耦合問題。目前對太陽帆的研究主要集中于太陽帆的軌道設(shè)計、姿態(tài)控制和精確受力模型三方面。在太陽帆傳統(tǒng)的軌道設(shè)計及姿態(tài)控制中，常忽略了太陽帆的實際結(jié)構(gòu)，即忽略太陽帆撓性結(jié)構(gòu)的變形及振動，將太陽帆作為質(zhì)點或剛體進行受力分析，并通過對姿態(tài)角的改變而調(diào)節(jié)推力的大小，從而實現(xiàn)對軌道的控制［4－5］。對太陽帆及類似地有撓性附件航天器來說，其彈性振動的處理方法有兩種：一是將航天器結(jié)構(gòu)的振動簡化為影響力矩引入軌道姿態(tài)方程，二是用支撐桿的模態(tài)替代整體模態(tài)進行姿態(tài)振動耦合分析［6－9］。這兩種方法都可得到與剛體模型不同的結(jié)果，說明彈性振動對太陽帆的姿態(tài)會有影響。

太陽帆在飛行過程中會不斷產(chǎn)生變形，因帆膜剛度低，其變形較復雜，所受光壓也由此變得復雜，故建立較精確的光壓力模型也成為研究的目標之一［10］。通過建立精確的受力模型，可根據(jù)太陽帆實時的彈性變形得到實時受力狀況［11－12］。但此種方法忽略了彈性振動對軌道、姿態(tài)的影響。文獻［13］采用混合坐標法建立了柔性太陽帆動力學模型，用預應力模態(tài)完成模型降階，建模中考慮了帆膜和支撐梁的中性面耦合變形，并通過ABAQUS／Explicit的瞬態(tài)動力學分析功能驗證了該降階模型的有效性。

對太陽帆的動力學進行了大量研究，但還沒有結(jié)合了軌道、姿態(tài)、彈性振動的動力學研究。由于三者間的強烈耦合，太陽帆軌道、姿態(tài)、彈性振動的聯(lián)合控制需作深入研究。本文根據(jù)文獻［13］提出的柔性太陽帆動力學降階模型，考慮天體引力及太陽光壓力的作用，對柔性太陽帆在空間運行時軌道、姿態(tài)以及彈性振動的耦合效應進行了研究。

1 柔性太陽帆動力學建模

目前，太陽帆的研究多針對方形太陽帆構(gòu)型，一般由支撐桿4根和三角形帆膜4塊組成，在太陽帆中心可搭載小衛(wèi)星等有效載荷。方形太陽帆的結(jié)構(gòu)較簡單，當有效帆膜面積相同時，所需支撐桿的長度更少，因此結(jié)構(gòu)質(zhì)量較輕，且方形太陽帆的姿態(tài)操控更簡易和精確，能在行星逃逸時提供很大的轉(zhuǎn)彎速度。

1.1 坐標系定義

定義地球赤道慣性坐標系O-XYZ、軌道坐標系o-xyz和太陽帆本體坐標系ob－xbybzb如圖1所示。圖中：Ω為升交點赤經(jīng)；ω為近地點幅角；f為真近點角；i為軌道傾角。具體如下。

圖1 太陽帆空間坐標系Fig.1 Coordinate systems definition for solar sail

a）O-XYZ系：坐標原點為地球質(zhì)心O；OX軸指向春分點；OZ軸垂直于赤道平面；OY軸由右手法則得到。

b）o-xyz系：坐標原點為太陽帆中心在軌位置o；ox軸由地心指向航天器；oy軸垂直于軌道平面，oz軸由右手法則得到。

c）ob－xbybzb系，坐標原點為太陽帆幾何中心ob；obxb，obyb軸在太陽帆平面內(nèi)；obzb軸沿太陽帆法向。

初始時刻，ob－xbybzb，o-xyz系重合；點ob，o重合；太陽帆帆面指向地心。

O-XYZ，o-xyz系間的轉(zhuǎn)換陣Ai2o取決于軌道所處的位置矢量r和軌道速度矢量v，有

O-XYZ系到ob－xbybzb系的轉(zhuǎn)換用四元數(shù)表示為

則O-XYZ系到ob－xbybzb系的轉(zhuǎn)換矩陣為

式中：

o-xyz系到ob－xbybzb系的轉(zhuǎn)換矩陣可表示為

通過轉(zhuǎn)換計算可得太陽帆本體系相對軌道坐標系的歐拉轉(zhuǎn)角。

1.2 柔性太陽帆動力學降階模型及方程

本文選擇五點連結(jié)方形太陽帆為研究對象，如圖2所示。

圖2 五點連結(jié)方形太陽帆Fig.2 Five-point connected square solar sail

太陽帆的主要柔性結(jié)構(gòu)包括支撐梁和帆膜。假設(shè)支撐梁為歐拉－伯努利梁，忽略其剪切變形，帆膜為薄板，跟隨支撐梁產(chǎn)生變形。支撐梁與帆膜通過拉索在太陽帆的中心及4個端點相互連接，作為整體建模。太陽帆上任一點相對慣性系O點的矢徑為

式中：RO為點ob相對O-XYZ系原點O的矢徑；rf為未變形時該點相對點ob的矢徑；uf為該點的變形矢量。僅考慮柔性結(jié)構(gòu)的橫向變形及由此引起的中性面耦合變形，uf在ob－xbybzb系中的表達式為

式中：w（x，y，t）為中性面的橫向變形；（x0，y0）為中性面上不發(fā)生變形的點，即為點ob。

對橫向變形w進行模態(tài)離散可得模態(tài)坐標表達式

式中：Φ為通過預應力模型分析得到的模態(tài)振型；q為模態(tài)坐標。數(shù)值仿真時僅保留前3階模態(tài)完成模型降階。用ABAQUS／Standard預應力分析，所得模型的前三階非剛體頻率及振型如圖3所示。

圖3 柔性太陽帆前三階非剛體頻率及振型Fig.3 Frst three non-rigid-body frequencies and mode shapes of solar sail

忽略各種天體攝動力和耗散阻力，僅考慮天體引力和太陽光壓力，根據(jù)虛功率原理，有

式中：ps為作用在帆面上的太陽光壓；pg為天體引力；Kf為柔性結(jié)構(gòu)模態(tài)剛度矩陣；，分別為太陽帆上任一點的速度變分和模態(tài)速度變分；dm，dA分別為該點的質(zhì)量與面積。

將R對時間的一、二次導代入式（8）可得柔性太陽帆軌道－姿態(tài)－振動耦合方程

式中：mtot為太陽帆的總質(zhì)量；Sf為太陽帆與彈性振動耦合后的靜矩；Jf為太陽帆與彈性振動耦合后的轉(zhuǎn)動慣量；ω為本體系相對慣性系的角速度矢量；Ct為平動與彈性振動耦合矩陣；Cr為轉(zhuǎn)動與彈性振動耦合矩陣；ΦNL為變形矢量uf相對模態(tài)坐標q求導的導數(shù)矩陣；Qft，Qfr，Qff分別為平動、轉(zhuǎn)動、振動耦合產(chǎn)生的力及力矩；Ft，F(xiàn)r，F(xiàn)f分別為太陽光壓力產(chǎn)生的力及力矩；Fgt，F(xiàn)gr，F(xiàn)gf分別為天體引力產(chǎn)生的力及力矩；（·），（··）分別為慣性系中對時間的一次導和二次導，（′），（″）分別為為本體坐標系中對時間的一次導和二次導。

用四元數(shù)Λ＝［Λ0Λ1Λ2Λ3］T描述太陽帆的姿態(tài)，則式（9）在本體系中可寫成

根據(jù)四元數(shù)導數(shù)與運動角速度間的關(guān)系

可將式（10）轉(zhuǎn)換成狀態(tài)空間方程

綜上，對式（16）可直接進行數(shù)值求解，本文采用MATLAB中的龍格庫塔四－五階積分函數(shù)ode45進行數(shù)值求解。計算結(jié)果可根據(jù)各坐標系間的轉(zhuǎn)換矩陣進行換算。

1.3 太陽帆模型參數(shù)

如圖2所示，太陽帆包含支撐梁4根和等腰直角三角形帆膜4個，每個帆膜通過繩索3根與支撐梁連接，太陽帆法線方向有控制桿。在太陽帆支撐梁的4個端點有質(zhì)量0.58kg的姿控小帆，在太陽帆中心有質(zhì)量54.84kg的盒艙，在控制桿末端有質(zhì)量228kg的有效負載。太陽帆的結(jié)構(gòu)參數(shù)見表1。

因控制桿末端存在較大質(zhì)量的有效負載，故太陽帆質(zhì)心與中心（本體系原點）不重合，太陽帆質(zhì)心在本體系中可表示為

太陽帆受到的引力梯度力矩

式中：μ為地球引力常數(shù)?？煽闯觯河捎谔柗钠馁|(zhì)量，太陽帆所受的引力梯度力矩有所增大，其轉(zhuǎn)動角位移也會明顯增大。同時，由式（9）也可發(fā)現(xiàn)：由于軌道與姿態(tài)的耦合性，由光壓力引起的軌道偏差亦會對太陽帆的姿態(tài)產(chǎn)生影響。

2 軌道－姿態(tài)－彈性振動耦合效應

文獻［14］研究了剛體帆在大偏心率地球橢圓軌道上太陽帆帆面對地指向及對日指向的姿態(tài)控制過程。由于地心軌道的引力較大，并且大偏心率橢圓軌道的引力及重力梯度力矩變化較大，更有利于檢驗太陽帆軌道－姿態(tài)－振動耦合效應，因此本文對柔性太陽帆在大偏心率橢圓軌道進行動力學分析，檢驗其軌道－姿態(tài)－振動耦合效應。

不考慮對太陽帆的控制，令太陽帆只依靠天體引力運行，地心橢圓軌道相應軌道要素為a＝46 428km，e＝0.819 6，i＝12°，f0＝0°；在地球軌道上的太陽光壓力大小為ps＝9.12×10－6Pa。假設(shè)太陽帆初始時刻本體坐標系與軌道坐標系重合，即帆面指向地心，太陽帆相對地球慣性坐標系初始角速度為0。

由于太陽帆耦合動力學方程中是對其軌道、姿態(tài)及模態(tài)坐標同時進行求解，而三者的數(shù)值量級存在很大差異，為保證姿態(tài)與模態(tài)坐標的求解精度，在方程求解時積分的絕對誤差與相對誤差分別取為10－12，10－10。MATLAB中ode45積分函數(shù)是變步長積分函數(shù)，為節(jié)約輸出內(nèi)存并完整體現(xiàn)太陽帆的振動特性，根據(jù)太陽帆的自然頻率可選定輸出步長為3s。

表1 太陽帆結(jié)構(gòu)參數(shù)Tab.1 Structural parameters of solar sail

2.1 剛性太陽帆動力學

先不考慮太陽帆的柔性，則式（9）可退化為軌道、姿態(tài)耦合方程

在地心橢圓軌道進行2個軌道周期的數(shù)值仿真。仿真結(jié)果如圖4～12所示。

太陽帆運行軌道在地球赤道慣性系XOY平面上的投影如圖4所示。由圖可知：在太陽光壓力的作用下，太陽帆的軌道出現(xiàn)了偏差，軌道半長軸最大偏移量369.206km。

圖4 太陽帆XY平面內(nèi)軌道Fig.4 Representation of orbit inXYplane for solar sail

對一般航天器，在無光壓力及其他攝動力作用時，航天器僅受軌道平面法線方向上的重力梯度力矩作用，將只存在軌道平面法線方向上的角速度和角位移。由于太陽帆在大偏心率的橢圓軌道上運動，引力梯度力矩隨軌道半徑變化出現(xiàn)較大變化，軌道角速度也發(fā)生明顯變化，且由于太陽帆質(zhì)心與中心不重合，會存在較大引力梯度力矩，太陽帆相對軌道坐標系y向會出現(xiàn)較大的角位移。同時由于太陽帆動力學方程的耦合性，在太陽光壓力作用下，太陽帆相對軌道坐標系x、z向也出現(xiàn)了較大角位移，如圖5～7所示。由圖可知：太陽帆相對x向歐拉角在±50°間振蕩。相對y、z向也存在大范圍的變化，由圖8～10可知：太陽帆在軌道坐標系三個方向的角速度的量級相同，且呈持續(xù)振蕩，其最大振幅均為0.1（°）／s。

由上述剛性太陽帆在大偏心率地心橢圓軌道上的動力學分析可知：太陽光壓力引起了太陽帆相對原軌道的偏差。由于太陽帆質(zhì)心與中心的不重合，在太陽光壓力的作用下，太陽帆會在3個方向上均出現(xiàn)大幅的角位移變化。

圖5 太陽帆相對軌道坐標系x軸角位移Fig.5 Representation of angular displacement atx-axis for solar sail

圖6 太陽帆相對軌道坐標系y軸角位移Fig.6 Representation of angular displacement aty-axis for solar sail

圖7 太陽帆相對軌道坐標系z軸角位移Fig.7 Representation of angular displacement atz-axis for solar sail

2.2 柔性太陽帆動力學

考慮太陽帆的柔性，動力學方程為式（9），對柔性太陽帆在地心橢圓軌道上進行2個軌道周期的動力學仿真。

太陽帆軌道在地球赤道慣性系XOY平面的投影如圖4所示。在太陽光壓力以及自身彈性振動的作用下，太陽帆運行軌道產(chǎn)生了偏差，其軌道半長軸的最大偏移量369.216km，與剛性太陽帆的結(jié)果近似。可見，太陽光壓力是太陽帆軌道產(chǎn)生偏差的主要原因，太陽帆柔性結(jié)構(gòu)振動對其軌道影響很小。

圖8 太陽帆相對軌道坐標系x軸角速度Fig.8 Representation of angular velocity atx-axis for solar sail

圖9 太陽帆相對軌道坐標系y軸角速度Fig.9 Representation of angular velocity aty-axis for solar sail

圖10 太陽帆相對軌道坐標系z軸角速度Fig.10 Representation of angular velocity atz-axis for solar sail

與剛性太陽帆相同，柔性太陽帆相對軌道坐標系x、y、z向均出現(xiàn)了較大的角位移變化，太陽帆相對x向歐拉角在±70°間振蕩，相對y、z向歐拉角變化范圍更大。與剛性太陽帆相比，柔性太陽帆雖在三個方向上同樣存在大幅轉(zhuǎn)角，但在結(jié)構(gòu)振動的作用下，其角位移隨著時間的增長逐漸出現(xiàn)差異，并隨著時間的增長偏差逐漸增大，說明太陽帆的柔性對其姿態(tài)的影響具有時間的累積效應，且這種影響不可忽視。

由太陽帆相對軌道坐標系三個方向的角速度可知：與剛性太陽帆相似，柔性太陽帆三個方向的角速度也有相同量級，其最大振幅均為0.1（°）／s。同樣的，柔性太陽帆與剛性太陽帆的角速度變化曲線類似，但是在太陽帆柔性的影響下，太陽帆角速度曲線逐漸出現(xiàn)偏差，并隨著時間的增長偏差逐漸增大。

柔性太陽帆端點1處的橫向變形及伸長變形分別如圖11～12所示。太陽帆端點1處最大橫向變形0.65m，最大伸長變形0.002 5m。由局部放大圖可知：在端點變形的增大與減小過程中，太陽帆的柔性結(jié)構(gòu)仍存在小幅高頻的彈性振動，其振動周期約38s。

圖11 太陽帆端點1處橫向變形Fig.11 Deflection deformation at tip 1for solar sail

由柔性太陽帆在大偏心率地心橢圓軌道上的動力學分析可知：太陽光壓力引起了太陽帆軌道的偏差，以及柔性結(jié)構(gòu)的彈性變形及振動。由于大偏心率橢圓軌道上引力梯度力矩的較大變化，以及太陽帆質(zhì)心與中心的不重合，在太陽光壓力的作用下太陽帆會產(chǎn)生較大幅度的轉(zhuǎn)動角位移和角速度，進而引起太陽帆柔性結(jié)構(gòu)較大幅度的彈性變形及振動。同時，太陽帆柔性結(jié)構(gòu)的振動也會導致其姿態(tài)產(chǎn)生與剛性太陽帆不同的姿態(tài)變化，而這種影響是一個長時間的積累效應，且其影響不可忽視。

圖12 太陽帆端點1處伸長變形Fig.12 Stretching deformation at tip 1for solar sail

3 結(jié)束語

本文建立了柔性太陽帆軌道－姿態(tài)－彈性振動耦合的動力學降階模型，研究了太陽帆軌道、姿態(tài)和結(jié)構(gòu)彈性振動的耦合效應。研究表明：太陽光壓力是引起太陽帆軌道產(chǎn)生偏差的主要原因，同時會引起太陽帆柔性結(jié)構(gòu)的彈性變形及振動；當太陽帆質(zhì)心與中心不重合時，較大的引力梯度力矩以及軌道和姿態(tài)的耦合特性會引起太陽帆姿態(tài)的大幅變化；太陽帆柔性結(jié)構(gòu)的振動會隨時間的增加引起姿態(tài)的偏差，這種隨時間積累的偏差在長期的深空探測任務中非常重要。同時，通過分析也可發(fā)現(xiàn)：由于在地球軌道上太陽帆所受的引力較大，會對太陽帆姿態(tài)和彈性變形產(chǎn)生較大的影響，因此，為更好地發(fā)揮太陽帆的推力優(yōu)勢，應將太陽帆用于深空探測任務中。為克服太陽帆偏心對其姿態(tài)的影響，可通過適當?shù)淖藨B(tài)控制方案保持其姿態(tài)的穩(wěn)定，進而實現(xiàn)軌道的穩(wěn)定。柔性太陽帆在空間飛行是一個軌道、姿態(tài)、彈性振動強烈耦合的復雜問題，其耦合效應的研究對其長期深空探測任務的完成具有重要意義。

［1］ 王偉志.太陽帆技術(shù)綜述［J］.航天返回與遙感，2007（2）：1-4.

［2］ 司洪偉，李東旭，陳衛(wèi)東.大撓性航天桁架結(jié)構(gòu)動力學及其主動控制研究進展［J］.力學進展，2008，38（2）：167-176.

［3］ 劉 軍，韓 潮.撓性航天器大角度機動的振動抑制控制［J］.系統(tǒng)仿真學報，2008，20（7）：1880-1883.

［4］ WIE B.Thrust vector control analysis and design for solar-sail spacecraft［J］.Journal of Spacecraft and Rockets，2007，44（3）：545-557.

［5］ GONG Seng－ping，LI Jun-feng，BAOYIN He-xi.Passive stability design for the solar sail on displaced orbits［J］.Journal of Spacecraft and Rockets，2007，44（5）：1071-1079.

［6］ 李爭學，王本利，馬興瑞.帶撓性附件航天器的魯棒姿態(tài)跟蹤控制［J］.宇航學報，2009，30（3）：974-981.

［7］ LIU Ying－ying，ZHOU Jun.Dynamics model and simulation for flexible satellite with orbit control force［C］／／Proceedings of the 2010 3rd IEEE International Conference on Computer Science and Information Technology （ICCSIT 2010）.Piscataway：IEEE Press，2010，9：227-230.

［8］ SMITH S W，SONG H B.Flexible models for solar sail control［C］／／46th AIAA／ASME／ASCE／AHS／ASC Structures，Structural Dynamics ＆ Materials Conference.Austin：American Inst Aeronautics and Astronautics Inc，2005，1：11-40.

［9］ LIU Jia-fu，RONG Si-yuan，LI Jian-guo，et al.Attitude dynamics modeling and control of large flexible solar sail spacecraft［C］／／ 2010 3rd International Symposium on Systems and Control in Aeronautics and Astronautics （ISSCAA 2010）.Piscataway：IEEE Press，2010：243-247.

［10］ WANG C G，MAO L N，DU X W，et al.Influence parameter analysis and wrinkling control of space membrane structures［J］.Mechanics of advanced materials and structures，2010，17（1）：49-59.

［11］ SAKAMOTO H，MIYAZAKI Y，PARKT K C.Finite element modeling of sail deformation under solar radiation pressure［J］.Journal of Spacecraft and Rockets，2007，44（3）：514-521.

［12］ DAVID M，THOMAS T.Solar sail propulsion modeling［C］／／ 48th AIAA／ASME／ASCE／AHS／ASC Structures，Structural Dynamics，and Materials Conference.Reston： American Institute of Aeronautics and Astronautics Inc，2007，2：1452-1468.

［13］ LI Qing，WANG Tian-shu，MA Xing-rui，Reduced model of flexible sail-boom interaction for solar sail dynamics［C］／／60th International Astronautical Congress 2009（IAC 2009）.Paris：International Astro-nautical Federation，2009，6：5045-5055.

［14］ WIE B，Solar sail attitude control and dynamics（part 1）［J］.Journal of Guidance，Control，and Dynamics，2004.27（4）：526-535.