一種高效的發(fā)動機輻射噪聲計算方法研究

馮仁華，劉敬平，付建勤，鄧幫林，張大鳴

（1．湖南大學 先進動力總成研究中心，長沙 410082；2．美國加州州立大學 弗雷斯諾分校，弗雷斯諾CA 93726）

隨著社會的發(fā)展，汽車的噪聲振動（Noise，Vibra-tion，Harshness）水平是車輛的一個重要指標，NVH性能直接反應車輛的舒適性品質，是否具有良好的NVH性能將是決定能否被消費者接受的重要條件。發(fā)動機作為汽車的主要振動噪聲源，其噪聲振動的控制［1］尤其是輻射噪聲非常重要［2－3］。控制和降低噪聲的方法很多，最有效的方法是設計出合理的結構，使其在影響噪聲的重要頻率范圍內產生較少的振動。因此，噪聲振動的計算和分析非常重要［4］。

數(shù)值模擬分析技術廣泛用于發(fā)動機噪聲振動的計算、分析和優(yōu)化。目前噪聲振動的計算和分析方法主要有有限元法（Finite Element Method）、邊界元法（Boundary Element Method）和快速多級邊界元法（Fast Multipole Boundary Element Method）。

有限元法把連續(xù)的彈性體劃分成有限個單元，通過劃分網格建立有限元模型，計算系統(tǒng)的應變和應力以及動力學特性［5］。有限元法在處理三維聲輻射問題時，其網格劃分、數(shù)據(jù)準備和計算工作量大，同時它采取邊界截斷來計算無限域的聲輻射問題會帶來不小的誤差［6］。

邊界元法是利用邊界積分方程解決邊界值或初始值問題的數(shù)值計算方法。它具有只在邊界離散、降低問題維數(shù)、建模效率高、求解精度高等特點。非常適合處理復雜的邊界問題和無限域問題［7－9］。但當求解問題的規(guī)模增大時，求解效率下降很快，存儲量和計算量已成為制約邊界元法發(fā)展的瓶頸［10］。

快速多級邊界元法是伴隨快速多級算法發(fā)展起來、利用快速多級算法加速求解的邊界元新型求解方法。雖然快速多級邊界元法能夠解決傳統(tǒng)邊界元法所存在的問題［11］，但是快速多級邊界元法存在低頻時計算效率低，高頻時精度低的問題。

表面振動速度法是一個簡單的聲學計算方法，結構表面噪聲聲壓與結構表面振動速度成正比，利用這個正比關系，可以建立結構表面振動與表面聲壓之間的關系，并可以通過表面振動速度求解輻射聲的聲壓和聲功率［12］。目前，表面振動速度法常用于發(fā)動機噪聲的測試［3］，但是用該簡單方法進行模擬計算的研究和成果較少。美國三大汽車制造商和供應商應用過一些簡單的方法和策略進行發(fā)動機輻射噪聲計算，但是他們的計算結果都沒有確切的物理意義，只是一種和聲輻射功率相關的“指標（Indicator）”，只是做相對對比。本文在表面振動速度法理論基礎上，利用Matlab程序，通過努力改進計算方法以控制計算誤差和加快計算速度，開發(fā)了一款簡單的計算表面輻射噪聲聲功率的軟件，該方法和軟件與傳統(tǒng)的計算方法及商用軟件在計算效率和精度上相比較，各有優(yōu)劣。商業(yè)軟件的計算方法復雜，計算速度慢，但能夠求解所有的聲學未知量。而該軟件和方法能夠大幅度提高計算速度，但只能計算“聲輻射功率”，不能計算其它聲學未知量。因此，如果分析的目的是改進結構設計，降低結構噪音，本文基于表面振動速度法開發(fā)的軟件可以在保證精度的基礎上大大加快計算時間，降低產品開發(fā)周期。

1 發(fā)動機輻射噪聲計算原理

1.1 發(fā)動機表面振動與輻射噪聲之間的關系

發(fā)動機表面振動和輻射噪聲之間存在一定的關系，在各種機械激勵和燃燒激勵的作用下，發(fā)動機及其各部件必然發(fā)生表面振動，這種振動向外傳播，使其發(fā)動機表面振動能量的一部分作為聲能輻射出去，這就是發(fā)動機表面輻射噪聲的能量來源［3，13］。

