飛輪擾動作用下衛星結構響應能量預示方法

王 坤，趙 陽，馬文來，王有懿

(哈爾濱工業大學 航天學院 航天工程系，哈爾濱 150001)

航天事業的不斷發展對航天器的指向精度及穩定性提出了更高的要求，在軌運行期間的各種干擾降低了天線和相機等有效載荷的工作性能，飛輪作為最主要的干擾源之一，嚴重影響衛星等航天器的姿態控制精度和穩定性，成為制約我國高精度航天器發展的主要因素［1］。研究航天器在飛輪作用下的動態響應特性，是抖動抑制的基礎，對我國研發高精度和高穩定度航天器具有重要意義。

研究飛輪擾動下航天器結構動態響應特性的方法主要是有限元法［2］，作為一種數值解法，為保證解的合理性，在劃分有限單元時需保證每個波長之內包含一定數目的單元［3］，這樣，網格密度和計算成本會隨著結構振動頻率的升高而增加，故實際工程中，有限元技術局限于低頻振動分析領域。

為解決航天器結構中高頻響應問題，田浩等［4］將波動方法用于衛星框式結構的抖動特性分析。波動法可對全頻域進行確定性描述，但它在高頻域對邊界條件的數學模型要求較高，也具有一定局限性。

有限元法和波動方法的一個主要缺陷在于其以位移為變量，力求得到結構所有位移和所有頻率的精確解，因而導致了響應預示結果對結構模型敏感和計算成本的增加。能量方法的出現在一定程度上解決了這一問題。Lyon［5－6］給出了基于能量方法的結構振動分析，這奠定了統計能量分析的基礎。能量方法具有計算量小、且變量對邊界條件不敏感等特點，因而在處理高頻域結構振動響應問題中得到了很大發展。統計能量分析的基本要求是高模態密度，對模態密度較低的中低頻域，要得到中低頻響應的能量預示結果，就需要突破模態密度的限制。Nefske等［7］在研究梁的振動問題時，得到了類似熱傳導方程的能量方程，該方法結合有限元技術，發展為能量有限元法，與統計能量分析一樣，用能量作為變量，計算量較小，可分析高頻域的響應。能量有限元由于不涉及模態密度，適用的頻域范圍比統計能量分析更廣，其結果可以精確到每個單元，且能得到更為精細的響應預示結果［8］，Cho［9］研究了典型耦合結構間的能量傳遞關系，使此方法更接近于實際工程應用。

在 Bouthier和 Cho基礎上，Wang［10］發展了簡化的能量有限元法，此方法在形式上與統計能量分析類似，可以方便的利用統計能量分析軟件進行響應分析。Zhang［11］則針對船舶結構，分析了流體載荷下加筋板的能量密度響應。Yan［12］進一步地推導了復合材料結構的能量密度控制方程，為能量有限元法應用于復合材料結構提供了理論依據。Lee［13］結合周期理論，將能量有限元法推廣到旋轉體結構的計算，并與試驗數據進行了對比驗證。

鑒于能量有限元法的優勢，本文將其引入航天器結構的振動響應分析。首先分析飛輪擾動模型的特性，而后基于能量密度控制方程，結合有限元技術對衛星框式結構進行抖動響應分析，為進一步研究抖動抑制提供理論依據。

1 飛輪擾動模型

文獻［1］認為轉子不平衡、飛輪結構撓性、軸承振動與摩擦以及電機的輸出轉矩等四方面是飛輪擾動主因。其中，轉子不平衡是最主要因素，它由飛輪轉子質量分布不均造成，可分為靜、動不平衡兩類，分別指飛輪轉子質心與慣性軸不重合、飛輪轉子慣性軸與旋轉軸不重合而產生的擾動。當飛輪高速旋轉時，這兩種擾動將作為高頻激勵作用于航天器結構，使其產生高頻響應，從而影響航天器的指向精度與穩定度。

根據作者前期研究［14］，通過對國內外文獻分析綜述，飛輪實驗擾動模型為:

其中:m(t)為飛輪產生的擾動力或力矩(單位:N或Nm);n為模型中包含的諧波數;Ci為第i次諧波幅值(單位:N/Hz2);Ω為飛輪的轉速(單位:Hz);hi為第i次諧波數;αi為隨機相角，假定均勻分布在［0，2π］內。

