基于malab的fe-aa混合方法的研究

基于malab的fe-aa混合方法的研究

基于波的fe-ea混合方法在高頻振動和噪聲負荷的作用下，每種復雜結構通常表現(xiàn)出相對復雜的動態(tài)特性。究其原因主要是由于各類不確定性引起的,如復雜結構高頻參數(shù)的不確定性、結構制造工藝誤差、振動載荷及聲場的不確定性等。根據(jù)結構特征尺寸與結構內波長關系的不同使得結構不確定性對結構響應的影響也有所不同[1~3]。結構特征尺寸小于或相當于結構內波長時,即結構的模態(tài)稀疏時,結構不確定性對響應的影響很小,甚至低頻時可忽略不計;反之,結構特征尺寸大于結構內波長時,即結構的模態(tài)密集時,不確定性對響應的影響較大,在高頻段結構響應往往呈現(xiàn)出隨機的特性。傳統(tǒng)的有限元方法(FiniteElementMethod,FEM)能很好地反映確定性對結構響應的影響而無法反映不確定性的影響,因此該方法僅適用于低頻段響應分析;而傳統(tǒng)的統(tǒng)計能量分析(StatisticalEnergyAnalysis,SEA)僅能反映不確定性對結構的影響,因此該方法僅適用于高頻段。對于復雜結構,各部組件特征尺寸與其結構波長的關系不同,如運載火箭上的桁條蒙皮結構,桁條的結構波長與其特征尺寸相當,而蒙皮的結構波長遠小于其特征尺寸,此時各部組件的不確定性對整體結構響應影響不一致,有的呈現(xiàn)出低頻特征,而有的呈現(xiàn)出高頻特征。單獨使用FEM或SEA均不能有效地對此類復雜結構的響應進行預示。對此類介于低頻FEM和高頻SEA的問題統(tǒng)稱為中頻問題。FE-SEA混合方法正是針對中頻問題的特點提出的一種方法。該方法可根據(jù)子結構特征尺寸與結構波長關系的不同選取不同方法進行建模,如子結構特征尺寸小于或相當于其結構內波長時,應用FEM或邊界元方法進行建模,稱其為確定性子結構(或確定性子系統(tǒng));而子結構特征尺寸大于結構內波長時,應用SEA進行建模,并稱其為隨機子結構(或隨機子系統(tǒng))。建模完成以后,通過確定性子結構與隨機子結構連接處直接場和混響場之間的互易原理將子結構重新連接,得到整體結構的響應。FE-SEA既是傳統(tǒng)FEM與SEA結合體,同時也是對傳統(tǒng)方法的改進,填補了傳統(tǒng)分析方法在中頻段應用的空白。而且,當頻率從中頻逐漸向低頻和高頻兩極趨近時,FE-SEA逐漸退化為FEM和SEA,也可以說FEM和SEA是FE-SEA的兩個特例。1999年,Langley和Bremner提出了FE-SEA混合方法的基本理論。在應用模態(tài)疊加法的同時,將振型及相應的主坐標分解為整體模態(tài)集和局部模態(tài)集,并對與兩個模態(tài)集相關的子系統(tǒng)分別應用有限元和統(tǒng)計能量分析進行建模。2005年,Langley和Shorter等在上述基于模態(tài)的FE-SEA混合法的基礎上,提出了基于波動理論的FE-SEA混合方法。在各類混合方法中,基于波動理論的FE-SEA混合法可以處理各類結構間載荷雙向傳遞問題,克服其他混合方法只能考慮能量單向傳遞的不足,使得該方法可應用于更為復雜的結構形式。同年,法國ESI集團也以基于波的FE-SEA混合法為理論基礎推出了PAM-VAOne商業(yè)軟件。此后,Cotoni和Shorter分別利用數(shù)值方法和試驗方法對FE-SEA方法進行了驗證。目前,伴隨VAOne商業(yè)軟件的推廣,基于波的FE-SEA混合法逐步應用到轎車隨機振動及噪聲、飛機艙內噪聲和火箭整流罩內噪聲分析及航天器聲振響應分析上[10~12],展現(xiàn)出良好的工程應用前景。本文首先闡述了基于波動理論的FE-SEA混合方法的基本理論,建立了該方法完整的分析流程。而后,應用Matlab開發(fā)了相應的計算程序并進行了校驗。最后,結合MonteCarlo仿真分析和基于有限元的能量流分析,通過實例對基于波的FE-SEA混合方法進行了仿真驗證。1基于波動理論的fe-ea混合方法1.1各子系統(tǒng)間邊界劃分考慮整體系統(tǒng)中的一個子系統(tǒng),如圖1所示,其區(qū)域為K,邊界為Γ。若穩(wěn)態(tài)諧波激勵加載在該子系統(tǒng)邊界區(qū)域上,則其邊界上的時域位移響應可表示為ub(t)=Re{ubexp(ikt)},其中ub為與廣義邊界自由度對應的幅值,通常為復數(shù)向量。按照邊界元法的思想,子系統(tǒng)邊界區(qū)域的響應可由格林函數(shù)積分求得。為簡化計算,可將邊界離散為一系列廣義坐標qb,此時邊界位移可由式(1)求得,即式中ub(x∈Γ)為邊界x處的位移,Okb為邊界上的基函數(shù),qkb為第k個邊界廣義坐標。只考慮子系統(tǒng)間的相互作用關系,忽略外載荷對邊界的作用,子系統(tǒng)邊界的動力學方程可表示為式中qb為子系統(tǒng)的邊界廣義坐標,fb為邊界廣義力(只表示其他子結構給予的反作用力)。若子系統(tǒng)中,一部分邊界上所有的物理性質已知,而另一部分邊界物理性質只是部分已知或完全無法確定,如圖2所示,則可將子系統(tǒng)的邊界劃分為確定性邊界Γd和隨機邊界Γr。子系統(tǒng)均可通過確定性邊界Γd和隨機邊界Γr與其他子系統(tǒng)相連并通過這兩類邊界進行能量交換,確定性邊界Γd和隨機邊界Γr也均可承受外部載荷。此外,由圖2可知,確定性邊界Γd無需滿足邊界連續(xù)的條件,任意子系統(tǒng)間的連接均可應用在這兩類邊界上。子系統(tǒng)邊界進行能量交換的實質就是通過邊界的位移場將能量傳遞到其他子系統(tǒng)上或接受其他子系統(tǒng)傳遞給自身的能量。對應子系統(tǒng)邊界的分類,可將邊界的位移場分為兩類:直接場(DirectField)與混響場(ReverberantField)。直接場僅滿足確定性邊界上的邊界條件,并不考慮隨機邊界上的邊界條件。通過直接場可精確表示出確定性邊界上的輸出位移場,即通過直接場子系統(tǒng)可向外輻射能量;但通過直接場表示的隨機邊界上的位移場一般很復雜,暫且記為O,如圖3所示。為滿足隨機邊界上的邊界條件,定義第二個位移場,即混響場?；祉憟龅淖饔镁褪菍⒒祉憟鲋械倪吔鐥l件與直接場中的邊界條件線性疊加后使確定性邊界條件和隨機邊界條件同時得到滿足。因此,混響場必須滿足兩個條件:1)混響場中確定性邊界上的位移為0;2)與直接場線性疊加后隨機邊界條件得到滿足,例如,若隨機邊界為固支,則混響場中確定性邊界上的位移為0,隨機邊界上的位移為-O。對應于邊界分類及直接場和混響場的定義,可將邊界廣義坐標qb分離為qdb和qrb,相應地將邊界廣義力fb分離為fdb和fbrev,則子系統(tǒng)邊界的動力學方式可分解為由于隨機邊界上部分信息未知,因此有必要對式(3)進行化簡。由式(3)的方程(1)可得式中Dbdir定義為直接場的動剛度陣并用于求解確定性邊界上的位移,其表達式為對于每一項均為實數(shù)的基函數(shù),矩陣Dbdir為復數(shù)對稱矩陣。Dbdir通?？捎蛇吔缭ㄇ蟮?但對于某些規(guī)則的連接方式,如點、線和面連接,與此類規(guī)則的連接方式相關的直接場的動剛度陣Dbdir可由理論公式直接推導出解析解,即利用邊界連接處的位移協(xié)調關系及波在結構中傳播的性質推導出解析解[14~18]。1.2混響場建模與約束在基于波動理論的FE-SEA混合方法中,一個復雜系統(tǒng)通常可分為若干子系統(tǒng)。據(jù)子系統(tǒng)的特征尺寸與其系統(tǒng)中波長的關系可將子系統(tǒng)分為兩類:若其特征尺寸與其系統(tǒng)中的波長相當,即該子系統(tǒng)的剛度很大、模態(tài)稀疏,系統(tǒng)中的一些不確定因素,如制造公差等,對其響應不產生影響,則該子系統(tǒng)為確定性子系統(tǒng),可用FEM建模;若其特征尺寸大于其系統(tǒng)中的波長,即該子系統(tǒng)的柔性很大、模態(tài)密集,其響應對系統(tǒng)中不確定因素的變化很敏感,如相同設計同一批次的兩個結構,由于制造公差等因素使兩個結構的響應變化很大,則此類子系統(tǒng)為隨機子系統(tǒng),可用SEA建模。同時,對應地在數(shù)值上對各系統(tǒng)做離散化處理,系統(tǒng)響應可由一系列位移廣義坐標q表示,其中確定性子系統(tǒng)的廣義坐標為qd。類比有限元頻響分析,經(jīng)傅里葉變換FE-SEA混合方法中系統(tǒng)的動力學方程為式中fext為確定性子系統(tǒng)所受外力的向量;frev(m)為混響場受擋力;m為與確定性子系統(tǒng)相連的隨機子系統(tǒng)的個數(shù);Dtot為系統(tǒng)的總動剛度陣,是確定性子系統(tǒng)的動剛度陣Dd與直接場動剛度陣Ddir(m)的線性疊加Dd可直接由有限元求得,而Ddir(m)可由相應的Dbdir進行坐標變換求得。在式(6)中,直接場動剛度陣Ddir(m)表示確定性子系統(tǒng)通過直接場向隨機子系統(tǒng)傳遞的能量以及隨機子系統(tǒng)間通過確定性子系統(tǒng)進行的間接能量傳輸,frev(m)表示隨機子系統(tǒng)通過混響場向確定性子系統(tǒng)傳遞的能量。由于結構不確定性的存在,類比統(tǒng)計能量分析中求集合平均的概念,將式(6)改寫為互譜的形式并求集合平均得其中符號·-H表示矩陣的共軛轉置并求逆的運算。假設系統(tǒng)具有最大熵(MaximumEntropy),即系統(tǒng)的平均值為所有可能出現(xiàn)樣本的平均,簡言之系統(tǒng)具有最大平均信息量。此時,混響場表示所有混響場樣本的集合平均,隨機子系統(tǒng)在混響場中的受擋力frev(m)與系統(tǒng)不確定性因素的變化無關,因此與frev(m)相關的集合平均趨近于下面的極限值式中Tm為與混響場振幅相關的比例常數(shù),其求解如下式所示式中Em和nm分別為第m個隨機子系統(tǒng)在混響場中的所具有的能量和第m個隨機子系統(tǒng)的模態(tài)密度。通過子系統(tǒng)間的功率平衡關系可求解各隨機子系統(tǒng)的能量Em。若系統(tǒng)中包含p+q個隨機子系統(tǒng),其中p個隨機子系統(tǒng)與確定性子系統(tǒng)相連,其余的q個子系統(tǒng)只與隨機子系統(tǒng)相連而與確定性子系統(tǒng)不相連,則隨機子系統(tǒng)間的功率平衡方程為式中htot,m表示第m個子系統(tǒng)在混響場中的能量損耗系數(shù),即隨機子系統(tǒng)通過混響場向確定性子系統(tǒng)輸出能量;hmn為隨機子系統(tǒng)m與n通過確定性子系統(tǒng)建立耦合關系的耦合損耗因子;Zmk為能量從第m個隨機子系統(tǒng)直接傳遞到第k個隨機子系統(tǒng)時的耦合損耗因子;Zm為第m個隨機子系統(tǒng)的內損耗因子;Pin(m,0)表示確定性子系統(tǒng)通過直接場傳遞到第m個隨機子系統(tǒng)上的功率;Pin(m,1)表示外界直接對第m個隨機子系統(tǒng)的輸入功率;Em和En分別為第m個隨機子系統(tǒng)和第n個隨機子系統(tǒng)的能量;nm和nn分別為第m個隨機子系統(tǒng)和第n個隨機子系統(tǒng)的模態(tài)密度。htot,m,hmn及Pin(m,0)分別為1.3隨機子系統(tǒng)相關參數(shù)分析綜上所述,對一個計算頻段基于波動理論的FE-SEA混合方法的分析過程可大體分為7步(如圖4所示)。(1)系統(tǒng)的劃分首先,將系統(tǒng)劃分為確定性子系統(tǒng)和隨機子系統(tǒng);而后,劃分子系統(tǒng)的確定性邊界和隨機邊界并確定邊界上的連接處及連接方式;最后,選取確定性廣義坐標qd。(2)邊界連接處的動剛度陣的獲取規(guī)則點、線和面連接處的直接場的動剛度陣Dbdir可直接得到解析解。據(jù)確定性廣義坐標qd的選取,動剛度陣Dbdir需進行相關的坐標變換。(3)與確定性子系統(tǒng)相關的參數(shù)分析應用FEM計算與確定性子系統(tǒng)相關的動剛度陣Dd,據(jù)選取的確定性廣義坐標qd,相應地對動剛度陣Dd進行坐標變換;然后,將動剛度陣Dd與Ddir(m)線性疊加求得系統(tǒng)的總動剛度陣Dtot。同時,將確定性子系統(tǒng)的外載荷fext變換為互功率譜Sffext并進行相應的坐標變換。(4)與隨機子系統(tǒng)相關的參數(shù)分析與隨機子系統(tǒng)相關的參數(shù)分析包括:計算確定性子系統(tǒng)對與其相連的隨機子系統(tǒng)的輸入功率;計算各類隨機子系統(tǒng)間的耦合損耗因子;計算隨機子系統(tǒng)在混響場中的能量損耗系數(shù);計算所有隨機子系統(tǒng)由于自身的內損耗而損失的功率;將直接作用在隨機子系統(tǒng)上的外力向量變換為相應輸入功率Pin(m,1)。(5)隨機子系統(tǒng)響應的求解據(jù)步驟(4)中的所得參數(shù)建立隨機子系統(tǒng)的功率平衡方程,求得各隨機子系統(tǒng)能量的集合平均。(6)確定性子系統(tǒng)響應的求解據(jù)步驟(5)中的所得隨機子系統(tǒng)的能量求得隨機子系統(tǒng)通過混響場作用在確定性子系統(tǒng)上的力譜;然后據(jù)步驟(3)建立的動力學求解確定性廣義坐標qd的集合平均值〈Sd,qq〉。(7)重復步驟(2)~(6),計算下一個頻段響應的集合平均。2計算機器的開發(fā)根據(jù)上述的FE-SEA混合方法分析流程,以Matlab為平臺開發(fā)了相應的計算程序。而后,通過實例將混合方法的結果與MonteCarlo仿真分析的結果相比較對基于波的FE-SEA混合方法進行仿真驗證。2.1響應系數(shù)和響應曲線如圖5所示的板梁組合結構,結構的具體參數(shù)如表1所示,具有16個圓孔的Z型梁分別通過兩個點和3個點與平板1和2相連。結構中Z型梁的一端施加固支邊界條件,在板1的板面垂直方向施加空間位置不相關的寬頻隨機分布式激勵。然后分別應用FE-SEA混合方法和MonteCarlo分析計算出板1受激勵時隨機子系統(tǒng)的能量影響系數(shù)(EnergyInfluenceCoefficient,EIC)并進行比較。其中,計算頻段取值為1Hz的等帶寬,中心頻率范圍為1~1500Hz。兩類方法的計算過程中均選用模態(tài)為廣義坐標,模態(tài)截斷位置為3500Hz,即選用3500Hz之前結構的所有模態(tài)進行計算。由于在3500Hz內,Z型梁僅有28階模態(tài),模態(tài)非常稀疏,因此將梁建為確定性子系統(tǒng);當中心頻率為200Hz時,如表2所示,由板1的波數(shù)推導出板1彎曲、拉伸和剪切方向的波長分別為0.222,27.165和15.198m,與板1的對角線長度相比,彎曲方向的波長約為對角線的一半,而拉伸和剪切方向遠大于對角線長度。對于該結構中的薄板,當某一方向的波長小于對角線長度的2/3時,該方向的位移可見為隨機子系統(tǒng)。因此,可將板1面外彎曲建為隨機子系統(tǒng),而面內的拉伸和剪切方向的位移建為確定性子系統(tǒng)。同理,對板2各方向的位移進行處理。綜合考慮該組合結構,則可將兩塊板彎曲方向的位移分別建為兩個隨機子系統(tǒng)1和2,面內位移與Z型梁組合建為確定性子系統(tǒng)。2.2不同集中質量對模擬結構的影響FE-SEA混合方法的結果是集合平均的結果,因此應用MonteCarlo仿真分別對300個樣本進行分析,然后對這300次的結果進行平均得到結構的隨機平均值,并與FE-SEA的結果比較來驗證FE-SEA方法。對于任意單次樣本,作用在板1上的輸入功率P1和兩個隨機子系統(tǒng)的能量均是通過基于有限元的能量流分析求得。在MonteCarlo仿真過程中,隨著結構中隨機樣本數(shù)的不斷增加,樣本平均結果的穩(wěn)定性越來越強,即當樣本具有最大包絡時仿真中的集合平均結果與結構的不確定性無關,集合平均結果唯一確定。由于結構中的不確定性的影響主要是由隨機子系統(tǒng)的不確定性引起的,因此在增加樣本的包絡方面主要是模擬隨機子系統(tǒng)不確定性的變化。為模擬結構樣本的最大包絡選擇在隨機子系統(tǒng)的隨機位置上添加集中質量。集中質量主要作用是增大隨機子系統(tǒng)低階模態(tài)對外界擾動的敏感性從而增大整個隨機子系統(tǒng)對外界擾動的敏感性(中高頻時,隨機子系統(tǒng)本身對外界擾動已經(jīng)很敏感)。建立MonteCarlo仿真中使用的FE模型(未添加集中質量),如圖6所示,整個組合結構上應用5148個Tria3殼單元,總的節(jié)點自由度數(shù)為16866。對組合結構進行模態(tài)分析并提取出3500Hz以內結構的前342階模態(tài),圖7所示為FE模型(未添加集中質量)的第100階模態(tài)?？紤]計算成本,仿真過程中選擇FEM中的模態(tài)法并對模型進行300次MonteCarlo仿真。仿真的每個樣本即對圖6所示FE模型,分別在兩平板上添加10個集中質量,所有集中質量的總和為整個組合結構質量的25%,所在位置為FE模型中兩板的節(jié)點位置。令每塊平板上集中質量所在節(jié)點的節(jié)點號的概率分布為均勻分布,應用Matlab分別生成每塊板上的300組節(jié)點號隨機數(shù)。據(jù)隨機數(shù)每次生成集中質量的節(jié)點位置并據(jù)該節(jié)點位置相應地將集中質量添加到組合結構的FE模型中,然后對樣本進行模態(tài)分析并提取出3500Hz以內樣本的所有模態(tài)。在此基礎上,對樣本FE模型,在板1每個節(jié)點上(除點連接所在節(jié)點)施加垂直于板面方向的外力自譜,其中外力自譜在300次仿真中保持不變,而后利用模態(tài)法求得樣本響應。最后利用基于有限元的能量流分析推導出兩個隨機子系統(tǒng)的能量E1和E2、輸入功率P1及相應的EIC。300次MonteCarlo仿真時間約為5天,計算周期較長。圖8所示為仿真過程中板2上的隨機子系統(tǒng)EIC的兩個單次樣本。圖9和10分別顯示了300次MonteCarlo仿真的結果,從圖中也可看出單次仿真結果間的差異明顯。圖9和10中,由300次仿真的平均結果可以看出在900Hz附近,即第100階模態(tài)附近,兩子系統(tǒng)間的能量傳輸明顯,板1(激勵板)的EIC下降即板1的能量向外傳輸時,板2(接收板)的EIC顯著增大即板2接受了外部傳入其子系統(tǒng)上的能量。2.3fe-ea混合法對FE-SEA混合方法,為得到確定性子系統(tǒng)的FE模型,將圖7中FE模型上兩塊板的殼單元換為膜單元,梁上的殼單元不變。為了保證板上面內位移的質量陣,只考慮板的面內質量需單獨建立膜單元的質量陣,因此在FE模型中可直接在節(jié)點相應的面內位移上添加集中質量自行建立膜單元的質量陣。求解之前通過Nastran輸出確定性子系統(tǒng)的剛度陣和質量陣,并導入Matlab程序作為輸入;同時還需將MonteCarlo仿真中所有的節(jié)點外力自譜通過點連接的動剛度陣將其轉化為FE-SEA混合法中隨機子系統(tǒng)的外部輸入功率。在此基礎上,應用FE-SEA混合法通過Matlab程序求得隨機子系統(tǒng)的能量和相應的EIC。FE-SEA混合方法的計算過程中使用的確定性子系統(tǒng)的模態(tài)為3500Hz內的共28階模態(tài),整個的分析過程約40s左右,兩項指標均遠小于MonteCarlo仿真。如圖9和10所示,由FE-SEA混合法所得的集合平均值與MonteCarlo仿真的樣本平均值結果基本一致,尤其是中高頻段,如900Hz附近,兩類結果吻合得很好。但同時,低頻段兩類結果間存在誤差,可能的原因包括:a)FE-SEA求解過程中多次求逆引入