當發(fā)動機處于穩(wěn)定工作狀態(tài)時，其表面振動速度和輻射聲功率可認為有如下關系［14－15］。

式中：Wrad為聲功率；ρ0c為聲輻射阻抗；S為噪聲輻射表面積為噪聲輻射表面法向振動速度的平方對時間及振動表面的平均值；σ為輻射比。

取基準聲功率W0＝10－12W，則A計權的聲功率為：

式中：Lw為聲功率級；ρ0c、S可通過計算得到；Δ為A計權網絡的衰減量〉可以通過模擬計算或者發(fā)動機表面振動試驗測量得到。

1.2 輻射比的計算

輻射比表示部件由于表面振動而輻射聲能的有效程度，它不僅與部件的形狀和邊界條件有關，還與振動頻率有關。一般可以認為在臨界頻率以上輻射比σ≈1，在臨界頻率以下輻射比為0～1，而且隨頻率的下降而變得越來越?。?3］。

一般情況下，發(fā)動機各零部件形狀復雜，輻射比的精確理論計算十分困難，常采用實驗的方法或經驗公式計算的方法得到輻射比。對于NVH的CAE計算和分析，首要目標是減振降噪，而不是單純的為了計算得到非常精確的數(shù)值，所以模擬計算相對對比結果是非常重要的。在采用A計權聲功率級時，可以假設在所有頻率上σ≈1而使得測量結果不致有較大誤差［16］，因為評價發(fā)動機聲功率時一般用A計權聲功率級，這種計權網絡在低頻段要衰減一定的分貝值，衰減后一般低頻段能量在總噪聲能量中所占比重較低，所以輻射比在低頻段上的誤差對總聲功率級測量結果影響并不是很大。車用發(fā)動機的輻射噪聲頻率范圍主要在500～3 000 Hz，而其主要噪聲輻射部件的臨界頻率大致在 500～800 Hz［17］。

2 M atlab程序軟件介紹及發(fā)動機輻射噪聲計算分析流程

2.1 M atlab程序軟件介紹

本文以式（1）～（2）為理論基礎，利用 Matlab程序，開發(fā)了發(fā)動機表面輻射噪聲計算分析軟件，在該Matlab程序軟件中，式（1）中的輻射比σ取值為1。該Matlab程序計算軟件界面如圖1所示。

圖1 Matlab程序計算分析軟件界面Fig．1 The interface of Matlab program software

按照Matlab程序軟件要求，需輸入要計算分析的發(fā)動機零部件或者動力總成的表面網格節(jié)點信息文件、表面網格文件和表面網格速度信息文件。表面網格節(jié)點信息文件、表面網格文件可以通過Hypermesh軟件獲得，表面網格速度信息文件可以通過NASTRAN計算得到。通過Matlab程序軟件計算可以得到輻射噪聲聲功率的頻譜積分曲線、1／3倍頻程曲線和窄頻帶頻譜曲線，如圖2所示。

圖2 Matlab程序計算分析得到的結果曲線Fig．2 Radiated sound power spectrum curves from Matlab program software

圖2中，頻譜積分曲線表示從結構輻射出的聲功率在頻率上的積分，該曲線在聲學優(yōu)化方面非常重要。1／3倍頻程頻譜曲線主要用于結構改進的分析和比較。窄頻帶曲線可以追蹤輻射噪聲能量峰值點，可為下一步結構優(yōu)化提供指導和依據(jù)。

2.2 發(fā)動機輻射噪聲計算分析流程

發(fā)動機輻射噪聲計算分析流程如圖3所示。

圖3 發(fā)動機輻射噪聲計算分析流程Fig．3 Powertrain radiated noise analysis produre

圖3中，首先利用發(fā)動機部件或動力總成的幾何模型建立網格模型，并施加邊界條件，如材料屬性、約束、載荷、阻力等。然后利用ABAQUS或者NASTRAN計算得到所需的振動信息，如速度、加速度或位移。最后利用Matlab程序軟件進行輻射噪聲聲功率計算和頻譜分析。噪聲計算和分析的主要目的是在指定的頻帶下降低輻射噪聲聲功率，以達到目標要求。如果計算和分析的結果能達到目標要求，那么整個輻射噪聲流程就結束，如果不能達到目標要求，就需要進行頻譜分析，然后通過不同的方法優(yōu)化結構以達到目標要求。

