基于ANCF的撓性太陽翼動力學建模及模態辨識研究

王慕昊, 舒通通, 魏鴻超, 趙亞濤, 夏 斌*

1. 中國星網網絡創新研究院有限公司,北京 100029

2. 哈爾濱工業大學,哈爾濱 150001

0 引 言

太空飛行任務的多樣化和規模化趨勢,正推動著當前太空航天器結構朝著更大型和更具柔性化的方向發展[1-4].太陽能帆板、大型桁架、天線和大型機械臂等都是航天器實際應用中常見的典型大撓性附件[5-8],獲取其動力學性能參數為航天器控制系統提供輸入變得至關重要.10 m級的薄膜結構尚能進行地面真空試驗,但目前大型薄膜結構展開尺度達到100 m,例如在地球靜止軌道上使用2.4 GHz作為工作頻率運行的紅外通訊衛星,其薄膜結構直徑需達到135 m[9].超大型太陽翼展開后的巨大尺寸成為限制其進行真空試驗的首要因素,而不考慮空氣阻尼測試獲取的薄膜結構動力學性能參數與真實在軌參數具有很大差異,因此,進行精細化的動力學建模和在軌模態參數辨識研究變得至關重要.通過仿真模擬在軌航天器工況來獲得更真實的模態參數,以評估薄膜結構的動力學性能,成為解決這一問題的有效途徑之一[10].一些學者已經對含預應力可展薄膜結構進行了相關的動力學分析研究.胡宇[11]研究了薄膜陣面結構剛度影響因素,試驗結果表明膜面振型與數值模擬結果基本吻合,但對應頻率有較大差異,原因可能為試驗條件未達到理想狀態.邱慧等[12]在研究中對含薄膜和彈性索張拉系統的航天器平面結構的動力學特性進行了分析,并通過試驗結果驗證了濕模態仿真分析的準確性,從而推理出干模態仿真分析的一致性,并未進一步作干模態的試驗驗證.國際上的相關機構已經進行了與模態在軌辨識相關的試驗驗證工作,如美國針對直徑8m的膜反射陣列天線的振動分析,利用雙變量參數(two-variable-parameter,2-VP)膜模型和分布式傳遞函數法(distributed transfer function method,DTFM)預測天線面內應力的分布和振動模式,但忽略了空氣阻尼對膜固有頻率的影響[13].ZHANG等[14]基于溫度-結構預應力導入方法,采用特征系統實現算法(eigensystem realization algorithm,ERA)對薄膜結構平面天線進行模態辨識,前兩階辨識準確率在5%以內.然而,上述相關研究多集中在試驗或有限元模型,對于利用絕對節點坐標法進行含預緊張力的薄膜結構平面天線動力學建模與分析的報道較少.

本文基于應變-預應力導入方法,采用ANCF索梁模型和ANCF薄膜單元對某低軌衛星的太陽翼柔性附件開展動力學建模及分析,綜合考慮了太陽光壓力、氣動阻力、重力梯度力矩和地磁力矩在內的復雜空間環境干擾,以高斯白噪聲曲線加速度作為激勵,施加于撓性結構連接處并獲取薄膜結構的均布響應,基于輸入輸出數據,采用復模態指示函數法辨識薄膜復合結構的動力學參數.通過辨識結果與有限元模態分析結果對比分析,研究結果表明CMIF方法可有效地識別出薄膜復合結構的低階固有模態,為該技術的工程實施提供了堅實的理論研究基礎.

1 問題描述

1.1 結構組成

低軌衛星的撓性部件主要由太陽翼邊框、太陽翼薄膜以及張力撐桿組成.星本體通過舉升機構連接大面積太陽翼,太陽翼薄膜兩側對稱布置張拉機構,太陽翼薄膜靠近剛性板的一端與剛性板布置有3處壓板.太陽翼邊框展開后,張力撐桿承受張拉機構的反作用力,有效降低邊框的預應力水平.單側太陽翼部件的具體約束形式見圖1.撓性太陽翼的長寬比達5:1,衛星姿態機動過程中極易激發耦合振動,因此獲取太陽翼部件動力學模態參數作為航天器控制系統的設計約束是必不可少的環節.

圖1 單側太陽翼部件約束示意圖Fig.1 Constraints on one-sided solar wing components

1.2 材料屬性