由此可見，擾動模型中的擾動力或力矩可分解成飛輪轉速的正弦波形式。其中諧波幅值Ci和諧波數hi需要利用振幅譜法和能量補償法對實驗數據進行參數辨識獲取。本文作為理論方法研究，根據飛輪擾動可以分解為正弦波激勵求和的特性，以單位正弦力作為激勵，研究頻域內結構的響應特性。在混響狀態下，飛輪擾動產生的能量密度響應遵循線性疊加原理，飛輪擾動的響應以諧波幅值和諧波數作為對應頻率響應的加權值求和即可得到。

2 能量密度控制方程

針對一控制體，根據能量平衡原理有:

其中:p為控制體內做功功率;πin為外界輸入功率;πdiss為損耗功率。

板的動力學方程為:

式中:D為板彎曲剛度，η為阻尼系數，w為板橫向位移，ρ為板密度，h為板厚度，t表示時間。Bouthier［8］將板的橫向位移表示為:

則方程可以寫為:

結構振動的解分為近場解和遠場解，近場解只在不連續處附近起作用，在遠離不連續處的位置，近場解很快衰減，只剩下遠場解的作用，在結構振動響應全局分析時，近場解可以忽略，并且考慮到阻尼系數η是個小量，方程的解可以寫為:

式中:kx，ky分別表示 x，y方向的波數。時間與空間平均能量密度為:

式中*表示求共軛。

功率為:

其中:

式中:Mxx，Mxy，Myx，Myy分別表示各方向彎矩，Qx，Qy表示各方向剪力。

將位移表達式代入，可得關系式:

代入能量平衡方程可得板的彎曲波能量密度控制方程:

其中:ω為角頻率;cg為彎曲波的群速度;e為能量密度。

由有限元理論［15］，根據加權余值法，其離散形

式為:

式中:

式中:Γ為單元邊界;n為單元邊界的單位方向矢量。NT為列向量，表示形函數向量的轉置向量。

工程中的絕大部分結構都非單獨的簡單結構，而是由簡單結構組合而成，因此需要進一步考慮簡單結構間互相耦合的情況。

由于結構不連續處的能量密度不連續，耦合結構間的能量傳遞需要通過耦合處的功率流表示。對于二維連接節點為:

式中:Q為節點功率流，e表示節點能量密度值。下標i，j表示連接單元編號，上標 m，m+1，n，n+1表示對應連接單元的節點數。KJ稱為連接矩陣，單元離散形式為:

得到考慮耦合情況的能量密度控制方程整體有限元形式為:

結構不連續處能量關系通過連接矩陣［KJ］實現。

3 算例及分析

本文設計的衛星框式承力結構如圖1所示，板的幾何尺寸均為1 m×1 m，模擬飛輪擾動的激勵作用于板1中部位置，各板的密度為ρ=2 700 kg/m3，阻尼系數 η =0.001，泊松比 ν=0.3，彈性模量 E=7.1 ×1010N/m2，板厚 h=0.001 m。

衛星框式承力結構的每個面劃分為24×24的四邊形有限元網格進行計算。

圖1 衛星框式結構圖Fig.1 satellite structure

板1是飛輪激勵直接作用的板，當載荷頻率為50 Hz、500 Hz 和 900 Hz時，對板1、板2的能量密度響應空間分布進行計算。

圖2中，板 1的中部能量密度出現峰值，此處為激勵作用位置，飛輪的擾動作為一種能量形式直接輸入到此位置，因此該位置的振動幅度最大，即對應的能量值也最高。在其它位置，由于阻尼的存在，能量密度值逐漸減小，由于能量密度為空間域內的平均值，所以整體呈現平滑的分布狀態。圖3所示的板2由于沒有激勵源輸入，整體呈現平滑狀態，并且隨著遠離激勵源，能量密度值逐漸降低。

由圖2和圖3所示各頻率下板1和板2的能量密度響應對比可以看出，在板的連接處存在能量密度值的不連續，這是由結構耦合引起的能量密度損失導致，Cho對耦合結構間能量傳遞關系進行了深入研究，用能量的透射和反射系數表征結構耦合間的能量關系，并給出了典型耦合時能量反射系數和透射系數的解析表達式［9］，但在工程中并不適用，實際中，耦合結構間的能量傳遞關系通常由實驗得到，本文利用有限元程序通過分析結構的位移響應得到。