3 數(shù)值算例

本文使用所編寫的Matlab程序軟件和LMS virtual lab 10中的傳統(tǒng)邊界元、快速多級邊界元模塊，利用傳統(tǒng)邊界元法、快速多級邊界元法和Matlab程序三種不同的計算方法分別對普通平板、發(fā)動機缸蓋罩、發(fā)動機動力總成的輻射噪聲功率進行計算分析。所有的計算在配置為Intel（R）Core（TM）i7－2600 CPU＠3．4處理器和32 GB內存的戴爾Vostro 460臺式機上進行。所有算例以空氣為聲學介質，其中空氣密度為1．29 kg·m－3，空氣中聲波波速為340 m·s－1。

3.1 平板

以受垂向點力激勵的矩形鋁板在空氣中的振動聲輻射為例進行計算，該平板的外形尺寸及計算邊界條件如圖4所示。

圖4 平板外形尺寸及振動邊界條件Fig．4 The dimensions and vibration boundary conditions of the plate

在圖4中，t為平板厚度，E為材料楊氏模量，μ為泊松比，ρ為鋁板密度，ξ為鋁板的臨界阻尼。強迫振動響應頻率為8～4 000 Hz，步長為8 Hz。

首先利用NASTRAN計算鋁板網格各節(jié)點的速度邊界條件，然后分別采用Virtual lab 10中的邊界元法、快速多級邊界元法和編寫的Matlab程序軟件進行輻射噪聲聲功率計算，計算結果如圖5所示。

不同計算方法得到的平板輻射聲功率曲線Fig．5 Radiated sound power spectrum curves of plate by different calculationmethods

從圖5中可以看出，用邊界元法、快速多級邊界元法和Matlab程序軟件計算得到的各頻率下輻射噪聲聲功率數(shù)值在低頻段有一定的差距，高頻段這種差距逐漸減小。整體上看，Matlab程序軟件計算得到的聲功率最高，邊界元法計算結果次之，快速多級邊界元法最低，這是因為在Matlab程序軟件中輻射比在所有頻率上取值為1。三種方法計算得到的峰值聲功率所對應的頻率點相同。

3.2 發(fā)動機缸蓋罩

發(fā)動機缸蓋罩是重要的輻射噪聲源，有時可占到發(fā)動機表面輻射噪聲總能量的30%以上［2］。計算和分析缸蓋罩的輸出聲功率頻譜，對于減小振動和降低發(fā)動機整機噪聲是非常有益的。本算例以某四缸發(fā)動機的缸蓋罩為例，計算缸蓋罩的輻射輸出聲功率頻譜。該發(fā)動機的動力總成網格模型如圖6所示。缸蓋罩表面網格模型如圖7所示。

圖6 發(fā)動機動力總成網格模型Fig．6 The engine powertrain grid model

圖7 缸蓋罩表面網格模型Fig．7 The surface grid model of engien head cover

圖6中的發(fā)動機動力總成網格模型主要包括：缸體、缸蓋、缸蓋罩、前端蓋、油底殼、軸承蓋、變速箱。圖7中缸蓋罩表面網格節(jié)點數(shù)為6 345。在圖6中的發(fā)動機動力總成模型中施加材料、載荷、連接等邊界條件，利用NASTRAN可以計算得到缸蓋罩8～2 000 Hz的表面速度邊界條件，最后計算缸蓋罩的聲輻射功率。用邊界元法、快速多級邊界元法以及Matlab程序軟件計算得到的缸蓋罩輻射噪聲輸出聲功率頻譜曲線，如圖8所示。

從圖8可以看出，用邊界元法、快速多級邊界元法和Matlab程序軟件計算得到的各頻率下輻射噪聲輸出聲功率數(shù)值在趨勢上和平板一樣，在低頻段有一定的差距，高頻段這種差別逐漸減小，Matlab程序軟件計算得到的聲功率最高，邊界元法計算結果次之，快速多級邊界元法最低，三種方法計算得到的峰值聲功率所對應的頻率點相同。

圖8 不同計算方法得到的缸蓋罩輻射聲功率曲線Fig．8 Radiated sound power spectrum curves of engine head

3．3 發(fā)動機動力總成

