彈性飛行器飛行動力學建模與剛柔耦合階次分析

安 陽，王天舒

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

0 引言

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

在目前的彈性飛行器動力學建模領域，對描述彈性體變形運動的變形場階次和模態截斷階次等問題已有較多研究［18，28-29］，而對在模態階數與彈性變形場階次選定后，在動力學模型推導過程中產生的各階彈性變形相關項取舍問題的研究較少。在傳統的中心剛體-柔性附件模型中，中心剛體的質量、轉動慣量等在飛行器整體中占有絕對優勢，柔性附件的彈性振動對飛行器整體動力學的影響較小，因此通常忽略動力學模型推導過程中產生的關于彈性變形的高階項，而近似后的模型仍能保有滿足工程要求的精度。但當飛行器主體呈現柔性特性時，這種高階項是否還能忽略及其對整體動力學的影響尚無定論。此外，在動力學模型建立后，通常選用無外力作用時系統的線動量和角動量作為守恒量，檢驗模型的正確性和對比模型的計算精度。系統動量表達式中彈性變形相關項的保留階次也可能對其計算結果的準確性產生影響。研究彈性變形相關項的保留程度對飛行器動力學及動量計算的影響，對動力學模型的降階與驗證具有重要的指導意義。本文首先基于虛功率原理與混合坐標法建立彈性飛行器的飛行動力學模型，并定義剛柔耦合階次以描述動力學模型推導過程中產生的彈性變形相關項的保留程度，然后以由柔性中心體和剛性附件組成的剛柔耦合系統為例，對比分析動力學模型及動量公式的剛柔耦合階次對動力學響應與守恒量的影響。

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

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

圖1 飛行器系統構型與坐標系定義示意圖Fig.1 Configuration of the flight vehicle system and coordinate system definition diagram

虛功率原理的表達式為：

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

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

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

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

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

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

飛行器平動動力學方程為：

式中：m為飛行器的總質量；S為飛行器相對本體系的總靜矩；ωc為本體相對慣性系的轉動角速度，即飛行器的姿態角速度；F為飛行器主動力的合力；Qt為飛行器平動耦合慣性力的合力；耦合系數陣Ct的表達式為：

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

飛行器姿態動力學方程為：

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

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

飛行器振動動力學方程為：

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

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

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

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

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

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

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

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

2 剛柔耦合階次的定義

上文推導的動力學模型是未忽略各階彈性變形相關項的完整模型，其動力學方程與系統角動量表達式中含有模態坐標的一階與二階項，將此動力學模型定義為二次剛柔耦合模型，動量公式定義為二次動量公式。在二次模型中，本體與舵面相對本體系的靜矩、線動量、平動與振動耦合慣性力、本體振動對飛行器轉動的耦合系數、舵面對飛行器轉動的附加耦合系數中含有模態坐標的一階項，本體與舵面相對本體系的慣量張量、角動量、轉動耦合慣性力中含有模態坐標的一階與二階項。如飛行器本體相對本體系x軸的轉動慣量Jc，xx、本體的轉動耦合慣性力Qcr的表達式分別為：

在實際的工程應用中，為了簡化模型和減少計算量，計算本體及舵面的速度、加速度與速度變分時通常根據彈性變形為小量而將模態坐標項忽略，只保留模態坐標的導數項。簡化后本體與舵面上任一質量單元相對慣性系的速度、加速度與速度變分分別為：

基于式（16）與式（17）推導的動力學模型中，慣性參數、彈性耦合系數、耦合慣性力與動量計算公式中均不含模態坐標項，將此動力學模型定義為零次剛柔耦合模型，動量公式定義為零次動量公式。如轉動慣量Jc，xx在此只保留式（14）中的第一項，耦合慣性力Qcr的表達式則簡化為：

在二次剛柔耦合模型的基礎上，忽略模態坐標的二階項而只保留一階項，將這樣簡化后的動力學模型定義為一次剛柔耦合模型，相應的動量公式定義為一次動量公式。如轉動慣量Jc，xx在此只保留式（14）中的前兩項。而對于耦合慣性力Qcr，由于式（15）僅在慣量張量Jc中含有模態坐標的二階項，故一次與二次模型中的Qcr僅存在慣量張量保留階次的不同。在這里，考慮到飛行器彈性振動的振幅雖然較小，但振動速度較快，模態坐標的一階導數的數量級高于模態坐標，故在一次模型中未將模態坐標與模態坐標的一階導數的乘積項作為高階小量忽略，在零次模型中也未將模態坐標的一階導數的一階項忽略。

3 模型剛柔耦合階次對動力學分析的影響

在本體振動的基頻f0=9.9 Hz，系統不受外力，舵面與本體間的相對轉角按αp=cost-1（t表示時間）規律變化的條件下，各次剛柔耦合模型動力學時域響應的部分結果對比如圖3所示。

圖3 各次剛柔耦合模型的時域響應對比Fig.3 Comparison of time domain response of rigid-flexible coupling models of each order

