空間太陽能電站姿態(tài)運(yùn)動(dòng)-結(jié)構(gòu)振動(dòng)耦合建模與分析

穆瑞楠，王藝睿，譚述君,4，吳志剛,4，齊朝暉

(1. 大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國家重點(diǎn)實(shí)驗(yàn)室，大連 116024；2. 中國科學(xué)院國家空間科學(xué)中心，北京 101408；3. 中國科學(xué)院大學(xué)，北京 101407；4. 大連理工大學(xué)航空航天學(xué)院，大連 116024)

0 引 言

空間太陽能電站(SSPS, Space Solar Power Station)是于1968年由美國科學(xué)家Peter Glaser提出的一種超大空間結(jié)構(gòu)概念[1]。因其可實(shí)現(xiàn)連續(xù)工作、能量利用率高等諸多優(yōu)點(diǎn)，近年來越來越受到各國和各研究機(jī)構(gòu)的關(guān)注。不僅相繼提出了數(shù)十種構(gòu)型設(shè)計(jì)方案，在能量傳輸、在軌組裝及動(dòng)力學(xué)分析與控制等各個(gè)相關(guān)領(lǐng)域開展了大量工作[2]。由于尺寸和剛度的限制，這類超大空間結(jié)構(gòu)很難在地面開展實(shí)驗(yàn)驗(yàn)證工作。因此，借助于仿真手段準(zhǔn)確建模并分析其在軌的動(dòng)力學(xué)特性至關(guān)重要。

絕大多數(shù)空間太陽能電站構(gòu)型尺寸巨大[3]，通常設(shè)計(jì)為模塊組裝的方式在軌構(gòu)建，包括桿件、薄膜、接口等基本組件，其建模難度高，模型包含的單元節(jié)點(diǎn)數(shù)目龐大，分析計(jì)算耗時(shí)較長，且不利于進(jìn)行大量的參數(shù)分析。特別是對(duì)于其在軌的動(dòng)力學(xué)分析，需要同時(shí)考慮其軌道運(yùn)動(dòng)和姿態(tài)運(yùn)動(dòng)，使得模型具有較高的非線性，進(jìn)一步加大求解難度并降低求解效率。因此，對(duì)于空間太陽能電站的在軌動(dòng)力學(xué)分析，基于其結(jié)構(gòu)特點(diǎn)，將其簡化為易于建模分析，并且具有一致動(dòng)力學(xué)特性的簡單結(jié)構(gòu)十分必要。由于模塊組裝的在軌構(gòu)建方式，其結(jié)構(gòu)具有周期特性，因此采用基于能量等效原理的連續(xù)體等效方法最為合適。這種方法基于結(jié)構(gòu)的周期性，不會(huì)因基本單元數(shù)量的擴(kuò)充而導(dǎo)致計(jì)算量大幅增加，又能較好刻畫結(jié)構(gòu)的動(dòng)力學(xué)特性?；谀芰康刃г淼倪B續(xù)體等效建模方法已有效應(yīng)用于國際空間站主支撐桁架結(jié)構(gòu)[4]以及大型環(huán)形可展天線的周邊桁架結(jié)構(gòu)[5]。

對(duì)于空間太陽能電站這種具有超大尺寸和超高柔性的空間結(jié)構(gòu)，其尺寸達(dá)到公里量級(jí)，其軌道運(yùn)動(dòng)、姿態(tài)運(yùn)動(dòng)和結(jié)構(gòu)振動(dòng)之間的耦合作用不可忽略。一些學(xué)者對(duì)空間太陽能電站的耦合動(dòng)力學(xué)開展了研究。Wie[6]基于Abacus電站構(gòu)型，分析了環(huán)境干擾力作用下的軌道運(yùn)動(dòng)與姿態(tài)運(yùn)動(dòng)，提出了軌道-姿態(tài)耦合控制方案。吳志剛[7]基于太陽塔式構(gòu)型，考慮地球扁率的引力場，在地球同步拉普拉斯軌道上研究了高階重力和力矩對(duì)空間太陽能電站姿軌運(yùn)動(dòng)的影響，發(fā)現(xiàn)高階力對(duì)軌道運(yùn)動(dòng)影響較大，而對(duì)姿態(tài)運(yùn)動(dòng)影響較小可忽略。同時(shí)，超大結(jié)構(gòu)的結(jié)構(gòu)振動(dòng)引起了更多學(xué)者的關(guān)注。Malla[8]、Ishimura[9]和穆瑞楠[10]建立了軸向振動(dòng)空間啞鈴的軌道-姿態(tài)-振動(dòng)耦合動(dòng)力學(xué)模型，分別討論了初始條件、系統(tǒng)參數(shù)以及結(jié)構(gòu)尺寸對(duì)空間啞鈴動(dòng)力學(xué)特性的影響。在穆瑞楠的研究工作中特別提到結(jié)構(gòu)振動(dòng)引起了姿態(tài)運(yùn)動(dòng)周期的改變，在較大初始姿態(tài)角下可能引起姿態(tài)運(yùn)動(dòng)的翻滾現(xiàn)象。Silva[11-12]和Liu[13]以柔性梁為研究對(duì)象，前者給出了柔性梁的非線性擾動(dòng)方程，并用擾動(dòng)分析方法討論了系統(tǒng)可能出現(xiàn)的共振現(xiàn)象；而后者考慮動(dòng)力剛化效應(yīng)，提出了參數(shù)激勵(lì)模型(PEM)，該模型可以更加精確地描述超大柔性空間結(jié)構(gòu)在持續(xù)轉(zhuǎn)動(dòng)時(shí)的結(jié)構(gòu)振動(dòng)。魏乙[14]將JAXA繩系構(gòu)型簡化為集中質(zhì)量-繩索-板模型，基于絕對(duì)節(jié)點(diǎn)坐標(biāo)法建立了耦合動(dòng)力學(xué)方程，提出了結(jié)合辛Runge-Kutta方法的微分-代數(shù)方程求解方法，分析了系統(tǒng)參數(shù)對(duì)耦合特性的影響并驗(yàn)證了求解方法的有效性。