太陽翼邊框采用碳纖維/環氧復合材料,截面形狀及鋪層設計參考文獻[14].張力撐桿為鋁合金材料,截面為管型.太陽翼薄膜采用聚酰亞胺樹脂材料,厚度為0.4 mm.拉索為 2 mm圓截面凱夫拉纖維,設計張緊力20 N.

2 動力學模型及分析結果

2.1 ANCF索梁單元與殼單元

ANCF全參數梁單元可以描述梁截面的變形,適合于大截面梁的建模.但對于橫截面很小、柔性較大的細長梁來說,截面變形不是關注的重點,從而可忽略沿截面方向的斜率矢量.這種節點坐標中忽略一個或幾個方向斜率矢量的單元稱作縮減單元,或稱作梯度不足單元.

BERZERI等[15]首先提出了這種類型的平面縮減梁單元,GERSTMAYR等[16]將該單元推廣為空間縮減梁單元(又稱為索單元).利用ANCF索梁單元建立拉索的動力學模型,如圖2所示.

圖2 空間縮減梁單元Fig.2 ANCF beam element

在絕對坐標系中,該單元中任意一點的位置可定義為

(1)

式中,S為3×12的形函數矩陣.單元的廣義坐標可定義為

(2)

顯然,該單元不能描述梁繞自身軸線的轉動.對于ANCF縮減曲梁單元,其單元動能可寫為

(3)

式中,單元的長度為L,ρ和A分別為索梁單元的密度和橫截面積,M為單元的質量矩陣,矩陣中各項均為常數.

基于歐拉-伯努利梁理論,對于ANCF縮減曲梁單元,其應變能分為兩部分:1)由中線的軸向變形引起;2)由彎曲變形引起,具體表達式為

(4)

式中,E為彈性模量,J為截面慣性矩,aε為軸向應變,aκ和aκ0分別為曲梁當前狀態與初始狀態的曲率.

如圖3所示,基于參考文獻[17]建立雙線性剪切變形ANCF殼單元.

圖3 雙線性剪切變形ANCF殼單元Fig.3 ANCF shell element

薄膜單元通過網格劃分,離散成一定數量的4節點單元,定義薄膜單元上任意物質點x的全局坐標br為

(5)

(6)

式中,向量bep、beg代表了單元上與節點全局位置和橫向梯度相關的矢量.對于單元b上節點k有bkep=bkrk,bkeg=?bkr/?bz.則ANCF全局坐標表示為

br(bx,by,bz)=bS(bx,by,bz)be

(7)

形函數矩陣和節點全局坐標向量bS和be的表達式分別為

(8)

2.2 平衡方程

通過拉格朗日方程導出的薄膜結構系統動力學方程,具體形式為

(9)

式中,q為廣義坐標,T為總動能,U為總彈性能,C(q,t)=0為約束方程,λ為約束方程對應的拉式乘子,Qf為外力的廣義力.

拉索預緊力可通過定義ANCF索梁單元的單元長度L來實現.結構由于單元長度的變化,其應變為

ε=(L-L0)/L0

(10)

式中,L0為無應力狀態下索梁單元的初始長度.

因此,給定單元參考長度L0,可以實現初始預應力的導入,同時在支撐邊框、薄膜結構和張力撐桿等多重作用下維持初始狀態的平衡.

2.3 邊界條件

如圖4所示,衛星包含4根太陽翼邊框和2塊太陽翼柔性基板,太陽翼邊框編號為①～④,太陽翼柔性基板編號分別為A和B.

圖4 柔性體布局Fig.4 Flex body layout

以+Y側大撓性太陽翼裝配體為例,簡述柔性體建模過程.如圖5所示,根據太陽翼幾何包絡尺寸抽取節點,其中每根太陽翼邊框選取12個節點,共組成11個ANCF索梁單元.太陽翼A抽取11×3個節點,每4個節點按照順時針順序組成一個ANCF薄膜單元,共計20個單元.拉索處建立ANCF索梁單元.張力撐桿、拉索和邊框等可展薄膜結構各個零件之間的連接采用耦合自由度的方式,以實現位移邊界的協調.

2.4 分析結果

由表1材料特性及式(10)可計算出單元參考長度,該預載荷可實現膜面張緊.編號3_2～3_10的張力機構對應圖6中1～9號拉索穩態拉力.基于薄膜結構的固支約束邊界和張拉力,利用有限元軟件可以計算出薄膜結構的主要頻率和振型(前6階見圖7).

表1 材料性能參數(常溫)Tab.1 Material performance parameters

