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

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

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

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

為解決航天器結(jié)構(gòu)中高頻響應(yīng)問題，田浩等［4］將波動方法用于衛(wèi)星框式結(jié)構(gòu)的抖動特性分析。波動法可對全頻域進(jìn)行確定性描述，但它在高頻域?qū)吔鐥l件的數(shù)學(xué)模型要求較高，也具有一定局限性。

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

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

鑒于能量有限元法的優(yōu)勢，本文將其引入航天器結(jié)構(gòu)的振動響應(yīng)分析。首先分析飛輪擾動模型的特性，而后基于能量密度控制方程，結(jié)合有限元技術(shù)對衛(wèi)星框式結(jié)構(gòu)進(jìn)行抖動響應(yīng)分析，為進(jìn)一步研究抖動抑制提供理論依據(jù)。

1 飛輪擾動模型

文獻(xiàn)［1］認(rèn)為轉(zhuǎn)子不平衡、飛輪結(jié)構(gòu)撓性、軸承振動與摩擦以及電機(jī)的輸出轉(zhuǎn)矩等四方面是飛輪擾動主因。其中，轉(zhuǎn)子不平衡是最主要因素，它由飛輪轉(zhuǎn)子質(zhì)量分布不均造成，可分為靜、動不平衡兩類，分別指飛輪轉(zhuǎn)子質(zhì)心與慣性軸不重合、飛輪轉(zhuǎn)子慣性軸與旋轉(zhuǎn)軸不重合而產(chǎn)生的擾動。當(dāng)飛輪高速旋轉(zhuǎn)時，這兩種擾動將作為高頻激勵作用于航天器結(jié)構(gòu)，使其產(chǎn)生高頻響應(yīng)，從而影響航天器的指向精度與穩(wěn)定度。

根據(jù)作者前期研究［14］，通過對國內(nèi)外文獻(xiàn)分析綜述，飛輪實(shí)驗擾動模型為:

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

由此可見，擾動模型中的擾動力或力矩可分解成飛輪轉(zhuǎn)速的正弦波形式。其中諧波幅值Ci和諧波數(shù)hi需要利用振幅譜法和能量補(bǔ)償法對實(shí)驗數(shù)據(jù)進(jìn)行參數(shù)辨識獲取。本文作為理論方法研究，根據(jù)飛輪擾動可以分解為正弦波激勵求和的特性，以單位正弦力作為激勵，研究頻域內(nèi)結(jié)構(gòu)的響應(yīng)特性。在混響狀態(tài)下，飛輪擾動產(chǎn)生的能量密度響應(yīng)遵循線性疊加原理，飛輪擾動的響應(yīng)以諧波幅值和諧波數(shù)作為對應(yīng)頻率響應(yīng)的加權(quán)值求和即可得到。

2 能量密度控制方程

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

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

板的動力學(xué)方程為:

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

則方程可以寫為:

結(jié)構(gòu)振動的解分為近場解和遠(yuǎn)場解，近場解只在不連續(xù)處附近起作用，在遠(yuǎn)離不連續(xù)處的位置，近場解很快衰減，只剩下遠(yuǎn)場解的作用，在結(jié)構(gòu)振動響應(yīng)全局分析時，近場解可以忽略，并且考慮到阻尼系數(shù)η是個小量，方程的解可以寫為:

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

式中*表示求共軛。

功率為:

其中:

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

將位移表達(dá)式代入，可得關(guān)系式:

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

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

由有限元理論［15］，根據(jù)加權(quán)余值法，其離散形

式為:

式中:

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

工程中的絕大部分結(jié)構(gòu)都非單獨(dú)的簡單結(jié)構(gòu)，而是由簡單結(jié)構(gòu)組合而成，因此需要進(jìn)一步考慮簡單結(jié)構(gòu)間互相耦合的情況。

由于結(jié)構(gòu)不連續(xù)處的能量密度不連續(xù)，耦合結(jié)構(gòu)間的能量傳遞需要通過耦合處的功率流表示。對于二維連接節(jié)點(diǎn)為:

式中:Q為節(jié)點(diǎn)功率流，e表示節(jié)點(diǎn)能量密度值。下標(biāo)i，j表示連接單元編號，上標(biāo) m，m+1，n，n+1表示對應(yīng)連接單元的節(jié)點(diǎn)數(shù)。KJ稱為連接矩陣，單元離散形式為:

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

結(jié)構(gòu)不連續(xù)處能量關(guān)系通過連接矩陣［KJ］實(shí)現(xiàn)。

3 算例及分析

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

衛(wèi)星框式承力結(jié)構(gòu)的每個面劃分為24×24的四邊形有限元網(wǎng)格進(jìn)行計算。

圖1 衛(wèi)星框式結(jié)構(gòu)圖Fig.1 satellite structure

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

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

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

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

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

4 結(jié)論

本文基于飛輪擾動模型特性，針對衛(wèi)星框式結(jié)構(gòu)，采用能量有限元方法進(jìn)行了結(jié)構(gòu)響應(yīng)預(yù)示研究，得出了以下結(jié)論與啟示:

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

(2)能量有限元法預(yù)示的是衛(wèi)星結(jié)構(gòu)響應(yīng)的平均效果，并且低于傳統(tǒng)方法的預(yù)示值，相對于傳統(tǒng)以位移為變量的方法，更適合中高頻響應(yīng)預(yù)示。

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

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

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

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

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

［4］田 浩，王有懿，劉春川.衛(wèi)星框式結(jié)構(gòu)擾動動態(tài)特性分析的行波方法［J］.哈爾濱工業(yè)大學(xué)學(xué)報，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］孫 杰，趙 陽，田 浩.改善航天器反作用輪擾動實(shí)驗?zāi)Ｐ蛥?shù)識別方法［J］.空間科學(xué)學(xué)報，2006，26(1):70－74.

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