本文基于能量等效的原則，將空間太陽能電站等效為柔性梁模型，并考慮重力梯度影響，建立姿態(tài)運(yùn)動(dòng)-結(jié)構(gòu)振動(dòng)耦合動(dòng)力學(xué)模型；提出Runge-Kutta法和Newmark法相結(jié)合的耦合動(dòng)力學(xué)方程的高效求解方法；最后對(duì)比不同參數(shù)下的仿真結(jié)果，并對(duì)超大柔性梁模型的在軌耦合特性進(jìn)行討論。

1 連續(xù)體等效建模

基于能量等效的連續(xù)體等效方法既可應(yīng)用于桁架結(jié)構(gòu)，也可應(yīng)用于框架結(jié)構(gòu)，只需要結(jié)構(gòu)具有周期性。多旋轉(zhuǎn)關(guān)節(jié)空間太陽能電站(MR-SPS)構(gòu)型的支撐部分為桁架梁組成的平面框架結(jié)構(gòu)，可分解為周期框架模塊和周期桁架模塊。首先針對(duì)周期桁架模塊，將桁架梁等效為連續(xù)梁模型；然后針對(duì)周期框架模塊，從而得到整體的等效連續(xù)梁模型?；谀芰康刃г淼倪B續(xù)體等效建模的理論方法如下。

周期模塊的各根構(gòu)件簡化為桿或梁模型。定義構(gòu)件局部坐標(biāo)系原點(diǎn)位于構(gòu)件中心，其x1軸方向沿構(gòu)件方向。每根構(gòu)件的變形為沿x1軸的拉伸變形以及繞y1軸和z1軸的彎曲變形，而其運(yùn)動(dòng)為質(zhì)心沿x1軸、y1軸和z1軸的平動(dòng)與繞質(zhì)心轉(zhuǎn)動(dòng)。周期模塊總的變形能和動(dòng)能可表示為

(1)

(2)

將周期模塊內(nèi)各點(diǎn)位移和應(yīng)變在周期模塊整體的幾何中心處Taylor展開，則可將周期模塊任一橫截面上各點(diǎn)的位移和應(yīng)變用周期模塊中心處的位移和轉(zhuǎn)角及應(yīng)變和曲率表示，即

u=Γuuc+Γεεc

(3)

其中uc表示周期模塊中心處的位移與轉(zhuǎn)角，而εc表示其應(yīng)變與曲率，Γu和Γε分別為對(duì)應(yīng)的展開系數(shù)矩陣。將式(3)代入式(1)和式(2)即可得到以中心點(diǎn)位移和應(yīng)變描述的周期模塊的變形能和動(dòng)能。

另一方面，可建立與模塊相同長度的等效梁模型的變形能及動(dòng)能。將等效梁模型在模塊長度內(nèi)的位移及應(yīng)變近似為均勻，則以中心點(diǎn)位移和應(yīng)變描述的變形能和動(dòng)能表達(dá)式分別為

(4)

(5)

其中L為模塊長度。

分別對(duì)比周期模塊和等效梁模型的變形能和動(dòng)能表達(dá)式，可得結(jié)構(gòu)的等效關(guān)系為

(6)

(7)

2 在軌柔性梁的動(dòng)力學(xué)建模