在頻率域內，用能量有限元方法與傳統波動方法得到的結果如圖4和圖5所示。經計算，板2中點的波動解平均值為79.6 dB，EFEM 平均值為77.3 dB;板3中點波動解平均值為81.5 dB，EFEM平均值為77.2 dB;平均值相差5%左右，并且能量有限元得到的值小于波動解的值，這種現象正是能量有限元法忽略近場解的體現。由圖4和圖5可見，近場解在低頻范圍內對結構振動的響應影響較大。以圖4為例，低于200 Hz時，波動解的上下波動幅度較大，隨著頻率的升高，波長變短，近場解作用的范圍隨之變小，對整體結構的動態響應的影響也逐漸弱化，波動解在能量有限元解附近波動的幅度逐漸減小。

由圖4和圖5顯示，傳統波動方法得到的能量密度響應值呈現出密集的模態，對頻率的變化敏感度很高，頻率的微小改變會引起振動預示結果急劇變化，而利用能量有限元技術得到的響應曲線呈平滑變化狀態，它預示的是在頻域范圍內平均的振動能量狀態。

4 結論

本文基于飛輪擾動模型特性，針對衛星框式結構，采用能量有限元方法進行了結構響應預示研究，得出了以下結論與啟示:

(1)飛輪擾動力在飛輪作用處的響應最大，響應值會隨著與激勵源距離的增大而減小，所以針對飛輪擾動的隔振在飛輪附近可以取得最好效果。

(2)能量有限元法預示的是衛星結構響應的平均效果，并且低于傳統方法的預示值，相對于傳統以位移為變量的方法，更適合中高頻響應預示。

(3)能量有限元法在預示飛輪擾動時對頻率變化不敏感。而實際工程中并不存在激勵的精確模型，因此，能量有限元法的預示效果更具有參考價值。

飛輪是引起航天器擾動的主要因素但不是唯一因素，而且航天器的結構十分復雜，為更好的解決航天器擾動引起的指向精度等問題，需要綜合考慮各種擾動的特征并結合更為精細的航天器模型進一步分析。

［1］王全武，虎 剛.飛輪擾動原因與測量技術現狀［J］.空間科學學報，2009，29(1):39 －44.

［2］白爭鋒，趙 陽，馬文來，等.反作用輪擾動對航天器結構動態特性的影響分析［J］.宇航學報，2009，30(5):2073－2079.

［3］韓增堯.有限元頻響分析中網格劃分技術研究［J］.強度與環境，2003，30(3):18 －22.

［4］田 浩，王有懿，劉春川.衛星框式結構擾動動態特性分析的行波方法［J］.哈爾濱工業大學學報，2010，42(5):696－699.

［5］Lyon R H，Maidanik G.Power flow between linearly coupled oscillators［J］.J Acoust.Soc.Am.，1962，34(5):623－639.

［6］Lyon R H，Scharton T D.Vibrational energy transmission in a three element structure［J］.J Sound and Vibration，1965，38:253－261.

［7］Nefske D J，Sung S H.Power flow finite element analysis of dynamic system:basic theory and application to beams［J］.Transactions of the ASME，1989，111(1):94－100.

［8］Bouthier O M，Bernhard R J.Simple models of the energetics of transversely vibrating plates［J］.J Sound and vibration，1995，182(1):149－164.

［9］Cho P E.Energy flow analysis of coupled structures［D］.West Lafayette:Purdue University，1992.

［10］Wang S.High frequency energy flow analysis methods:numerical implement-tation， applications， and verification［D］.West Lafayette:Purdue University，2000.

［11］Zhang W G.Energy finite element method for vibration analysisofstiffened platesunderfluid loading［D］.Michigan:University of Michigan，2003.

［12］Yan X Y.Energy finite element analysis development for high frequency vibration analysis of composite structures［D］.Michigan:University of Michigan，2008.

［13］Lee S M.Energy finite element method for high frequency vibration analysis of composite rotorcraft structure［D］.Michigan:University of Michigan，2003.

［14］孫 杰，趙 陽，田 浩.改善航天器反作用輪擾動實驗模型參數識別方法［J］.空間科學學報，2006，26(1):70－74.

［15］Zienkiewicz O C.The finite element method［M］.New York:McGraw-Hill，1977.