圖6 穩態拉索拉力Fig.6 Steady cable force

圖7 薄膜結構前6階頻率及振型(有限元)Fig.7 Membrane structure vibration modes(FEM)

3 參數辨識算法

本文基于CMIF[18]進行模態參數辨識.

步驟1.利用仿真實測數據繪制的頻響函數矩陣H(jω)∈CNo×Ni(No、Ni分別表示輸出、輸入數目)在固有頻率附近進行奇異值分解,有如下公式:

[H]=[U][Σ][V]H

(11)

式中:U為左奇異矩陣,對應于模態振型向量矩陣;Σ為對角線型奇異值矩陣;V為右奇異矩陣,對應于模態參與因子向量矩陣.

步驟2.構造增強頻響函數,其表達式為

Hen(jω)=uHH(jω)v

(12)

式中,u∈CNo×1、v∈CNi×1分別是U、V的第一列向量.

步驟3.對式(12)得到的增強頻響函數在每階固有頻率附近進行單自由度系統擬合

(13)

式中,λ為極點,未知數b1、b0和λ可通過最小二乘擬合計算得到

Xθ=Y

(14)

式中,X、Y和θ的表達式為

(15)

式中,Nf為給定頻段內所包含的譜線數.

步驟4.計算出b1、b0和λ后,可計算固有頻率f與阻尼比ξ

(16)

4 數值仿真及分析

本文基于多體系統動力學軟件(MBDyn)[19]建立大撓性太陽翼多體模型,見圖8.

圖8 太陽翼仿真示意Fig.8 Solar wing simulation

利用第2.3節建立的薄膜結構平面太陽翼模型,綜合考慮太陽光壓力、氣動阻力、重力梯度力矩和地磁力矩在內的復雜空間環境干擾,在可展結構連接點施加高斯白噪聲加速度激勵作為輸入數據,通過瞬態非線性動力學分析,在薄膜陣面和邊框上取圖5所示共計54個點的位移響應作為輸出數據.由于不同模態的貢獻沿頻率軸變化,因此在特定頻率下,兩種模態的貢獻可以大致相等.在此頻率下,這兩條特征值曲線相互交叉.由于頻率分辨率和CMIF繪制方式的限制,較低的特征值曲線出現峰值,而較高的特征值曲線出現下降,這樣就會導致虛假模態的出現.本文采用模態相位共線性法(modal phase collinearity,MPC)根據辨識得到的各階模態矢量實部和虛部之間的關系來判斷模態的真假.MPC定義如下[20]:

(17)

(18)

式中,第j個模態的振型為φj,λ1、λ2分別表示Scov特征值的最大值和最小值.如果僅一個特征值與0不同,則在復平面中呈現完美的共線性,MPC會等于1,如果兩個特征值相等,則不存在共線性,MPC會等于0.

瞬態動力學仿真過程中步長取0.001 s,仿真時間為100 s,延長原始數據至300 s并作加窗處理.如圖9所示,在2 Hz內6個頻率下,對應的奇異值存在極大值.如表2所示,由仿真數據辨識出的前3階模態的MPC值接近1,可排除虛假模態.但是第4階模態的MPC值顯著降低,消除環境干擾后進行模態辨識,第4階模態的MPC顯著趨向于1,可見環境干擾會影響模態之間的相位關系,從而導致較低的共線性矩陣值.圖10為辨識結果,可以看出辨識結果與有限元分析的振型吻合較好.

表2 模態MPC值對比Tab.2 Modal MPC values comparison

表3 固有頻率辨識結果Tab.3 Natural frequency identification results

圖9 CMIF峰Fig.9 CMIF peaks

圖10 模態辨識結果Fig.10 Modal identification results

5 結 論

結合可展薄膜結構平面太陽翼纖薄和輕柔的特性,本文建立了考慮拉索預應力的大型空間可展結構的動力學模型,利用CMIF對復雜空間環境干擾下的太陽翼模型進行模擬在軌模態辨識.仿真結果表明,CMIF方法可以有效辨識出薄膜平面結構的低階模態,辨識結果與有限元分析的比對誤差小于4%,可為后續開展工程化實施和在軌振動控制提供有效指導.

但是,在輸出數據采集過程中沒有考慮信號噪聲和結構非線性等因素的影響,體現為頻率辨識結果可能與實際傳感器數據分析結果存在較大差距.因此,下一步的工作目標是探究CMIF識別算法中模態參數不確定性的計算方法.