基于哈密頓原理，考慮重力梯度影響，建立柔性梁的姿態(tài)運(yùn)動(dòng)和結(jié)構(gòu)振動(dòng)的耦合動(dòng)力學(xué)模型。將地球視為理想球體，用Oe表示地心，柔性梁在軌道平面內(nèi)繞地球運(yùn)動(dòng)，如圖1所示。

圖1 柔性梁結(jié)構(gòu)在軌運(yùn)行示意圖Fig.1 Schematic diagram of flexible beam structure in orbit

以柔性梁幾何中心的運(yùn)動(dòng)作為其軌道運(yùn)動(dòng)，軌道半徑R由地心指向柔性梁幾何中心，指向近地點(diǎn)的長半軸方向與軌道半徑矢量方向之間的夾角θ為軌道角，并用ωo表示軌道角速度。選取柔性梁幾何中心O為原點(diǎn)，柔性梁幾何中心處切線方向?yàn)閤軸方向，建立柔性梁的固聯(lián)浮動(dòng)坐標(biāo)系，則軌道半徑矢量方向與浮動(dòng)坐標(biāo)系x軸方向的夾角φ為面內(nèi)姿態(tài)角，梁上各點(diǎn)偏移浮動(dòng)坐標(biāo)系x軸的橫向位移為彎曲變形。坐標(biāo)量ρP由幾何中心O到點(diǎn)P變形前位置；橫向變形量wP由點(diǎn)P變形前位置指向變形后位置。則梁上點(diǎn)P相對(duì)于地心的位置rP為

(8)

(9)

其中mP表示點(diǎn)P的處的結(jié)構(gòu)質(zhì)量，即線密度。

對(duì)點(diǎn)P的動(dòng)能式(9)在結(jié)構(gòu)長度上積分，得到整體柔性梁的動(dòng)能為

(10)

對(duì)于在軌柔性梁結(jié)構(gòu)而言，其勢能包含重力勢能和變形能兩部分。這里考慮的柔性梁為純彎曲梁模型，其應(yīng)變能與教材中給出的歐拉梁應(yīng)變能一致[15]。因此柔性梁的總勢能為

(11)

其中μ為地球引力常數(shù)，EI為柔性梁的彎曲剛度，r為梁上某點(diǎn)到地心的距離，表示為

r2=(R2+x2+w2+2xRcosφ-2wRsinφ)

(12)

無非保守力做功的Hamilton原理為

(13)

將柔性梁的動(dòng)能式(10)和勢能式(11)代入式(13)，得到柔性梁的在軌姿態(tài)-振動(dòng)耦合動(dòng)力學(xué)方程為

(14)

(15)

注意到耦合動(dòng)力學(xué)方程中與重力梯度有關(guān)的項(xiàng)含有(r-3)項(xiàng)，該非線性項(xiàng)難以用有限元方法處理，因此需要將其近似展開。由式(12)得(r-3)的二項(xiàng)式展開為

(16)

其中參數(shù)ε=x/R,表示結(jié)構(gòu)尺寸與軌值半徑之比上式近似保留至參數(shù)ε的二次項(xiàng)以及w的線性項(xiàng)。將式代入式和式得到近似耦合動(dòng)力學(xué)方程，其中含有參數(shù)ε的一階項(xiàng)和二階項(xiàng)分別對(duì)應(yīng)于重力梯度力/力矩的一階項(xiàng)和二階項(xiàng)。

由于柔性梁的在軌耦合動(dòng)力學(xué)方程既包括偏微分方程，又含有在結(jié)構(gòu)長度上的積分項(xiàng)，因此采用有限元方法將結(jié)構(gòu)離散，引入只與時(shí)間t有關(guān)的節(jié)點(diǎn)位移變量和只與位置變量x有關(guān)的形函數(shù)，從而消去方程中關(guān)于位置變量x的微分和積分。梁上某點(diǎn)彎曲變形重新表述為

w(x,t)=N(x)ae(t)

(17)

其中N(x)表示單元內(nèi)的形函數(shù)，ae(t)表示單元內(nèi)的節(jié)點(diǎn)位移，這里梁單元采用經(jīng)典的二節(jié)點(diǎn)Hermite單元?；谟邢拊碚揫15]，利用梁橫向變形的新形式，可將在軌耦合動(dòng)力學(xué)方程式(14)和式(13)分別轉(zhuǎn)化為

(18)

(19)

其中J為柔性梁面內(nèi)轉(zhuǎn)動(dòng)慣量，M和K分別為柔性梁質(zhì)量陣和剛度陣，TG1和TG2分別為重力梯度力矩的一階項(xiàng)和二階項(xiàng)，PG1和PG2分別為重力梯度橫向分力的一階項(xiàng)和二階項(xiàng)。各項(xiàng)的表達(dá)式如下

其中Φ0、Φ1和Φ2是用于簡化標(biāo)記的矩陣和向量，其表達(dá)式為

