基于能量有限元法和虛擬模態綜合法的高頻沖擊響應分析方法

陳兆林，楊智春,*，王用巖，張新平

1. 西北工業大學 航空學院，西安 710072 2. 航空工業成都飛機設計研究所，成都 610091 3. 航空工業陜西飛機工業(集團)公司設計院，漢中 723213

在航空航天等工程領域中，沖擊載荷普遍存在。例如導彈發射[1]，艦載機彈射起飛和降落[2]，航天器采用爆炸螺栓進行級間分離[3-4]等，都會對結構施加短時、高幅值的沖擊載荷，使結構產生包含高頻成分在內的瞬態響應，造成結構尤其是敏感儀器和電子設備的故障或損傷，影響航空航天飛行器的正常工作。因此，發展高效的分析方法以求解結構在高頻沖擊載荷作用下的瞬態響應，具有明確的工程應用背景和意義。

在現有的沖擊響應分析方法中，有限元法(Finite Element Method, FEM)最為成熟且應用廣泛[5-6]，常見的商業有限元軟件基本都可以進行沖擊響應分析。但是對于高頻沖擊響應分析問題，FEM卻具有局限性。這是因為FEM模型需要滿足一個波長內至少6個單元，而在高頻段結構振動的波長很小，FEM模型就需要劃分非常細密的網格，導致模型規模過大，計算成本過高[7]。

為了進行高頻沖擊響應分析，Dalton[8-9]將統計能量(Statistical Energy Analysis, SEA)法和虛擬模態綜合(Virtual Mode Synthesis and Simulation, VMSS)法相結合，提出了SEA-VMSS法，并應用該方法分析了某裝甲車在沖擊載荷作用下控制艙的高頻響應[10]，通過與試驗測得的加速度沖擊響應譜進行對比，驗證了該方法的有效性。近年來，國內也開展了SEA-VMSS法的相關研究工作。楊博[2]采用SEA-VMSS法預測了艦船在沖擊載荷作用下結構的響應和艙室內部的噪聲，發現沖擊載荷的類型和作用時間對結構響應幅值和艙室噪聲都有很大影響。王軍評等[3]采用SEA-VMSS法對某航天器的部分艙段在爆炸分離時的沖擊響應進行了分析，發現結構上各個位置的加速度響應譜都含有較高的頻率成分，且加速度的幅值隨著與沖擊源距離的增加而減小。彭志剛[4]以某衛星平臺為研究對象，采用SEA-VMSS法分析了星箭分離時衛星不同部位的沖擊響應譜，研究了結構的材料屬性、爆炸分離方式等因素對沖擊響應譜的影響。從上述研究中可以看出，SEA-VMSS法綜合了SEA適用于高頻穩態分析而VMSS適用于瞬態分析的優點，成為高頻沖擊響應分析的高效途徑。但是，由于SEA只能獲得一個子系統的平均響應，不能得到響應在子系統內部的空間變化，導致SEA-VMSS法只能獲得一個子系統的平均沖擊響應，不能分析子系統內部某一點的沖擊響應。

為了克服SEA-VMSS法的上述缺點，發展新的高頻沖擊響應分析方法是必要的。能量有限元法(Energy Finite Element Method, EFEM)是近年來發展的另外一種適用于高頻穩態分析的新方法。EFEM以時空平均的能量密度為變量，具有模型簡單、求解迅速的優點，而且能夠得到子系統內部任一點的響應。對該方法的研究始于20世紀80年代，Nefske和Sung[11]提出時間平均能量強度正比于時間平均能量密度的梯度，并根據此假設推導得到類似于熱傳導方程的能量密度控制方程。Wohlever和Bernhard[12-13]從桿的縱向振動和梁的橫向振動方程出發，證明了Nefske和Sung[11]的假設，成功將EFEM應用于求解桿、梁模型的高頻振動響應。Bouthier和Bernhard[14-16]以時間和空間平均的能量密度為變量，推導得到薄膜、平板結構的能量密度控制方程并進行求解。上述研究中，載荷均為單點簡諧激勵，Han等[17-18]計算了分布載荷作用下梁和板的輸入功率，并發展了相應的EFEM。在考慮熱效應的EFEM研究方面，Zhang等[19]建立了受熱梁結構的能量密度控制方程。Wang等[20]將Zhang等[19]的方法推廣到受熱平板模型中，并考慮了非均勻溫度場的影響。經過科研工作者的努力，EFEM在高頻響應分析領域的應用越來越廣泛。但是，EFEM僅適用于穩態分析，不能進行沖擊載荷作用下的高頻瞬態分析。

