彈性飛行器飛行動力學(xué)建模與剛?cè)狁詈想A次分析

安 陽，王天舒

（清華大學(xué)航天航空學(xué)院，北京 100084）

0 引言

隨著飛行器飛行速度的提高、結(jié)構(gòu)尺寸的增大以及結(jié)構(gòu)重量的降低，柔性結(jié)構(gòu)在飛行器整體結(jié)構(gòu)中的占比越來越大［1-3］，彈性振動對飛行器動力學(xué)與控制的影響也隨之增強［4-6］。例如，高超聲速飛行器為了滿足高速和大射程的要求，一般具有較大的長細(xì)比和較輕的結(jié)構(gòu)重量，導(dǎo)致其本體彈性模態(tài)的固有頻率較低，表現(xiàn)出明顯的柔性特性［7-8］。空間站太陽能電池陣列等大型網(wǎng)架式空間結(jié)構(gòu)由于含有桿、板、繩索等柔性構(gòu)件以及它們之間的非線性約束［9］，不再符合傳統(tǒng)柔性航天器模型中具有中心剛體平臺、柔性結(jié)構(gòu)僅在附件中存在的假設(shè)［10］。此外，近年來提出的太陽帆航天器［11］、空間太陽能電站［12］等新概念航天器，大型柔性結(jié)構(gòu)也占其整體的絕大部分。針對彈性飛行器的動力學(xué)建模及模型降階問題，國內(nèi)外學(xué)者已進行了大量研究［13-17］。楊輝等［18］指出了混合坐標(biāo)法的不足，并基于Hamilton 變分原理和有限元法導(dǎo)出了剛?cè)狁詈舷到y(tǒng)一次近似的動力學(xué)模型，該模型充分考慮了中心剛體大范圍運動與柔性附件彈性振動的耦合效應(yīng)。郭東等［19］基于一種瞬態(tài)坐標(biāo)系，通過拉格朗日方程和有限元法，推導(dǎo)了能夠描述彈性飛行器結(jié)構(gòu)、流場與控制等多學(xué)科耦合特性的動力學(xué)模型。Cao等［20］基于哈密頓原理建立了具有側(cè)向太陽電池陣的三軸穩(wěn)定航天器的剛?cè)狁詈蟿恿W(xué)模型，采用全局模態(tài)法求解了模型的固有頻率和全局模態(tài)振型，并通過與有限元軟件的結(jié)果對比，說明了該方法可以準(zhǔn)確描述航天器的剛?cè)狁詈线\動。袁秋帆等［21］通過定義廣義全局模態(tài)振型，基于哈密頓原理推導(dǎo)了一種單翼大撓性航天器的全局模態(tài)動力學(xué)模型，并采用瑞利瑞茲法得到了非約束模態(tài)頻率和振型的計算方法。閆玉龍等［22］基于D’Alembert 原理推導(dǎo)了航天器柔性附件的振動方程，通過準(zhǔn)坐標(biāo)Lagrange 方程建立了剛-柔-晃耦合的航天器狀態(tài)方程，并指出柔性附件的安裝位置對于航天器的動力學(xué)特性具有重要影響。張恒浩等［23］基于正交理論建立了撓性航天器剛?cè)狁詈舷到y(tǒng)的離散動力學(xué)模型，闡明了動力學(xué)剛化現(xiàn)象產(chǎn)生的原因和對航天器產(chǎn)生干擾的干擾源，并提出了一種能夠抑制這種干擾的一階動力學(xué)模型。曹芊等［24］采用李群SE（3）上的指數(shù)坐標(biāo)描述了航天器的位置與姿態(tài)，在SE（3）框架下推導(dǎo)了柔性航天器姿-軌-結(jié)構(gòu)一體化的動力學(xué)模型，進而得到了其相對動力學(xué)模型。Huang等［25］采用有限元法建立了柔性桁架兩端附著剛體類航天器的動力學(xué)模型，通過分析得到了柔性航天器的無約束模態(tài)。朱尊紅等［26］提出了一種自旋撓性航天器帆板的梁式簡化模型，基于歐拉方程和哈密頓原理建立了航天器的動力學(xué)方程，并解釋了剛性模態(tài)和彈性模態(tài)之間的耦合。Kim等［27］將一種基于單元連通性的剛度評估方法與絕對節(jié)點坐標(biāo)法相結(jié)合，提出了一種基于剛度評估的柔性多體系統(tǒng)動力學(xué)模型降階方法，提高了傳統(tǒng)建模方法的計算效率。