3 模型參數(shù)及求解方法

多旋轉(zhuǎn)關(guān)節(jié)空間太陽能電站示意圖以及其中作支撐骨架的結(jié)構(gòu)力學(xué)模型如圖2(a)所示，其整體結(jié)構(gòu)由具有周期性的框架結(jié)構(gòu)組成，而周期框架結(jié)構(gòu)中的單根構(gòu)件又是由具有周期性的桁架結(jié)構(gòu)組成。桁架周期模塊和框架周期模塊以及它們的整體坐標(biāo)系如圖2(b)和(c)中藍(lán)色虛線(見電子版)框所示。

多旋轉(zhuǎn)關(guān)節(jié)空間太陽能電站構(gòu)型的幾何尺寸如表1所示，各桿件選用T700碳纖維復(fù)合材料，其材料參數(shù)在表1中給出。

利用基于能量等效的連續(xù)體等效方法，將該構(gòu)型等效為矩形等截面的均勻質(zhì)量柔性梁模型，其彎曲剛度EI= 1.56×109N·m2，密度與橫截面積的乘積(線密度)ρA= 6.20 kg/m。等效梁模型的結(jié)構(gòu)基頻為4.058 × 10-4Hz，而多旋轉(zhuǎn)關(guān)節(jié)構(gòu)型的有限元模型的基頻為3.841 × 10-4Hz，二者不僅具有相同量級(jí)，且差別僅為5%，能夠反應(yīng)其動(dòng)力學(xué)現(xiàn)象。

圖2 桁架和框架周期模塊示意圖Fig.2 Diagram of truss and frame periodic modules

桁架模塊尺寸1 m×1 m×1 m框架模塊尺寸210 m×310 m總長度11.8 km桿件截面積2.359 cm2桿件彈性模量248 GPa桿件密度1750 kg/m3

涉及到的常數(shù)有：地球平均半徑RE=6378 km，地球引力常數(shù)為μ=3986 00 km3/s2。這里主要研究大尺寸結(jié)構(gòu)在低軌道(LEO)和地球同步軌道(GEO)下的動(dòng)力學(xué)現(xiàn)象。選取柔性梁所處的低軌道的高度為322 km(即幾何中心軌道半徑R=6700 km)，此時(shí)的軌道頻率為1.833×10-4Hz，而地球同步軌道的高度為35 786 km(即幾何中心軌道半徑R=42 164 km)，此時(shí)的軌道頻率為1.161×10-5Hz. 同時(shí)，結(jié)構(gòu)尺寸是結(jié)構(gòu)基頻的重要影響因素，因此在50 m～10 000 m之間選取一系列值作為對(duì)比。初始姿態(tài)角對(duì)姿態(tài)運(yùn)動(dòng)幅度有很大影響，這里取兩種情況：小姿態(tài)角(φ0=0.1 rad)和大姿態(tài)角(φ0=π/4 rad)；其他初始條件均設(shè)置為零。

在求解耦合方程式和式時(shí)，由于其已被轉(zhuǎn)化為二階常微分方程組，因此可采用經(jīng)典的Runge-Kutta 法求解。但由于有限元離散后的節(jié)點(diǎn)較多，使得方程維數(shù)較大，因此直接采用這種方法求解效率低。由于姿態(tài)和振動(dòng)方程中系數(shù)相互包含，耦合求解時(shí)具有強(qiáng)非線性，因此難以使用常用的有限元方程求解方法Newmark法。對(duì)于非線性微分代數(shù)方程缺乏有效且高效的求解方法是超大柔性空間結(jié)構(gòu)耦合動(dòng)力學(xué)研究所面臨的問題之一。本文將兩種方法結(jié)合，提出了Runge-Kutta 法和Newmark法同時(shí)迭代的方法，具體過程如圖3所示。與應(yīng)用Runge-Kutta法求解耦合方程相比，本文提出的方法在保證求解精度的同時(shí)效率大幅提升。在低軌道、小初始姿態(tài)角、結(jié)構(gòu)尺寸為1000 m的條件下，仿真時(shí)長1000 min，時(shí)間步長為0.1 s，本文提出的方法在效率上提高約752倍，兩種方法得出的響應(yīng)相差5.36%. 后續(xù)仿真結(jié)果均基于本文提出的改進(jìn)方法得出，其中時(shí)間步長取為0.1 s。

圖3 改進(jìn)的耦合模型求解算法流程圖Fig.3 Flow chart of improved algorithm for coupling model

在利用提出的方法求解動(dòng)力學(xué)響應(yīng)時(shí)，姿態(tài)運(yùn)動(dòng)方程式(18)中含有未知的初始加速度響應(yīng)，無法開始迭代。因此利用振動(dòng)方程式(19)簡化姿態(tài)方程得