為了分析結構在高頻沖擊載荷作用下的瞬態響應，本文通過理論推導，將EFEM和VMSS相結合，提出了EFEM-VMSS法，并通過算例驗證了該方法在高頻沖擊響應分析方面的有效性及其相對于FEM和SEA-VMSS法的優點。

1 EFEM-VMSS法的建立

1.1 能量有限元法和頻響函數

能量有限元法的核心是以能量密度為基本變量，描述系統內部的能量分布。下面以一維結構為例進行推導。

處于穩態振動的一維結構中，能量由傳播方向相反的2列行波(分別記為左行波和右行波)進行傳遞。在左行波中，任取一微元體，存在如圖1所示的能量平衡關系，該能量平衡關系的表達式為[21]

(1)

由于系統處于穩態振動，所以一個周期內的平均能量密度為常數，即?〈e-〉/?t=0，其中符號〈·〉代表對變量在一個周期內進行時間平均。對式(1)在一個振動周期內進行時間平均，得到

(2)

同理，在右行波中也存在類似的能量平衡關系，如圖2所示。由于右行波和左行波的能量強度方向相反，所以，右行波中的能量平衡關系為

圖1 左行波能量平衡關系示意圖Fig.1 Sketch of energy balance relationship in left-traveling wave

圖2 右行波能量平衡關系示意圖Fig.2 Sketch of energy balance relationship in right-traveling wave

(3)

(4)

式中：η為結構阻尼比；ω為振動的圓頻率。

另外，在每一列行波中，能量強度正比于能量密度[11]，即

〈I±〉=cg〈e±〉

(5)

式中：cg為能量傳播的速度，即群速度。

將式(4)和式(5)代入式(2)和式(3)可得

(6)

(7)

將左行波和右行波的能量密度進行疊加，可以得到結構各點的能量密度

〈e〉=〈e+〉+〈e-〉

(8)

由于左行波和右行波的能量強度方向相反，所以結構各點指向x正方向的凈能量強度為

〈I〉=〈I+〉-〈I-〉

(9)

需要說明的是，式(8)和式(9)忽略了左行波和右行波之間的干涉作用，如果考慮干涉作用的影響，會導致能量密度和能量強度的分布出現空間諧波分量，此時需要對能量密度和能量強度在一個波長內進行空間平均，消除空間諧波分量，式(8)和式(9)才能適用[12-13]。

根據式(8)和式(9)，用式(7)減式(6)可得

(10)

用式(3)減式(2)，并結合式(4)和式(10)，可得能量密度控制方程為

(11)

對式(11)可以采用有限元方法進行離散求解。采用伽遼金法，根據格林第二公式，可得每個單元內的能量有限元方程為

Κeee=Qe

(12)

式中：Κe為單元系數矩陣；ee為單元節點能量密度向量；Qe為進入單元的能量強度向量。

(13)

(14)

式中：Φ為單元形函數矩陣。

將各單元的能量有限元方程進行組裝，可以得到結構的總體能量有限元方程，然后進行求解，得到結構各點的能量密度。其中，單元向量Qe組裝成的總體向量Q的表達式為

Q=[0,…,0,P,0,…,0]T

(15)

式中：P為激勵點處的輸入功率。

輸入功率P的計算可以采用阻抗法[11,17-20]，即

(16)

式中：Frms為激勵力的均方根值；Z為輸入阻抗。在高頻分析中，常用無限大結構的阻抗來近似有限大結構的阻抗[17-18]，這樣得到的輸入功率并不是每一個激勵頻率下精確的輸入功率，而是輸入功率在每一個頻段內的平均值。