利用圖6中發(fā)動機動力總成網格模型提取發(fā)動機表面網格，發(fā)動機動力總成表面網格包括缸蓋罩、缸蓋、缸體、前端蓋、油底殼和變速箱。動力總成表面網格節(jié)點數(shù)為23 662。利用NASTRAN計算得到發(fā)動機動力總成8～2 000 Hz的表面速度邊界條件，最后計算發(fā)動機動力總成輻射聲功率。用邊界元法、快速多級邊界元法以及Matlab程序軟件計算得到的發(fā)動機動力總成表面輻射噪聲輸出聲功率頻譜曲線，如圖9所示。

圖9 不同計算方法得到的動力總成輻射聲功率曲線Fig．9 Radiated sound power spectrum curves of engine powertrain by different calculationmethods

圖9中，用三種不同方法計算得到的發(fā)動機動力總成輻射噪聲輸出聲功率頻譜曲線與平板算例和缸蓋罩算例在趨勢上一致，Matlab程序軟件得到的聲功率最高，快速多級邊界元計算得到的聲功率最低。在1 000 Hz以下，三種方法計算得到的輸出功率有一定的差別，但是隨著頻率的增加，這種差別越來越小。另外，用Matlab程序軟件計算得到的峰值聲功率所對應的頻率與邊界元法和快速多級邊界元法相比，其一致性和關聯(lián)性很好。

4 計算效率

4．1 快速邊界元法計算中存在的問題

在利用LMS virtual lab 10中的快速多級邊界元模塊計算平板、發(fā)動機缸蓋罩和發(fā)動機動力總成輻射噪聲聲功率時，出現(xiàn)了如圖10所示的問題。

圖10 快速多級邊界元法計算輻射噪聲時的問題Fig．10 The problem use fastmultipole boundary elementmethod to calculate radiated sound power

圖10中，用快速多級邊界元計算輻射噪聲聲功率時，不同的計算模型存在一個有效的頻率區(qū)間，這個有效的頻率區(qū)間和模型網格尺寸以及模型的結構有一定的關系。當?shù)陀谶@個頻率范圍時，求解效率會下降，高于這個頻率范圍時，求解精度會下降。

4．2 計算效率對比

計算效率是數(shù)值計算中評價計算方法好壞的一個重要指標，用三種不同的計算方法計算平板、發(fā)動機缸蓋罩和發(fā)動機動力總成輻射噪聲聲功率所用的時間如表1所示。

表1 不同計算方法的計算效率對比Tab．1 Calculation efficiency com parison of different calculatemethods

從表1中可以看出，利用Matlab程序軟件在計算平板、缸蓋罩和發(fā)動機動力總成輻射噪聲聲功率所用的時間都是最少的。平板算例中，模型節(jié)點數(shù)為77，邊界元法和快速多級邊界元法所用的時間分別是Matlab程序軟件的4倍和26倍。在計算發(fā)動機缸蓋罩輻射噪聲聲功率時，模型節(jié)點數(shù)為6 345，邊界元法和快速多級邊界元法所用的時間分別是Matlab程序軟件的25倍和123倍。而在發(fā)動機動力總成算例中，模型節(jié)點數(shù)為23 662，邊界元法和快速多級邊界元法所用的時間分別是Matlab程序軟件的141倍和108倍。

5 結 論

（1）本文以表面速度振動原理為基礎提出了一種簡單的輻射噪聲聲功率計算方法，用Matlab程序編寫了計算輻射噪聲聲功率的軟件，并結合該Matlab程序軟件提出了完整的發(fā)動機零部件及動力總成輻射噪聲計算分析流程。

（2）利用邊界元法、快速多級邊界元法和Matlab程序軟件對平板、發(fā)動機缸蓋罩和發(fā)動機動力總成進行了輻射噪聲聲功率的計算，對比結果顯示Matlab程序軟件計算結果與邊界元法和快速多級邊界元法相比，在低頻段有一定的差距，隨著頻率的增加，差距越來越小，三種方法計算得到的峰值聲功率所對應的頻率點相同。在計算效率上，Matlab程序軟件是邊界元法的25倍和141倍，是快速多級邊界元法的123倍和108倍。