(20)

其中

4 姿態(tài)-振動(dòng)耦合特性

4.1 結(jié)構(gòu)振動(dòng)幅值影響分析

從結(jié)構(gòu)振動(dòng)方程式(19)可以看出，姿態(tài)運(yùn)動(dòng)對(duì)結(jié)構(gòu)振動(dòng)既有慣性力的直接影響，也有改變重力梯度分力的間接影響。若忽略結(jié)構(gòu)變形對(duì)姿態(tài)運(yùn)動(dòng)的影響，則姿態(tài)運(yùn)動(dòng)方程式(18)可簡化為

(21)

將式(21)代入結(jié)構(gòu)振動(dòng)方程(19)中發(fā)現(xiàn)，方程右端第一項(xiàng)(角加速度項(xiàng))與第二項(xiàng)(重力梯度一階項(xiàng))的一部分相互抵消，則結(jié)構(gòu)振動(dòng)方程整理為

(22)

由上式可以看出，重力梯度的二階項(xiàng)PG2是激發(fā)結(jié)構(gòu)振動(dòng)的主要因素，姿態(tài)運(yùn)動(dòng)通過重力梯度的二階項(xiàng)PG2對(duì)結(jié)構(gòu)振動(dòng)的間接影響。圖4分別給出了不同軌道下的大幅度姿態(tài)運(yùn)動(dòng)下的結(jié)構(gòu)振動(dòng)響應(yīng)。由重力梯度的二階項(xiàng)PG2的表達(dá)式可知，該項(xiàng)大小隨軌道半徑的增加反比例變化，同時(shí)以軌道角速度的平方量級(jí)減小，因此不同軌道下的振動(dòng)量級(jí)近似相差3個(gè)量級(jí)。

圖4 大幅度姿態(tài)運(yùn)動(dòng)下的結(jié)構(gòu)振動(dòng)響應(yīng)Fig.4 Structural vibration response under attitude motion with large magnitude

由前文中的重力梯度的二階項(xiàng)PG2表達(dá)式可以看出，該項(xiàng)與質(zhì)量陣的比值以結(jié)構(gòu)尺寸的二次方量級(jí)增加。同時(shí)，剛度陣與質(zhì)量陣的比值以結(jié)構(gòu)尺寸的四次方量級(jí)減小。因此當(dāng)結(jié)構(gòu)基頻遠(yuǎn)大于姿態(tài)運(yùn)動(dòng)頻率時(shí)，振幅以結(jié)構(gòu)尺寸的6次方量級(jí)增加。

圖5 結(jié)構(gòu)尺寸與振動(dòng)幅度變化關(guān)系(LEO)Fig.5 Relationship between structural vibration magnitude and structure size(LEO)

圖5給出了低軌道下的結(jié)構(gòu)振幅隨結(jié)構(gòu)尺寸變化結(jié)果，在尺寸較小時(shí)，結(jié)構(gòu)振幅量級(jí)變化規(guī)律與理論結(jié)果一致，即按結(jié)構(gòu)尺寸的6次方量級(jí)增加。但仿真發(fā)現(xiàn)，結(jié)構(gòu)尺寸過大會(huì)引起系統(tǒng)不穩(wěn)定現(xiàn)象。以本文的參數(shù)為例，低軌道下結(jié)構(gòu)尺寸接近9450 m時(shí)，結(jié)構(gòu)振幅增長速度加快，最終系統(tǒng)運(yùn)動(dòng)發(fā)散。這種不穩(wěn)定現(xiàn)象對(duì)應(yīng)的臨界尺寸與結(jié)構(gòu)柔性和軌道高度有關(guān)，軌道高度增加則對(duì)應(yīng)的臨界尺寸增大。

4.2 結(jié)構(gòu)振動(dòng)頻率影響分析

由簡化的結(jié)構(gòu)振動(dòng)方程式(22)可以得到結(jié)構(gòu)振動(dòng)的第一階頻率為

(23)

其中ωs為結(jié)構(gòu)基頻。從式(23)可以看出，結(jié)構(gòu)振動(dòng)頻率受到轉(zhuǎn)動(dòng)耦合項(xiàng)(等式右端第二項(xiàng))以及重力梯度項(xiàng)(等式右端第三項(xiàng))的影響，使結(jié)構(gòu)振動(dòng)周期出現(xiàn)波動(dòng)現(xiàn)象。

考慮到在無控情況下姿態(tài)運(yùn)動(dòng)接近于正弦振動(dòng)，因此將姿態(tài)運(yùn)動(dòng)近似表示為φ=φocos(κωot)，其中φ0為初始姿態(tài)角，kωo為姿態(tài)運(yùn)動(dòng)頻率，k與φ0的取值相關(guān)，可通過與仿真得到的姿態(tài)運(yùn)動(dòng)響應(yīng)擬合得到。將近似的姿態(tài)運(yùn)動(dòng)代入式(23)并無量綱化得