在目前的彈性飛行器動力學(xué)建模領(lǐng)域，對描述彈性體變形運動的變形場階次和模態(tài)截斷階次等問題已有較多研究［18，28-29］，而對在模態(tài)階數(shù)與彈性變形場階次選定后，在動力學(xué)模型推導(dǎo)過程中產(chǎn)生的各階彈性變形相關(guān)項取舍問題的研究較少。在傳統(tǒng)的中心剛體-柔性附件模型中，中心剛體的質(zhì)量、轉(zhuǎn)動慣量等在飛行器整體中占有絕對優(yōu)勢，柔性附件的彈性振動對飛行器整體動力學(xué)的影響較小，因此通常忽略動力學(xué)模型推導(dǎo)過程中產(chǎn)生的關(guān)于彈性變形的高階項，而近似后的模型仍能保有滿足工程要求的精度。但當(dāng)飛行器主體呈現(xiàn)柔性特性時，這種高階項是否還能忽略及其對整體動力學(xué)的影響尚無定論。此外，在動力學(xué)模型建立后，通常選用無外力作用時系統(tǒng)的線動量和角動量作為守恒量，檢驗?zāi)Ｐ偷恼_性和對比模型的計算精度。系統(tǒng)動量表達式中彈性變形相關(guān)項的保留階次也可能對其計算結(jié)果的準(zhǔn)確性產(chǎn)生影響。研究彈性變形相關(guān)項的保留程度對飛行器動力學(xué)及動量計算的影響，對動力學(xué)模型的降階與驗證具有重要的指導(dǎo)意義。本文首先基于虛功率原理與混合坐標(biāo)法建立彈性飛行器的飛行動力學(xué)模型，并定義剛?cè)狁詈想A次以描述動力學(xué)模型推導(dǎo)過程中產(chǎn)生的彈性變形相關(guān)項的保留程度，然后以由柔性中心體和剛性附件組成的剛?cè)狁詈舷到y(tǒng)為例，對比分析動力學(xué)模型及動量公式的剛?cè)狁詈想A次對動力學(xué)響應(yīng)與守恒量的影響。

1 基于虛功率原理的彈性飛行器動力學(xué)建模

假設(shè)飛行器本體為彈性體，舵面為剛體，二者通過單自由度旋轉(zhuǎn)副連接。飛行器系統(tǒng)的構(gòu)型及坐標(biāo)系定義示意圖如圖1 所示。其中，慣性系Oxyz固定在地面，本體系Ocxcyczc為與本體固連的浮動坐標(biāo)系，第i(i=1，…，np)個舵面的隨體系Op，ixp，iyp，izp，i與舵面固連且原點位于舵面與本體的連接處。

圖1 飛行器系統(tǒng)構(gòu)型與坐標(biāo)系定義示意圖Fig.1 Configuration of the flight vehicle system and coordinate system definition diagram

虛功率原理的表達式為：

式中：δPI為慣性力虛功率，δPa為主動力虛功率，δPe為虛應(yīng)變能變化率，δPd為阻尼力虛功率，符號Σ 表示對系統(tǒng)求和。

基于虛功率原理推導(dǎo)彈性飛行器飛行動力學(xué)方程的主要步驟為：首先，計算飛行器本體與舵面等附件的速度、加速度與速度變分；然后，計算本體及各附件的虛功率；最后，將各部件虛功率代入虛功率原理的表達式，根據(jù)各變分相互獨立得到飛行器的動力學(xué)方程。

