朱衛紅,馬興瑞,韓增堯,鄒元杰
(中國空間技術研究院,北京100094)
準確預示航天器的力學環境是指導航天器設計和地面試驗方案與條件制定的重要依據,但是航天器在發射過程中經歷的力學環境頻率寬,一般在10 Hz~10 kHz,很難采用單一的分析方法實現精確的力學環境預示[1]。在低頻段航天器結構和聲腔的模態稀疏,有限元和邊界元等基于單元離散技術的方法最為常用;而在高頻段,結構和聲腔的模態密集且隨機特性影響突出,因此,統計能量分析和能量有限元分析等都是比較有效的預示手段。但是在低頻和高頻中間存在一個過渡頻段(一般對于航天器系統級模型大約在300 Hz~1 kHz),在該頻段航天器各子結構或部組件之間的動力學特性差異大,此時傳統低頻方法和高頻方法都無法準確地進行預示,這就是航天器中頻段的力學環境預示問題。目前中頻段的響應預示是航天器結構力學環境預示的難點,也是國內外研究的熱點[2],因此展開航天器中頻段的預示方法研究具有重要的理論價值和工程應用價值。
為了解決中頻段的力學環境預示問題,2005年Shorter等[3]提出了一種混合FE-SEA中頻預示方法(Hybrid Finite Element-Statistic Energy Method)。該方法首先將系統根據動力學特性劃分為確定性子系統和隨機子系統,然后確定性子系統采用有限元方法建模,隨機子系統采用SEA方法建模,最后耦合求解得到系統不同子系統的響應。該方法提出后受到了廣泛關注,目前國外已將其廣泛應用于航空航天領域[4-6]。我國航天工業部門也在“十一五”初期啟動了該方法的理論和應用研究,目前已完成了該方法的基本理論推導和部分仿真分析及試驗驗證[7-9],但是結合目前我國航天工程中頻預示的實際需求,仍然需要針對該方法從建模理論、工程應用和試驗驗證等方面展開更深入的研究[2]。
對于大型復雜航天器,不同部組件的連接形式通常復雜多樣,總體可歸納為點、線、面三種基本連接形式,其中線連接是航天器結構中一種比較普遍的連接形式,如星體壁板的交界處,又如衛星與火箭連接的星箭界面等。如果線連接的兩類結構剛度差異大,這就是典型的中頻預示問題。根據混合FESEA理論,需采用有限元方法對剛性較大的部件建模,而柔性組件則采用統計能量方法建模,然后通過混合線連接耦合進行求解分析。目前國內外公開可查閱的混合線連接建模的資料非常少,因此本文針對航天器結構中頻預示的混合線連接建模理論展開研究,最后以典型結構為例驗證了本文方法的有效性。
混合FE-SEA基本理論認為統計能量子系統的能量由直接場和混響場兩部分能量組成。直接場能量為通過連接輸入到子系統后未經過任何邊界反射的能量,混響場能量為經過多次反射后的能量,圖1為線連接結構的直接場和混響場示意圖。
假設統計能量子系統的混合線連接m的位移自由度為q(m),則混合線連接的動力學響應可表示為:


圖1 線連接的直接場和混響場Fig.1 Direct field and diffuse reverberant field for line junction


上標-H為共軛轉置求逆運算符,力互譜矩陣Sff包括外載荷互譜和混響受擋力互譜:


式中:Im(·)為取虛部,E為統計能量子系統的混響能量,ω為分析頻率,n為模態密度。方程描述了統計能量子系統的能量響應與其在混合線連接m處產生的反力之間的關系。同時基于該關系可建立耦合后的統計能量模型的功率流平衡方程[11]:


式中:式中:Ddir,j為第j個統計能量子系統所有混合線連接直接場動剛度矩陣之和,ηd,j為有限元模型對統計能量模型j的等效損耗因子,ηje為SEA子系統j和SEA子系統e間的等效耦合損耗因子,Pjin為施加在統計載荷輸入到SEA子系統j的等效輸入功率,nj為SEA子系統j的模態密度。分析時首先求解耦合后的功率流平衡方程(6),求解得到統計能量模型的能量響應,然后由式(5)得混合連接邊界處的混響受擋力,代入方程(3)求得有限元模型的響應。
從上述求解流程可以看出,混合FE-SEA分析關鍵在于混合連接的建模,即混合線連接的直接場動剛度矩陣。下面將針對如何建立混合線連接的直接場動剛度矩陣展開研究。
航天器結構以具有大面質比的板殼結構為主,如星體壁板結構和太陽翼等,因此本文重點研究板殼結構的混合線連接模型。假設板j通過線連接與其他結構相連,定義其局部坐標系為oxjyjzj,如圖2所示。線連接沿xj軸,其運動可由沿三個方向的平動[uj,vj,wj]T與繞 xj軸的轉角 θj描述,這里主要考慮面外運動[wj,θj]T,面內運動的建模將不再贅述。根據各向同性薄板理論,板的橫向自由運動的控制方程為:


圖2 線連接的局部坐標系Fig.2 The local coordinate of line junction
波從線連接處向板結構內部傳遞能量,波的形式可寫為exp(-i kxj+μyj+iωt),負號表示波沿xj軸正方向傳播,k為該方向波數,則橫向位移可寫為:

式中:α為橫向位移的幅值,將式(8)代入式(7)則可得出:

式中:板的彎曲波數 kB=(ρjhjω2/Dj)1/4。由方程可求出4個μ值,由于yj方向的波只能沿正方向傳播,因此μ只能取負值。當k≥kB時,μ取兩個負實數,這表明此時yj方向的波為近場波,不傳遞能量;當k<kB時,μ取一個負實數和一個負虛數,此時yj方向的波為一個近場波和一個遠場波。假設μ的兩個可能值分別為μB1和μB2,相對應的兩類波的橫向位移幅值分別為α1和α2,則板的橫向位移可以表示為:

xj方向的轉角θj可由橫向位移的偏導給出:

考慮到線連接位于xj軸,省略時間項,則由式(10)和式(11)可求出:

線連接上xj軸的力矩與橫向力可表示為:

式中:v0為板的泊松比。將式(10)~(12)代入式(13)(14)中即可求得波數空間下線連接的動剛度矩陣為[12]:

式(15)即為應用波動理論得到的波數空間下線連接處的動剛度矩陣。注意式(15)為基于波數坐標k建立的,因此無法應用到混合FE-SEA分析中。另一方面從推導過程中可以看出,在基于波動理論推導線連接動剛度矩陣時,線連接沿xj方向為無限長,所以需要進一步處理。
由于有限元模型是基于節點坐標系的,所以需要將波數空間下的動剛度矩陣變換到節點坐標系下。設基于廣義坐標an的線連接的面外運動可表示為:

式中:U(x)=[w(x),θ(x)]T,φn(x)線連接的面外位移與轉角的形函數向量。假設在線連接外的位移為0,即線連接具有固支邊界條件,則線連接上直接場動剛度矩陣的第mn項{Ddir}mn可表示為[13]:

式中:L為實際結構中線連接的長度,fn(x)為線連接處的外力 f(x)在 φn(x)上的分量,對形函數φm(x)應用空間傅立葉變換:

在波數空間下線連接處位移形函數與其引起的相應分布力之間的關系可表示為:

式中:D∞(k)為波數空間下半無限板線連接的動剛度陣,由式(15)給出。對式(19)應用傅立葉逆變換可得:

將式(20)代入式(17)有:

式(21)給出了線連接直接場動剛度矩陣的第mn項的計算格式,可以看出線連接的直接場動剛度矩陣與形函數相關,形函數的構造方法是本文研究的重點。
由式(21)可以看出,形函數是直接場動剛度矩陣的計算的關鍵參數。本文提出了三種位移形函數構造方法,下面將對三種方法逐一進行介紹。
形函數的選取可基于線連接處的有限元模態分析結果。首先對混合模型中的有限元模型進行模態分析,提取線連接上與其節點相關的N階模態,然后對模態進行線性插值獲得第m個形函數φm(x):
式中:系數am和bm通過對第m個模態進行線性插值獲得,對式(22)進行傅里葉變換:
式(23)即為由模態j線性插值構造得到波數空間下的形函數。可以看出模態插值形函數的計算比較復雜,首先必須進行有限元分析提取線連接節點的模態,然后進行線性插值后求其傅里葉變換,最后才能代入式(21)中求解節點坐標系的動力學剛度矩陣。
模態函數的插值和傅里葉變換計算繁瑣復雜,形函數的數目取決于模態的數目,因此當模型模態密集且分析頻率范圍非常寬時,采用該方法建立線連接模型就非常困難了。為了解決這個問題,本節構造一種三角波線性插值形函數。構造方法如下:第m個形函數定義為線連接上節點m處的位移為1,其余節點位移為0,節點間的位移采用線性插值,如圖3所示。

圖3 三角波形函數Fig.3 Triangle wave shape function
若節點的位置為x=0時,形函數可寫為:

式中:Δl為節點間的距離。該三角波形函數的空間傅里葉變換可寫為:

根據傅里葉變換的時移特性,若第m個形函數的位置為x=xm,則有:

式(26)給出了三角波形函數的表達式。從推導過程可以看出,基于三角波插值的形函數特點在于形函數的數目由線連接上的節點數目決定,這是其優于模態插值形函數的一個重要特點,尤其是在寬頻帶分析時。但是注意到式(21)是奇異積分,因此采用上述兩種形函數構造方法可能存在積分不收斂的問題。
為了解決式積分由可能存在不收斂的問題,本節提出了一種基于Shannon小波 (Sinc函數)的線連接位移形函數構造方法,該函數的傅里葉變換在波數空間下是窗函數(如圖4所示),可有效簡化線連接建模并解決積分奇異問題。

圖4 Shannon小波函數及其傅里葉變換Fig.4 Shannon wavelet and it's Fourier transform
定義以Shannon小波函數為基的形函數為:

式中:ks與節點相關,本文取ks=2π/Δl。對上式進行傅里葉變換有:

由傅里葉變換的時移定理,第m個形函數φm(x)的傅里葉變換為:

將式(28)和式(29)代入式(21)可以得到基于Shannon小波形函數的線連接動剛度矩陣為:

式(30)表明利用Shannon小波函數的傅里葉變換特性,直接場的動剛度矩陣的積分變為有界積分,且積分只與節點間的相對位置有關,因此簡化了計算,可采用常規的積分算法求解式(30)。
為了校驗本文提出的混合線連接模型的有效性,設計如圖5所示的星體壁板與連接結構的組合模型對本文提出的混合線連接模型進行校驗。模型為星體壁板連接的簡化模型,兩塊星體壁板結構通過復雜的剛性結構連接,每個壁板結構與剛性連接結構之間通過線連接進行耦合。星體壁板結構楊氏模量為71 GPa,密度2700 kg/m3,泊松比 0.3296,剛性連接結構楊氏模量為 210 GPa,密度為 7800 kg/m3,泊松比0.3125,整個結構的損耗因子為1%,模型的幾何參數如表1所示。載荷為點載荷,垂直施加在其中一個星體壁板結構上,分析頻率范圍為1 Hz~1 kHz。考慮到整個系統的動力學特性,柔性較大的板面外位移采用SEA建模,剛性較大的板的面內位移和連接結構采用FEM建模,兩種模型通過兩個混合線連接耦合,如圖6所示。參考算法采用基于有限元的能量流分析與蒙特卡洛仿真(Monte Carlo Simulation)相結合的方法[14],并與國外商用軟件 VA One[15]的計算結果進行對比。

圖5 組合結構模型Fig.5 Combined structure model

圖6 混合線連接模型Fig.6 Hybrid line junction model

表1 結構的幾何參數Table 1 Geometrical parameters of the structure
為了描述中高頻響應對參數攝動的敏感性,在壁板結構中引入不確定性參數,本文中在每塊板上附加20個隨機分布的集中質量塊,集中質量塊總質量為壁板結構的15%,以有限元能量流分析—蒙特卡洛仿真得到系統的平均響應作為參考。圖7為兩個統計能量子系統間的耦合損耗因子計算結果。仿真結果表明,本文提出的基于Shannon小波的線連接建模方法和VA One的計算結果均與Monte Carlo仿真的結果都相吻合。但是VA One計算結果偏于保守,尤其是在非共振峰位置。
圖8給出了壁板子系統的無量綱響應結果(子系統能量與系統輸入能量的比值)。從圖中可以看出,對于能量輸入的壁板結構1,本文方法和 VA One均與Monte Carlo試驗相吻合,但是對壁板結構2,本文提出的方法在低頻處與Monte Carlo試驗存在誤差,其余頻段都比較吻合,而VA One誤差在非共振峰位置相對較大。低頻處的差異主要是由于在混合線連接建模中引入的固支邊界假設造成的。

圖7 基于Shannon小波方法的耦合損耗因子Fig.7 The coupling loss factor based on the Shannon wavelet method

圖8 基于Shannon小波方法的壁板響應Fig.8 The plates response based on the Shannon wavelet method

圖9 不同形函數的壁板之間的耦合損耗因子計算結果Fig.9 The coupling loss factor of two plates with different shape function