1-3sin2(φ0coskωot)

(24)

注意到，式(24)中的轉(zhuǎn)動(dòng)耦合項(xiàng)和重力梯度項(xiàng)中的軌道角速度處于三角函數(shù)內(nèi)，因此軌道角速度主要對(duì)頻率比的變化周期有影響，對(duì)幅值影響較小。

圖6 姿態(tài)運(yùn)動(dòng)和重力梯度對(duì)振動(dòng)頻率的影響(φ0 = π/4 rad)Fig.6 Effects of attitude motion and gravity gradient on frequency of structural vibration (φ0 = π/4 rad)

圖6給出了低軌道大幅姿態(tài)運(yùn)動(dòng)下的轉(zhuǎn)動(dòng)耦合項(xiàng)和重力梯度項(xiàng)對(duì)結(jié)構(gòu)振動(dòng)頻率隨時(shí)間變化的影響曲線。對(duì)于大幅姿態(tài)運(yùn)動(dòng)(φ0=π/4 rad)，對(duì)應(yīng)的k=1.5003，從圖6可以看出，重力梯度項(xiàng)波動(dòng)幅度增大至一倍軌道頻率，時(shí)而使結(jié)構(gòu)振動(dòng)頻率增大，時(shí)而使其減小。轉(zhuǎn)動(dòng)耦合項(xiàng)仍使結(jié)構(gòu)振動(dòng)頻率降低，其量級(jí)增大至兩倍軌道頻率。小幅姿態(tài)運(yùn)動(dòng)下的結(jié)果類似，但幅值較小。

當(dāng)結(jié)構(gòu)基頻極低時(shí)，由式(23)可知結(jié)構(gòu)可能發(fā)生屈曲的失穩(wěn)現(xiàn)象，文獻(xiàn)[13]中對(duì)慣性穩(wěn)定的柔性梁的分析中也提及了這種現(xiàn)象。將令式(23)左側(cè)最小值等于0的結(jié)構(gòu)基頻稱為臨界結(jié)構(gòu)基頻。則對(duì)于不同幅度姿態(tài)運(yùn)動(dòng)，對(duì)應(yīng)的臨界結(jié)構(gòu)基頻也不同。臨界基頻ωb滿足

1+3sin2(φ0coskωot)}

(25)

利用正弦函數(shù)的泰勒展開，以及二項(xiàng)式展開公式，式(25)右端括號(hào)內(nèi)的表達(dá)式可以轉(zhuǎn)化為參數(shù)λ的四階多項(xiàng)式形式。

(26)

其中λ=sinkωot。對(duì)式(20)求導(dǎo)并分析發(fā)現(xiàn)其導(dǎo)數(shù)表達(dá)式只有一個(gè)實(shí)根，同時(shí)該式的最高階的四次項(xiàng)系數(shù)為負(fù)，因此在全局有且只有一個(gè)極大值。對(duì)其導(dǎo)數(shù)表達(dá)式應(yīng)用一元三次方程求根公式，得到該極值點(diǎn)為

(27)

同時(shí)注意到λ在[-1,1]范圍內(nèi)取值，若極值點(diǎn)不在取值范圍內(nèi)則最大值點(diǎn)在邊界取值。經(jīng)數(shù)值驗(yàn)證，在初始姿態(tài)角φ0在[0 , 0.476π]時(shí)，最大值在λ = -1處取得，則臨界振動(dòng)頻率ωb滿足

(28)

圖7 臨界振動(dòng)頻率隨初始姿態(tài)角變化曲線Fig.7 Variation of critical frequency of structural vibration with respect to initial attitude angle

圖7根據(jù)式(28)給出了低軌道和地球同步軌道下的臨界振動(dòng)頻率隨初始姿態(tài)角變化曲線，發(fā)現(xiàn)不論在低軌道還是地球同步軌道，耦合作用對(duì)結(jié)構(gòu)振動(dòng)頻率的影響都存在。初始姿態(tài)角在0.386π rad時(shí)臨界振動(dòng)頻率最大，此時(shí)影響幅值約為2.11倍軌道頻率。因此，在設(shè)計(jì)梁式空間結(jié)構(gòu)時(shí)，考慮耦合作用對(duì)結(jié)構(gòu)振動(dòng)頻率的影響，才能保證結(jié)構(gòu)不會(huì)發(fā)生屈曲不穩(wěn)定現(xiàn)象。

4.3 結(jié)構(gòu)振動(dòng)對(duì)姿態(tài)運(yùn)動(dòng)的影響