采用有限元法將飛行器本體離散為l個質(zhì)量單元，則本體上第k(k=1，…，l)個質(zhì)量單元相對慣性系原點的矢徑為：

式中：Rc為本體系原點相對慣性系原點的矢徑，表示飛行器的平動位移為本體未發(fā)生變形時，其第k個質(zhì)量單元相對本體系原點的矢徑為本體第k個質(zhì)量單元的彈性變形矢量，利用模態(tài)展開法表示為：

式中：rcp，0為本體未發(fā)生變形時舵面系原點相對本體系原點的矢徑；rp為舵面上任一質(zhì)量單元相對舵面系原點的矢徑。式（5）體現(xiàn)了本體彈性振動對舵面運動的影響，以下記rcp=rcp，0+

基于式（2）與式（5）計算飛行器本體及舵面的速度、加速度、速度變分與虛功率，代入式（1）并根據(jù)變分、δωc、相互獨立，可得：

飛行器平動動力學(xué)方程為：

式中：m為飛行器的總質(zhì)量；S為飛行器相對本體系的總靜矩；ωc為本體相對慣性系的轉(zhuǎn)動角速度，即飛行器的姿態(tài)角速度；F為飛行器主動力的合力；Qt為飛行器平動耦合慣性力的合力；耦合系數(shù)陣Ct的表達式為：

式中：Cct為本體彈性振動對飛行器平動的耦合系數(shù)陣；Crt，i為第i個舵面受本體振動影響而產(chǎn)生的對飛行器平動的附加耦合系數(shù)陣；mc，k為本體上第k個質(zhì)量單元的質(zhì)量；mp，i為第i個舵面的質(zhì)量；np為舵面的數(shù)量。

飛行器姿態(tài)動力學(xué)方程為：

式中：J為飛行器相對本體系的總慣量張量；T為飛行器主動力矩的合力矩；Qr為飛行器轉(zhuǎn)動耦合慣性力的合力；耦合系數(shù)陣Cr的表達式為：

式中：rpp，i為第i個舵面的質(zhì)心相對其隨體系原點的矢徑；Ccr為本體彈性振動對飛行器轉(zhuǎn)動的耦合系數(shù)陣；Crr，i為第i個舵面受本體振動影響而產(chǎn)生的對飛行器轉(zhuǎn)動的附加耦合系數(shù)陣。需要說明的是：為了方便表示，在式（9）中仍使用了空間矢量的叉乘符號來描述位置矢量與模態(tài)矢量陣之間的運算關(guān)系。它表示將位置矢量分別與模態(tài)矢量陣中的各階模態(tài)矢量進行矢量叉乘，再將叉乘得到的結(jié)果矢量按順序排列成一行，得到一個新的矢量陣。

飛行器振動動力學(xué)方程為：

式中：Cc為飛行器本體的模態(tài)阻尼陣，Kc為模態(tài)剛度陣；Ff為飛行器的模態(tài)力；Qf為飛行器振動耦合慣性力的合力；模態(tài)質(zhì)量陣Mf的表達式為：

式中：Mcf為飛行器本體的模態(tài)質(zhì)量陣，當(dāng)本體的模態(tài)矢量陣采用模態(tài)質(zhì)量規(guī)一時，Mcf近似為單位陣；它是第i個舵面受本體振動影響而產(chǎn)生的附加模態(tài)質(zhì)量陣，一般不是對角陣。則矩陣Mf的非對角線位置也將出現(xiàn)非零元素，使方程（10）中經(jīng)主坐標(biāo)變換后的廣義坐標(biāo)也不再完全解耦。

飛行器的線動量表達式為：

式中：mc為飛行器本體的質(zhì)量；Sc為飛行器本體相對本體系的靜矩，Sp，i為第i個舵面相對舵面系的靜矩；ωp，i為第i個舵面相對本體的轉(zhuǎn)動角速度。

飛行器的角動量表達式為：