圖10 不同形函數的壁板響應預示結果對比Fig.10 The plates response results with different shape function
圖9 和圖10為基于不同形函數的線連接模型計算得到的兩個壁板之間的耦合損耗因子和壁板結構響應的對比曲線。從圖中可以看出,三種形函數模型下得到的結果基本吻合,但是模態形函數的誤差相對較大,這是由于模態截斷和線性插值引入的誤差導致的。從結果的對比可以看出,本文提出的三種形函數都可以有效建立混合線連接模型,對結構的響應進行正確預示。但是基于Shannon小波形函數的線連接模型構造最為簡單,能夠有效簡化計算過程,解決奇異積分不收斂的問題,因此更適合于航天器工程應用。
本文研究了航天器結構中頻段力學環境預示的混合線連接建模方法,并針對典型結構展開了仿真校驗,主要內容和結論如下:
(1)采用傅里葉變換方法,將波數空間下的線連接動剛度矩陣變換為有限元節點坐標系下的線連接直接場動剛度矩陣,傅里葉變換技術同樣也實現了將無限線連接模型轉化為符合實際結構的有限線連接模型。
(2)提出的三種形函數構造方法中,Shannon小波形函數線連接模型具有建模簡單、可解決動剛度矩陣積分奇異的問題等特點;模態插值線連接模型建模相對比較復雜,特別是當分析頻帶非常寬且模態密集時,該模型的求解效率急劇下降。
(3)算例仿真分析表明,本文提出的混合線連接建模方法具有較好的預示精度。同時對比三種形函數線連接模型的預示結果可以看出預示結果基本吻合,其中模態線性插值函數誤差相對較大。
(4)本文的驗證對象是簡化的星體壁板連接結構,下一步工作應當從算法優化、更加復雜的航天器結構的應用和試驗驗證兩方面展開深入研究。
[1] 馬興瑞,于登云,韓增堯,等.星箭力學環境分析與試驗技術研究進展[J].宇航學報,2006,27(3):323-331.[Ma Xing-rui, Yu Deng-yun, Han Zeng-yao, et al. Research evolution on the satellite-rocket mechanical environment analysis& test technology[J].Journal of Astronautics,2006,27(3):323 -331.]
[2] 馬興瑞,韓增堯,鄒元杰,等.航天器力學環境分析與條件設計研究進展[J].宇航學報,2012,33(1):1-12.[Ma Xing-rui,Han Zeng-yao,Zou Yuan-jie,et al. Review and assessment of spacecraft mechanical environment analysis and specification determination[J].Journal of Astronautics,2012,33(1):1 -12.]
[3] Shorter P J,Langley R S.Vibro-acoustic analysis of complex systems[J].Journal of Sound and Vibration,2005,288(3):669-699.
[4] Jeffrey M L,Cotoni V.Vibro-acoustic response of the NASA ACTS spacecraft antenna to launch acoustic excitation[R].Washington D C,Glenn Research Center,NASA/TM -2008-215168,1-15,2008.
[5] Chen SM,Wang D F,Zan JM.Interior noise prediction of the automobile based on hybrid FE-SEA method[J].Mathematical Problems in Engineering,2011,2011(4):1-20.
[6] Prasanth S,Charpentier A,Fukui K.Using the hybrid FE-SEA model of a trimmed full vehicle to reduce structure born noise from 200Hz to 1kHz[C].12th Symposium on International Automotive Technology,Pune,India,2011.
[7] 張瑾.FE-SEA方法在航天器力學環境預示中的應用研究[D].北京:中國空間技術研究院,2011.[Zhang Jin.Spacecraft mid-frequency dynamic predication using FE-SEA hybrid method[D]. Beijing:China Academy of Space Technology,2011.]
[8] 鄒元杰,韓增堯,張瑾.航天器全頻域力學環境預示技術研究進展[J].力學進展,2012,42(4):445-454.[Zou Yuanjie,Han Zeng-yao,Zhang Jin.Research progress on fullfrequency predication techniques of spacecraft’s mechanical environment[J].Advances in Mechanics,2012,42(4):445 -454.]
[9] 朱衛紅,馬興瑞,韓增堯.航天器中頻力學環境預示技術研究進展[J].航天器工程,2014,23(1):110-117.[Zhu Wei-hong,Ma Xing-rui,Han Zeng-yao.Research evolution on mid-frequency mechanical environment predication of spacecraft[J].Spacecraft Engineering,2014,23(1):110 -117.]
[10] Shorter PJ,Langley R S.On the reciprocity relationship between direct field radiation and diffuse reverberant loading[J].The Journal of the Acoustic Society of America,2005,117(1):85-95.
[11] Langley R S,Cordioli J A.Hybrid deterministic– statistical analysis of vibro-acoustic systems with domain coupling on statistical components[J].Journal of Sound and Vibration,2009,321(1):893-912.
[12] Langley R S,Heron K H.Elastic wave transmission through plate/beam junction[J].Journal of Sound and Vibration,1990,143(2):241-253.
[13] Cotoni V,Shorter P J.Numerical and experimental validation of a hybrid finite element-statistical energy analysis method[J].The Journal of the Acoustic Society of America,2007,122(1):259-270.
[14] Mace B R,Shorter P J.Energy flow models from finite element analysis[J].Journal of Sound and Vibration,2000,233(3):369-389.
[15] VA One 2010.5 user’s guide[Z].ESI Group,Paris,France,2010.