若引入小角度假設(shè)條件，并忽略結(jié)構(gòu)振動(dòng)對(duì)轉(zhuǎn)動(dòng)慣量的影響，則姿態(tài)運(yùn)動(dòng)方程簡化為

(29)

其中

僅考慮柔性梁一階振型的振動(dòng)，令

a=w0Ω(x)sin(ωst)

(30)

其中w0為最大振動(dòng)幅值，Ω(x)為一階振型。將式(30)代入式(29)，得到姿態(tài)運(yùn)動(dòng)有系統(tǒng)阻尼的自振頻率ωa為

(31)

其中阻尼比為

(32)

阻尼比的最大值可表示為

(33)

其中α是系統(tǒng)常數(shù)，與材料屬性、梁振型等參數(shù)有關(guān)。由式可以看出，姿態(tài)運(yùn)動(dòng)頻率受結(jié)構(gòu)振動(dòng)幅度和結(jié)構(gòu)尺寸的影響。結(jié)構(gòu)振動(dòng)幅度增大時(shí)，姿態(tài)運(yùn)動(dòng)頻率降低。同時(shí)注意到結(jié)構(gòu)振動(dòng)幅度與結(jié)構(gòu)尺寸的6次方量級(jí)關(guān)系，因此結(jié)構(gòu)尺寸增加也導(dǎo)致姿態(tài)運(yùn)動(dòng)頻率降低。

圖8分別給出了小幅度姿態(tài)運(yùn)動(dòng)(φ0=0.1 rad)和大幅度姿態(tài)運(yùn)動(dòng)(φ0= π/4 rad)下的剛體與柔性體姿態(tài)角響應(yīng)對(duì)比。為了顯示出柔性引起的姿態(tài)周期變化，圖中截取了仿真時(shí)間3200～3300 min的結(jié)果。從圖中可以發(fā)現(xiàn)，隨著仿真時(shí)間增加，柔性體姿態(tài)運(yùn)動(dòng)越來越偏離剛體的姿態(tài)運(yùn)動(dòng)，即柔性體姿態(tài)運(yùn)動(dòng)周期較大。對(duì)比兩圖發(fā)現(xiàn)，初始姿態(tài)角越大，這種結(jié)構(gòu)柔性引起的姿態(tài)運(yùn)動(dòng)周期增大現(xiàn)象越明顯，這是由于初始姿態(tài)角越大，姿態(tài)運(yùn)動(dòng)引起的結(jié)構(gòu)振動(dòng)幅度越大，反過來使得姿態(tài)運(yùn)動(dòng)阻尼比越大，因此姿態(tài)運(yùn)動(dòng)周期增大就越明顯。

圖8 柔性體與剛體的姿態(tài)角響應(yīng)比較(LEO)Fig.8 Comparison of attitude angle between flexible and rigid model(LEO)

圖9給出了地球同步軌道下的大幅度姿態(tài)運(yùn)動(dòng)時(shí)的剛?cè)峤Y(jié)構(gòu)姿態(tài)角響應(yīng)結(jié)果對(duì)比，結(jié)構(gòu)柔性引起的姿態(tài)運(yùn)動(dòng)周期增大現(xiàn)象仍然存在。分析式(33)可知，軌道角速度降低原本是會(huì)增大姿態(tài)運(yùn)動(dòng)阻尼比，但是由圖4可知此時(shí)結(jié)構(gòu)振動(dòng)幅度大幅下降，因此在地球同步軌道下剛體和柔性體的姿態(tài)運(yùn)動(dòng)周期差別較小。

圖9 柔性體與剛體的姿態(tài)角響應(yīng)比較(GEO)Fig.9 Comparison of attitude angle between the flexible and rigid model(GEO)

5 結(jié)論與展望

本文將MR-SPS電站構(gòu)型等效為柔性梁模型，并考慮重力梯度影響，建立了姿態(tài)運(yùn)動(dòng)和結(jié)構(gòu)振動(dòng)耦合動(dòng)力學(xué)模型，利用改進(jìn)算法分析了不同參數(shù)取值下兩者耦合關(guān)系。得到主要結(jié)論如下：1)所提出的基于Runge-Kutta 法和Newmark法的改進(jìn)算法適用于姿態(tài)運(yùn)動(dòng)-結(jié)構(gòu)振動(dòng)耦合方程求解，大幅提高了計(jì)算效率；2)在軌柔性梁模型中的重力梯度項(xiàng)至少近似保留至ε的二階展開項(xiàng)，尺寸過大可能導(dǎo)致結(jié)構(gòu)屈曲不穩(wěn)定；3)姿態(tài)運(yùn)動(dòng)和結(jié)構(gòu)振動(dòng)耦合效應(yīng)導(dǎo)致結(jié)構(gòu)振動(dòng)頻率降低，姿態(tài)運(yùn)動(dòng)周期增大，在GEO軌道仍有這種現(xiàn)象。為保證系統(tǒng)在軌穩(wěn)定性，超大型空間結(jié)構(gòu)設(shè)計(jì)時(shí)應(yīng)考慮耦合因素影響，并需要進(jìn)一步開展與結(jié)構(gòu)振動(dòng)抑制協(xié)同的姿態(tài)控制策略研究。