式中：pc為飛行器本體的線動量，pp，i為第i個舵面的線動量；Jc為飛行器本體相對本體系的慣量張量，Jp，i為第i個舵面相對舵面系的慣量張量。

為了檢驗所建立的動力學(xué)模型的正確性，采用MATLAB 語言編寫其仿真程序，并與ADAMS 軟件進行仿真結(jié)果對比。當(dāng)本體與舵面均不受外力作用，僅對兩者的相對轉(zhuǎn)角施加驅(qū)動使其按照正弦規(guī)律變化時，自編程序與ADAMS 軟件的部分仿真結(jié)果對比如圖2 所示。圖中兩組曲線吻合較好，驗證了動力學(xué)模型及其仿真程序的正確性。

圖2 自編程序與ADAMS的仿真結(jié)果對比Fig.2 Comparison of simulation results between self-programmed program and ADAMS

2 剛?cè)狁詈想A次的定義

上文推導(dǎo)的動力學(xué)模型是未忽略各階彈性變形相關(guān)項的完整模型，其動力學(xué)方程與系統(tǒng)角動量表達式中含有模態(tài)坐標(biāo)的一階與二階項，將此動力學(xué)模型定義為二次剛?cè)狁詈夏Ｐ停瑒恿抗蕉x為二次動量公式。在二次模型中，本體與舵面相對本體系的靜矩、線動量、平動與振動耦合慣性力、本體振動對飛行器轉(zhuǎn)動的耦合系數(shù)、舵面對飛行器轉(zhuǎn)動的附加耦合系數(shù)中含有模態(tài)坐標(biāo)的一階項，本體與舵面相對本體系的慣量張量、角動量、轉(zhuǎn)動耦合慣性力中含有模態(tài)坐標(biāo)的一階與二階項。如飛行器本體相對本體系x軸的轉(zhuǎn)動慣量Jc，xx、本體的轉(zhuǎn)動耦合慣性力Qcr的表達式分別為：

在實際的工程應(yīng)用中，為了簡化模型和減少計算量，計算本體及舵面的速度、加速度與速度變分時通常根據(jù)彈性變形為小量而將模態(tài)坐標(biāo)項忽略，只保留模態(tài)坐標(biāo)的導(dǎo)數(shù)項。簡化后本體與舵面上任一質(zhì)量單元相對慣性系的速度、加速度與速度變分分別為：

基于式（16）與式（17）推導(dǎo)的動力學(xué)模型中，慣性參數(shù)、彈性耦合系數(shù)、耦合慣性力與動量計算公式中均不含模態(tài)坐標(biāo)項，將此動力學(xué)模型定義為零次剛?cè)狁詈夏Ｐ停瑒恿抗蕉x為零次動量公式。如轉(zhuǎn)動慣量Jc，xx在此只保留式（14）中的第一項，耦合慣性力Qcr的表達式則簡化為：

在二次剛?cè)狁詈夏Ｐ偷幕A(chǔ)上，忽略模態(tài)坐標(biāo)的二階項而只保留一階項，將這樣簡化后的動力學(xué)模型定義為一次剛?cè)狁詈夏Ｐ停鄳?yīng)的動量公式定義為一次動量公式。如轉(zhuǎn)動慣量Jc，xx在此只保留式（14）中的前兩項。而對于耦合慣性力Qcr，由于式（15）僅在慣量張量Jc中含有模態(tài)坐標(biāo)的二階項，故一次與二次模型中的Qcr僅存在慣量張量保留階次的不同。在這里，考慮到飛行器彈性振動的振幅雖然較小，但振動速度較快，模態(tài)坐標(biāo)的一階導(dǎo)數(shù)的數(shù)量級高于模態(tài)坐標(biāo)，故在一次模型中未將模態(tài)坐標(biāo)與模態(tài)坐標(biāo)的一階導(dǎo)數(shù)的乘積項作為高階小量忽略，在零次模型中也未將模態(tài)坐標(biāo)的一階導(dǎo)數(shù)的一階項忽略。

3 模型剛?cè)狁詈想A次對動力學(xué)分析的影響