在穩態振動下，經過時間平均和空間平均后，動能密度等于勢能密度[12-13]。因此，通過能量有限元法獲得結構各點的能量密度后，可認為動能密度是能量密度的1/2，所以可根據式(17)求得結構各點的速度均方根值

(17)

式中：ρ為密度；S為截面面積。進而可以得到各響應點與激勵點之間的速度頻響函數為

Hv=vrms/Frms

(18)

需要說明的是，式(18)所得速度頻響函數Hv并不是精確的頻響函數，而是頻響函數在每一個頻段內的平均值。這是因為能量有限元法的推導過程中進行了時間平均、空間平均，而且計算輸入功率時采用了無限大結構的阻抗，所以能量有限元法的解并不是精確解，而是具有統計意義的結果，這一點與統計能量法類似。這種解在低頻段會產生較大誤差，但隨著分析頻率的升高，誤差會逐漸減小。同時，進行頻段平均等統計處理，需要保證頻段內有足夠多的模態數，才能具有足夠的精度。對于統計能量法，要求一個頻段內至少有5階模態，對于能量有限元法，目前還沒有確定的判據，但一般來說模態數越多精度越好。另外，在高頻分析中，精確地分析每一個頻率點處結構的響應會消耗大量的計算資源，而采用頻段平均等統計思想，既可以提高計算效率，又可以在初步設計階段較為準確地預示結構的高頻響應[1]，因此統計能量法和能量有限元法在高頻振動響應分析中獲得了廣泛的應用。

1.2 虛擬模態綜合法和瞬態響應分析

虛擬模態綜合法的基本思想是通過已知的頻響函數頻段平均值，獲得虛擬模態振型系數，來構造虛擬模態綜合系統。

在小阻尼假設下，結構上i點和j點之間在頻率ω處的速度頻響函數的幅值可以按照模態疊加的形式表示為

(19)

式中：pm為第m階模態頻率；ζm=η/2為第m階模態阻尼比；Ψim和Ψjm分別為第m階模態振型向量的第i個和第j個元素；n為分析頻帶內的模態數。

在每個頻段內取足夠多的頻率點，通過梯形積分，可以得到|Hij|在每個頻段內的平均值為

(20)

式中：Δω為頻段的帶寬；s為該頻段內頻率點的個數。

式(20)可以寫成矩陣的形式，即

BΨij=Hij

(21)

式中：Hij為速度頻響函數頻段平均值向量，可由能量有限元法獲得；Ψij為虛擬模態振型系數向量，是各階虛擬模態振型向量中與響應點和激勵點對應的元素的乘積，即

(22)

由于頻段數一般遠小于模態數，所以矩陣B并不是方陣，則需要通過求矩陣B的廣義逆矩陣的方式來求解式(21)，得到

(23)

需要說明的是，Ψij不是結構精確的模態振型系數，而是虛擬模態振型系數。這是因為，在高頻分析中結構模態密集，要精確計算結構的各階模態參數變得困難，但可以計算結構的模態密度，進而得到結構在每個頻段內的模態數，然后認為所有模態在每個頻段內是均布的，從而得到各階虛擬模態固有頻率。Ψij正是這些虛擬模態的振型系數。雖然虛擬模態不是結構真實的模態，但能夠確保虛擬模態綜合系統的頻響函數與結構精確的頻響函數相比具有相同的頻段平均值。因為虛擬模態綜合法采用了這種頻段平均的思想，所以只有在模態密集的高頻段才能保證足夠的精度。

獲得虛擬模態振型系數Ψij后，便可以通過模態疊加法得到結構上任意一點i在k個激勵下的位移時域響應為

(24)

式中：Rmj為第m階模態力等于Qj(t)時第m階模態位移的時域響應。

(25)

式中：Qj(t)為施加在j點的激勵力的時間歷程。

Rmj可以通過Duhamel積分進行計算，即

(26)

式中：h(·)為第m階模態的位移脈沖響應函數。