一次模型與二次模型的時域響應基本吻合，而零次模型的時域響應隨著時間的增長逐漸與一二次模型的結果出現偏差，絕對值較小的響應值的偏差較為顯著?？紤]到系統不受外力，本體與舵面間的相對轉角按正弦規律驅動，故理論上動力學時域響應值應在一固定的平衡位置上周期性變化，不會呈現出發散的變化趨勢，所以零次模型的長時間計算結果是不合理的。

在沒有外力作用且初始速度為零時，飛行器總的線動量和角動量的計算值都接近理論值0，則認為動力學模型是合理的，計算值的絕對值越小，則認為模型的計算精度越高。為了比較動力學模型的剛柔耦合階次對系統動量計算結果的影響，本節中線動量與角動量的計算均使用二次動量公式。各次剛柔耦合模型的線動量與角動量計算結果對比如圖4 所示，各次模型的線動量與角動量計算結果均接近0，二次模型計算結果的量級最低，可以認為其計算精度最高。

圖4 各次剛柔耦合模型的動量對比Fig.4 Comparison of linear momentum and angular momentum of rigid-flexible coupling models of each order

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

圖5 一、二次剛柔耦合模型的動量對比Fig.5 Comparison of linear momentum and angular momentum of one-order and two-order rigid-flexible coupling models

一次模型的角動量計算結果比二次模型的高3 個數量級，且隨著時間的增長總體上呈近似線性地單調增大，而兩者線動量的計算結果較為接近?？梢缘贸觯合到y角動量的計算精度對模型剛柔耦合階次的變化比線動量更敏感。

為了分析剛柔耦合階次對動力學分析的影響與中心柔性體振動頻率的關系，將本體有限元模型的材料彈性模量改為原來的1/100，其他條件保持不變，本體振動的基頻變為f′0=0.99 Hz。本體基頻降低后的各次模型動力學時域響應的部分結果對比如圖6所示。

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

與本體振動基頻為9.9 Hz 時的動力學時域響應結果對比，基頻降低后零次模型與高次模型的計算結果偏差顯著增大。此外，零次模型的計算結果中平動位移、歐拉角等位移量與高次模型的偏差比平動速度、角速度等速度量的偏差更明顯。可以得出：剛柔耦合階次對動力學分析的影響受中心柔性體振動頻率的影響，中心柔性體振動的基頻越低，零次模型動力學時域響應的結果與高次模型的偏差越大，計算誤差越大。

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

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

本體振動基頻降低后，一次模型的角動量增大4個數量級，零次模型的線動量與角動量增大2個數量級?？梢缘贸觯褐行娜嵝泽w振動頻率降低會使由于模型剛柔耦合階次降低而導致的線動量與角動量的計算誤差增大。

4 動量公式剛柔耦合階次對動量計算的影響

與本體振動基頻f0=9.9 Hz 時的條件保持相同，采用各次動量公式計算零次剛柔耦合模型的線動量與角動量的結果對比如圖8所示。

圖8 零次剛柔耦合模型的各次動量對比Fig.8 Comparison of linear momentum and angular momentum of each order of zero-order rigid-flexible coupling model

各次動量公式計算所得的零次模型的線動量與角動量分別在同一個數量級。同樣可以驗證，采用二次動量公式計算一次模型的線動量與角動量，其結果與一次動量公式的計算結果在同一個數量級。可以得出：當動力學模型的剛柔耦合階次較低時，采用較高階次的動量公式并不能使線動量與角動量的計算精度提高。

采用各次動量公式計算二次剛柔耦合模型的線動量與角動量的結果對比如圖9所示，僅將其中的一次和二次動量公式的計算結果進行對比如圖10所示。

圖9 二次剛柔耦合模型的各次動量對比Fig.9 Comparison of linear momentum and angular momentum of each order of two-order rigid-flexible coupling model

圖10 二次剛柔耦合模型的一、二次動量對比Fig.10 Comparison of the one-order and two-order linear momentum and angular momentum of two-order rigid-flexiblecoupling model

采用零次動量公式計算二次模型的線動量與角動量，結果比二次動量公式的高8 個數量級。采用一次動量公式計算二次模型的角動量，結果比二次動量公式的高2個數量級。由于線動量的計算公式中不含模態坐標的二階項，一次與二次線動量計算公式相同，故采用一次與二次動量公式計算的二次模型的線動量結果相同。同樣可以驗證，采用零次動量公式計算一次模型的線動量與角動量，其結果的數量級也比高次動量公式的計算結果高。可以得出：除一次與二次動量公式計算各次模型線動量的精度相同外，采用剛柔耦合階次較低的動量公式計算階次較高模型的線動量與角動量，誤差將比高次動量公式的大。

5 結論

本文研究了彈性飛行器飛行動力學建模的問題，提出了飛行器的彈性耦合系數和廣義模態質量的修正方法，可以精確描述連接在柔性本體上的部件對飛行器動力學的影響。分析了剛柔耦合階次對計算精度的影響，發現隨著剛柔耦合階次的降低，動力學模型的時域響應和動量計算結果的誤差可能會增大，且中心柔性體振動的基頻越低誤差越顯著。進一步指出了在利用動量守恒檢驗系統的計算精度時，動量計算公式的剛柔耦合階次應不低于動力學模型的剛柔耦合階次。