在本體振動的基頻f0=9.9 Hz，系統(tǒng)不受外力，舵面與本體間的相對轉(zhuǎn)角按αp=cost-1（t表示時間）規(guī)律變化的條件下，各次剛?cè)狁詈夏Ｐ蛣恿W(xué)時域響應(yīng)的部分結(jié)果對比如圖3所示。

圖3 各次剛?cè)狁詈夏Ｐ偷臅r域響應(yīng)對比Fig.3 Comparison of time domain response of rigid-flexible coupling models of each order

一次模型與二次模型的時域響應(yīng)基本吻合，而零次模型的時域響應(yīng)隨著時間的增長逐漸與一二次模型的結(jié)果出現(xiàn)偏差，絕對值較小的響應(yīng)值的偏差較為顯著。考慮到系統(tǒng)不受外力，本體與舵面間的相對轉(zhuǎn)角按正弦規(guī)律驅(qū)動，故理論上動力學(xué)時域響應(yīng)值應(yīng)在一固定的平衡位置上周期性變化，不會呈現(xiàn)出發(fā)散的變化趨勢，所以零次模型的長時間計算結(jié)果是不合理的。

在沒有外力作用且初始速度為零時，飛行器總的線動量和角動量的計算值都接近理論值0，則認(rèn)為動力學(xué)模型是合理的，計算值的絕對值越小，則認(rèn)為模型的計算精度越高。為了比較動力學(xué)模型的剛?cè)狁詈想A次對系統(tǒng)動量計算結(jié)果的影響，本節(jié)中線動量與角動量的計算均使用二次動量公式。各次剛?cè)狁詈夏Ｐ偷木€動量與角動量計算結(jié)果對比如圖4 所示，各次模型的線動量與角動量計算結(jié)果均接近0，二次模型計算結(jié)果的量級最低，可以認(rèn)為其計算精度最高。

圖4 各次剛?cè)狁詈夏Ｐ偷膭恿繉Ρ菷ig.4 Comparison of linear momentum and angular momentum of rigid-flexible coupling models of each order

零次模型的線動量與角動量計算結(jié)果比一、二次模型的高6～10 個數(shù)量級，且隨著時間的增長總體上呈近似線性地單調(diào)增大，該誤差并非由數(shù)值積分過程引入，而是模型建立時在式（16）和式（17）中忽略了部分彈性變形相關(guān)項導(dǎo)致的。僅將一次和二次模型的動量計算結(jié)果進行對比，如圖5所示。

圖5 一、二次剛?cè)狁詈夏Ｐ偷膭恿繉Ρ菷ig.5 Comparison of linear momentum and angular momentum of one-order and two-order rigid-flexible coupling models

一次模型的角動量計算結(jié)果比二次模型的高3 個數(shù)量級，且隨著時間的增長總體上呈近似線性地單調(diào)增大，而兩者線動量的計算結(jié)果較為接近。可以得出：系統(tǒng)角動量的計算精度對模型剛?cè)狁詈想A次的變化比線動量更敏感。

為了分析剛?cè)狁詈想A次對動力學(xué)分析的影響與中心柔性體振動頻率的關(guān)系，將本體有限元模型的材料彈性模量改為原來的1/100，其他條件保持不變，本體振動的基頻變?yōu)閒′0=0.99 Hz。本體基頻降低后的各次模型動力學(xué)時域響應(yīng)的部分結(jié)果對比如圖6所示。

圖6 本體基頻降低后的時域響應(yīng)對比Fig.6 Comparison of time domain response after the fundamental frequency reduction of the body

與本體振動基頻為9.9 Hz 時的動力學(xué)時域響應(yīng)結(jié)果對比，基頻降低后零次模型與高次模型的計算結(jié)果偏差顯著增大。此外，零次模型的計算結(jié)果中平動位移、歐拉角等位移量與高次模型的偏差比平動速度、角速度等速度量的偏差更明顯。可以得出：剛?cè)狁詈想A次對動力學(xué)分析的影響受中心柔性體振動頻率的影響，中心柔性體振動的基頻越低，零次模型動力學(xué)時域響應(yīng)的結(jié)果與高次模型的偏差越大，計算誤差越大。