通過類似式(24)～式(26)的方法，還可以得到結構各點的速度和加速度響應。

1.3 EFEM-VMSS法

能量有限元法可以進行高頻穩態分析，得到結構頻響函數的頻段平均值，但無法進行瞬態響應分析。而虛擬模態綜合法則能根據頻響函數的頻段平均值，進行瞬態響應分析。因此本文提出將能量有限元法穩態分析的輸出結果，作為虛擬模態綜合法的輸入，從而完成高頻沖擊響應分析。如圖3所示，為EFEM-VMSS法的分析流程。

圖3 EFEM-VMSS法分析流程Fig.3 Flow chart of EFEM-VMSS

2 算例驗證

為了驗證本文建立的EFEM-VSMM法的正確性，并說明該方法的優點，下面以一個兩端簡支梁模型為研究對象，進行沖擊響應分析，并與傳統FEM和SEA-VSMM法的結果進行對比。簡支梁模型如圖4所示，梁長度l=4 m，密度ρ=2 700 kg/m3， 彈性模量E=71 GPa，結構阻尼比η=0.02，截面面積S=4×10-4m2，截面慣性矩Ib=1.33×10-10m4。圖4中b1、b2和b3為3個響應點，到梁左端的距離分別為3l/8、l/4和l/8，a1和a2為2個激勵點，到梁左端的距離分別為l/2和3l/4。

分析頻帶為891～14 130 Hz。工程中常按照倍頻程、1/3倍頻程或恒定帶寬將分析頻帶劃分為若干頻段，本文按照表1將分析頻帶劃分為12個1/3倍頻程。雖然第1個頻段內的模態數最少，為7個，但所有頻段均滿足統計能量法中一個頻段內至少5階模態的要求。

圖4 簡支梁示意圖Fig.4 Sketch of a pinned-pinned beam

表1 1/3倍頻程的中心頻率和帶寬

Table 1 Center frequencies and bandwidths of one-thirdoctave bands

頻段號中心頻率/Hz帶寬/Hz11 000891～1 12221 2501 122～1 41331 6001 413～1 77842 0001 778～2 23952 5002 239～2 81863 1502 818～3 54874 0003 548～4 46785 0004 467～5 62396 3005 623～7 079108 0007 079～8 9131110 0008 913～11 2201212 50011 220～14 130

2.1 穩態分析結果

圖5所示為激勵頻率f=10 kHz時梁上各點的橫向振動速度均方根值。從圖5中可以看出，SEA只能得到結構平均的速度均方根值，不能反映速度均方根值在結構內部的空間分布。

由于EFEM的推導過程中采用了局部空間平均，使得EFEM的結果并不表征某一點的速度均方根值，而是表征該點附近一個波長內平均的速度均方根值。因此，從圖5中可以看出，EFEM的結果相當于FEM的結果在一個波長內進行了空間平均，消除了其空間諧波分量，但可以表征速度均方根值在空間的總體變化趨勢。也正是因為EFEM只需要計算速度均方根值的總體變化趨勢，不需要計算空間諧波分量，所以EFEM模型所需的單元數量很少。

圖6和圖7所示為分析頻帶內各響應點與激勵點a1之間的速度頻響函數(FRF)。其中圖6為FEM和EFEM的結果對比，圖7為SEA和EFEM的結果對比。從圖6中可以看到，EFEM所得速度頻響函數與FEM所得速度頻響函數在每一個頻段內的平均值吻合良好，對于不同位置的響應點，EFEM的結果都可以表征FEM的結果在頻域的總體變化趨勢。需要說明的是，正如1.1節中所述，EFEM的結果是時間、空間和頻域上的統計結果，其求解精度隨著頻率增大而提高，因此，如圖6所示，除個別頻段外，總體上EFEM結果與FEM結果頻段平均值之間的誤差也隨頻率增大而減小。從圖7中EFEM的結果可以看出，隨著響應點與激勵點之間的距離逐漸增大，速度頻響函數逐漸減小，而SEA無法反映速度頻響函數在結構內部的空間變化。

圖5 速度均方根值分布Fig.5 Distribution of RMS of velocity