參 考 文 獻(xiàn)

[1] Glaser P E. Power from the sun: its future[J]. Science, 1968, 162(3856):857-861.

[2] 侯欣賓, 王立. 空間太陽能電站技術(shù)發(fā)展現(xiàn)狀及展望[J]. 中國航天, 2015, (02):12-15. [Hou Xin-bin, Wang Li. Present situation and prospect for the solar power station [J]. Space System and Technology, 2015, (02):12-15.]

[3] 侯欣賓, 王立, 張興華等. 多旋轉(zhuǎn)關(guān)節(jié)空間太陽能電站概念方案設(shè)計(jì)[J]. 宇航學(xué)報(bào), 2015, (11): 1332-1338. [Hou Xin-bin, Wang Li, Zhang Xing-hua, et al. Concept design on multi-rotary joints SPS [J]. Journal of Astronautics, 2015, (11): 1332-1338.]

[4] Noor A K, Mikulas M M. Continuum modeling of large lattice structures: status and projections[M]// Atluri S N, Anthony K A. Large space structures: dynamics and control. Heidelberg, Berlin: Springer, 1988:1-34.

[5] 劉福壽, 金棟平. 環(huán)形桁架結(jié)構(gòu)徑向振動(dòng)的等效圓環(huán)模型[J]. 力學(xué)學(xué)報(bào), 2016, 48(5):1184-1191. [Liu Fu-shou, Jin Dong-ping. Equivalent circular ring model for the radial vibration analysis of hoop truss structures[J]. Chinese Journal of Theoretical and Applied Mechanics, 2016, 48(5):1184 - 1191.]

[6] Wie B, Roithmayr C M. Attitude and orbit control of a very large geostationary solar power satellite[J]. Journal of Guidance, Control, and Dynamics. 2005, 28(3):439-451.

[7] 吳志剛, 劉玉亮, 張開明, 等. 高階重力和力矩對(duì)空間太陽能電站運(yùn)動(dòng)的影響[J]. 空間控制技術(shù)與應(yīng)用, 2016, 42(4):1-5. [Wu Zhi-gang, Liu Yu-liang, Zhang Kai-ming, et al. The influences of higher order gravitational force and torque to the motion of space solar power station [J]. Aerospace Control and Application, 2016, 42(4):1-5.]

[8] Malla R B. Structural and orbital conditions on response of large space structures[J]. Journal of Aerospace Engineering, 1993, 6(2):115-132.

[9] Ishimura K, Higuchi K. Coupling among pitch motion, axial vibration, and orbital motion of large space structures[J]. Journal of Aerospace Engineering, 2008, 21(2):61-71.

[10] 穆瑞楠, 譚述君, 吳志剛. 重力梯度對(duì)超大柔性空間結(jié)構(gòu)在軌動(dòng)力學(xué)特性的影響[J]. 國防科技大學(xué)學(xué)報(bào), 2017, 39(3):7-14. [Mu Rui-nan, Tan Shu-jun, Wu Zhi-gang. Effect of gravity gradient on dynamical characteristics of very large flexible space structures in orbit[J]. Journal of National University of Defense Technology, 2017, 39(3):7- 14.]

[11] Dasilva M M, Zaretzky C L. Nonlinear dynamics of a flexible beam in a central gravitational field—I. equations of motion[J]. International Journal of Solids and Structures. 1993, 30(17):2287-2299.

[12] Dasilva M M, Zaretzky C L. Nonlinear dynamics of a flexible beam in a central gravitational field—II. nonlinear motions in circular orbit[J]. International Journal of Solids and Structures. 1993, 30(17): 2301-2316.

[13] Liu Y, Wu S, Zhang K, et al. Parametrical excitation model for rigid-flexible coupling system of solar power satellite[J]. Journal of Guidance, Control, and Dynamics. 2017:1-8.

[14] 魏乙, 鄧子辰, 李慶軍,等. 繩系空間太陽能電站動(dòng)力學(xué)響應(yīng)分析[J]. 宇航學(xué)報(bào), 2016, 37(9): 1041-1048. [Wei Yi, Deng Zi-Chen, Li Qing-jun, et al. Analysis of dynamic response of tethered space solar power station [J]. Journal of Astronautics, 2016, 37(9):1041-1048.]

[15] 王勖成. 有限單元法[M]. 北京: 清華大學(xué)出版社, 2003:306-311.