（3）本文中所提出的輻射噪聲聲功率簡單計算方法有效地解決了傳統(tǒng)邊界元存在不適合進行大規(guī)模聲學問題仿真的問題，也解決了快速多級邊界元存在低頻段計算效率低、高頻段精度低的問題，它能滿足各種規(guī)模網格模型輻射噪聲的計算要求，可為發(fā)動機部件及整機輻射噪聲分析和優(yōu)化提供依據(jù)。而且其計算效率比邊界元法和快速多級邊界元法高很多，這樣就提高了工作效率，節(jié)約了計算成本。

［1］張袁元，李舜酩，賀巖松，等．摩托車整車噪聲分析與消聲器改進［J］．振動與沖擊，2011，30（2）：263－266．ZHANG Yuan-yuan，LI Shun-ming，HE Yan-song，et al．Analysis of whole body noise and muffer improvement for a motorcycle［J］．Journal of Vibration and Shock，2011，30（2）：263－266．

［2］舒歌群，劉寧．車輛及發(fā)動機噪聲聲音品質的研究與發(fā)展汽車工程，2002，24（5）：403－407．SHU Ge-qun，LIU Ning．Advance in sound quality research for noise of vehicle and its engine［J］． Automotive Engineering，2002，24（5）：430－407．

［3］岳東鵬，郝志勇，劉月輝，等．柴油機表面輻射噪聲源識別的研究［J］．汽車工程，2004，26（5）：613－616．YUE Dong-peng，HAO Zhi-yong，LIU Yue-hui．A study on surface noise sources identification for diesel engine［J］．Automotive Engineering，2004，5（26）：613－618．

［4］Ide T，Otomori M，Kitamura Y，et al．Improvement of noise performance for automatic transmission using sound pressure optimization［C］．SAE Technical Paper，2010－01－0394．

［5］林逸，馬天飛，姚為民，等．汽車NVH特性研究綜述［J］．汽車工程，2002，24，（3）177－181．LIN Yi，MA Tian-fei，YAOWei-ming，et al．The summary of study on vehicle NVH performance［J］．Automotive Engineering，2002，24，（3）：177－181．

［6］王秀峰．機器輻射聲場分析的新型計算方法研究［D］．合肥：合肥工業(yè)大學，2001．

［7］Chen Z S，Hofstetter G，Mang H B．A symmetric Galerkin formulation of the boundary element method for acoustic radiation and scattering［J］． Journal of Computational Acoustics，1997，5（2）：219－241．

［8］Wu W C．Boundary element acoustics：fundamentals and computer codes［M］．WIT Press，Southampton，UK，2000．

［9］Chew W C，Chao H Y，Cui T J，etal．Fast integral equation solvers in computational electromagnetics of complex structures［J］．Engineering Analysis with Boundary El-ements，2003，27（8）：803823．

［10］戴呈豪，何锃，彭偉才，等．快速多級邊界元進行大規(guī)模數(shù)值計算［J］．機械強度．2008，30（2）：229－234．DAI Cheng-hao， HE Zeng，PENG Wei-cai， et al． Fast multipole boundary elementmethod for large-scale numerical computation［J］．Journal of Mechanical Strength，2008，30（2）：229－234．

［11］Barnes J， Hut P． A hierarchical O （NlogN ） force calculation algorithm［J］．Nature，1986，324：446－449．

［12］龐劍，諶剛，何華．汽車噪聲與振動：理論與應用［M］．北京：北京理工大學出版社，2006．

［13］劉月輝，郝志勇，付魯華，等．車用發(fā)動機表面輻射噪聲的研究［J］．汽車工程，2002，24（3）：213－216．LIU Yue-hui，HAO Zhi-yong，F(xiàn)U Lu-hua，et al．A study on engine radiated noise［J］Automotive Engineering，2002，24（3）：213－216．

［14］馬大猷．噪聲控制學［M］．北京：科學技術出版社．1987．

［15］張俊紅，鄭勇．內燃機振動、噪聲的多體動力學分析［J］．中國機械工程，2006，17（1）：25－28．ZHANG Jun-hong， ZHENG Yong． Multi-body dynamics analysis of internal combustion engine vibration and noise［J］．China Mechanical Engineering，2006，17（1）：25－28．

［16］張重超，楊玉致，朱夢周，等．機電設備噪聲控制工程學［M］．北京：輕工業(yè)出版社，1989．

［17］楊慶佛．內燃機噪聲控制［M］．山西：山西人民出版社，1991．