本體基頻降低前后零次模型的線動量與角動量、一次模型的角動量計算結(jié)果對比如圖7所示。

圖7 本體基頻降低前后的動量對比Fig.7 Comparison of linear momentum and angular momentum before and after the fundamental frequency reduction of the body

本體振動基頻降低后，一次模型的角動量增大4個數(shù)量級，零次模型的線動量與角動量增大2個數(shù)量級。可以得出：中心柔性體振動頻率降低會使由于模型剛?cè)狁詈想A次降低而導(dǎo)致的線動量與角動量的計算誤差增大。

4 動量公式剛?cè)狁詈想A次對動量計算的影響

與本體振動基頻f0=9.9 Hz 時的條件保持相同，采用各次動量公式計算零次剛?cè)狁詈夏Ｐ偷木€動量與角動量的結(jié)果對比如圖8所示。

圖8 零次剛?cè)狁詈夏Ｐ偷母鞔蝿恿繉Ρ菷ig.8 Comparison of linear momentum and angular momentum of each order of zero-order rigid-flexible coupling model

各次動量公式計算所得的零次模型的線動量與角動量分別在同一個數(shù)量級。同樣可以驗證，采用二次動量公式計算一次模型的線動量與角動量，其結(jié)果與一次動量公式的計算結(jié)果在同一個數(shù)量級。可以得出：當(dāng)動力學(xué)模型的剛?cè)狁詈想A次較低時，采用較高階次的動量公式并不能使線動量與角動量的計算精度提高。

采用各次動量公式計算二次剛?cè)狁詈夏Ｐ偷木€動量與角動量的結(jié)果對比如圖9所示，僅將其中的一次和二次動量公式的計算結(jié)果進行對比如圖10所示。

圖9 二次剛?cè)狁詈夏Ｐ偷母鞔蝿恿繉Ρ菷ig.9 Comparison of linear momentum and angular momentum of each order of two-order rigid-flexible coupling model

圖10 二次剛?cè)狁詈夏Ｐ偷囊弧⒍蝿恿繉Ρ菷ig.10 Comparison of the one-order and two-order linear momentum and angular momentum of two-order rigid-flexiblecoupling model

采用零次動量公式計算二次模型的線動量與角動量，結(jié)果比二次動量公式的高8 個數(shù)量級。采用一次動量公式計算二次模型的角動量，結(jié)果比二次動量公式的高2個數(shù)量級。由于線動量的計算公式中不含模態(tài)坐標(biāo)的二階項，一次與二次線動量計算公式相同，故采用一次與二次動量公式計算的二次模型的線動量結(jié)果相同。同樣可以驗證，采用零次動量公式計算一次模型的線動量與角動量，其結(jié)果的數(shù)量級也比高次動量公式的計算結(jié)果高。可以得出：除一次與二次動量公式計算各次模型線動量的精度相同外，采用剛?cè)狁詈想A次較低的動量公式計算階次較高模型的線動量與角動量，誤差將比高次動量公式的大。

5 結(jié)論

本文研究了彈性飛行器飛行動力學(xué)建模的問題，提出了飛行器的彈性耦合系數(shù)和廣義模態(tài)質(zhì)量的修正方法，可以精確描述連接在柔性本體上的部件對飛行器動力學(xué)的影響。分析了剛?cè)狁詈想A次對計算精度的影響，發(fā)現(xiàn)隨著剛?cè)狁詈想A次的降低，動力學(xué)模型的時域響應(yīng)和動量計算結(jié)果的誤差可能會增大，且中心柔性體振動的基頻越低誤差越顯著。進一步指出了在利用動量守恒檢驗系統(tǒng)的計算精度時，動量計算公式的剛?cè)狁詈想A次應(yīng)不低于動力學(xué)模型的剛?cè)狁詈想A次。