圖6 EFEM和FEM所得的速度頻響函數Fig.6 Mobility FRF obtained by EFEM and FEM

圖7 EFEM和SEA所得的速度頻響函數Fig.7 Mobility FRF obtained by EFEM and SEA

2.2 沖擊響應分析

本節利用2.1節中穩態分析的結果和虛擬模態綜合法，進行沖擊響應分析，在本節的分析中，僅在激勵點a1處施加如圖8所示的三角脈沖激勵。對于沖擊響應分析的結果，采用加速度沖擊響應譜(SRS)進行表達，加速度沖擊響應譜的計算采用的是遞歸數字濾波法[23]。

圖9和圖10所示為分析頻帶內各響應點頻段平均的加速度沖擊響應譜。其中圖9對比了FEM和EFEM-VMSS法的結果，圖10對比了SEA-VMSS法和EFEM-VMSS法的結果。從圖9中可以看出，通過EFEM-VMSS法得到的各響應點的加速度沖擊響應譜與FEM的結果吻合良好，從而說明了本文建立的EFEM-VMSS法的正確性。需要說明的是，由于EFEM-VMSS法是在EFEM所得頻響函數基礎上進行沖擊響應分析的，所以EFEM所得頻響函數的精度總體上決定著EFEM-VMSS法所得沖擊響應譜的精度，這一點結合圖6和圖9即可以看出。例如，8 kHz附近，圖6(a)和圖6(b)中EFEM和FEM所得的頻響函數頻段平均值誤差很小，所以圖9(a)和圖9(b)中2種方法所得沖擊響應譜峰值吻合良好。而圖6(c)中EFEM所得的頻響函數較大，從而導致圖9(c)中EFEM-VMSS法所得的沖擊響應譜峰值較大。

圖8 三角脈沖Fig.8 Triangular pulse

圖9 EFEM-VMSS法和FEM所得沖擊響應譜Fig.9 SRS obtained by EFEM-VMSS and FEM

圖10 EFEM-VMSS法和SEA-VMSS法 所得沖擊響應譜Fig.10 SRS obtained by EFEM-VMSS and SEA-VMSS

從圖10中EFEM-VMSS法的結果可以看出隨著響應點與激勵點之間的距離逐漸增大，加速度沖擊響應譜逐漸減小，而SEA-VMSS法無法反映結構內部各響應點沖擊響應譜的分布。

2.3 激勵位置的影響

改變激勵位置，僅在激勵點a2處施加如圖8所示的三角脈沖激勵，通過EFEM-VMSS法進行沖擊響應分析，并與激勵位置為點a1時的結果進行對比，如圖11所示。可以看出，隨著激勵點位置向右平移，即各響應點與激勵點之間的距離均增大，導致各響應點的加速度沖擊響應譜均減小。這表明EFEM-VMSS法可以分析不同激勵位置下結構內部各點的沖擊響應。而SEA-VMSS法受制于SEA理論，不能考慮激勵點位置在結構內部的變化，在一個子系統內部不同位置進行激勵，SEA-VMSS法所得結果是不變的。

圖11 不同激勵位置下EFEM-VMSS法所得沖擊響應譜Fig.11 SRS obtained by EFEM-VMSS at different exciting positions

3 結 論

1) 詳細推導了能量有限元法和虛擬模態綜合法的基本方程及其求解過程，通過將能量有限元法和虛擬模態綜合法相結合，提出了結構高頻沖擊響應分析的EFEM-VMSS法。通過算例分析，將EFEM-VMSS法的分析結果與FEM、SEA-VMSS法的結果進行了對比，驗證了本文方法能夠分析結構在高頻沖擊載荷作用下的瞬態響應。

2) 與FEM相比，EFEM-VMSS法具有模型簡單、求解速度快的優點。與SEA-VMSS法相比，EFEM-VMSS法能夠分析結構內部各點的沖擊響應，而且能夠考慮沖擊載荷的作用位置。因此，本文提出的EFEM-VMSS法在高頻沖擊響應分析方面具有更好的適用性。