圓柱繞流噪聲預報的流場與聲場模擬方法對比研究

張翰欽，陳　明，孫國倉

圓柱繞流噪聲預報的流場與聲場模擬方法對比研究

張翰欽，陳明，孫國倉

（武漢第二船舶設計研究所，武漢 430205）

以三維圓柱為研究對象，使用Lighthill聲類比法研究其繞流發聲問題。第一步進行流場計算，分別用大渦模擬（LES）、脫體渦模擬（DES）和瞬態雷諾平均法（URANS）模擬聲源區流場，通過對比流場壓力和渦量等參數，據此選取合適的流場仿真方法；第二步用基于Lighthill方程FW-H積分法和邊界元法預報遠場直發聲，通過和Revell試驗結果比較，分析各種計算方法差別。研究表明：進行流場仿真時LES計算結果最好，IDDES法在保證計算精度條件下能有效減少流場網格數量，URANS法誤差很大；進行輻射噪聲預報時，FW-H積分法和邊界元法基本相同。

聲學；流噪聲；大渦模擬；脫體渦模擬；Lighthill聲類比

圓柱結構作為一種基本的結構形式，在交通運載工程、海洋工程和流體機械等領域都很常見，流體介質在流經圓柱結構后，在一定雷諾數下，圓柱尾流會交替出現脫落渦，這些渦街引起圓柱表面的脈動壓力，從而產生噪聲。圓柱繞流發聲問題很早就有學者作了研究。1977年，Revell通過風洞試驗，以直徑19 mm和38 mm的圓柱為研究對象，測試馬赫數從0.1到0.5之間，雷諾數從45 000到450 000間變化，Revell測量了該雷諾數范圍圓柱繞流的阻力系數及總聲壓級等數據，研究了圓柱繞流阻力和流噪聲之間的關系［1］。該實驗已作為一個基準算例，提供給眾多的研究者作為計算參考。數值仿真方面，針對該問題，Kenneth S.Brentner和Cox等對二維的圓柱繞流進行了仿真計算，所采用的大都是雷諾平均法［2-3］。2002年，Osamu Inoue為了研究均勻流中二維圓柱流噪聲的產生和傳播機理，采用了直接數值模擬的方法，但所計算的雷諾數只有150［4］。隨著計算條件的提高，越來越多的人開始使用大渦模擬（LES）的方法預報流噪聲，2004年，Tang用大渦模擬方法來獲取圓柱繞流水動力數據，再用解FW-H方程的方法預報遠場噪聲［5］。

雖然大渦模擬的計算精度得到了很多學者的認可，但其所需的計算網格數量巨大，對于實尺度艦船噪聲問題仍存在使用限制，而雷諾平均法（URANS）又存在難以準確捕捉湍流的問題。一種混合LES/URANS的方法被提了出來，即脫體渦方法（DES），這種方法在近壁區采用URANS模擬，而在遠離壁面區域采用LES，網格數較小，且能夠保證計算精度。預報輻射噪聲有兩種方法，一種是基于FW-H積分形式方程的解法，另一種用邊界元法求解聲學Helmholtz方程。用積分法求解FW-H方程的優點是它將聲的產生和傳播分別計算，計算量和計算格式要求相對較低，但也存在只能用于計算遠場輻射和不能考慮結構和聲學裝置的影響等缺點。而邊界元方法則在計算流激振動噪聲方面有其獨特的優勢。江文成對比了FW-H積分法和邊界元法在近場和遠場的流噪聲計算上的區別［6］。

為探尋一種保證精度、資源配置合理的計算方法，本文首先對圓柱的低馬赫數下的氣動噪聲進行了數值模擬，采用不可壓縮流動，對比了不同的流場仿真方法和聲場預報方法的區別。分別采用大渦模擬、脫體渦模擬和雷諾平均法進行了圓柱流場仿真，對比了各種方法的差別，從中挑選適用于流噪聲計算的流場仿真方法。流場計算后進行聲場預報，使用了FW-H積分法和邊界元法，對比了兩種方法的差異。通過這些研究工作，本文的計算方法可為大尺度艦船流噪聲預報提供支撐。

1　數學模型及計算方法

流體流動最基本的控制方程為Navier-Stokes方程（NS方程），可表示為如式（1）和式（2）形式

理論上流噪聲分析和流場分析使用同樣的控制方程，但由于聲場和流場在尺度、能量上的差異，用流場控制方程求解聲場物理量會非常復雜。為此Lighthill［7，8］提出聲類比理論，Lighthill沒有做任何簡化及假設，直接對NS方程和連續性方程進行變換，得到非齊次波動方程

2　流場預報

預報氣動噪聲之前首先要對流場進行精確計算，CFD方法一般分為以下幾類：大渦模擬（LES）、脫體渦模擬（DES）和瞬態雷諾平均法（URANS），幾種方法的介紹如表1所示。

表1　流場仿真方法概述

具體的模型設置中，LES的亞格子應力模型采用動態的Smagorinsky-Lilly模型，這種模型可以考慮能量在亞網格上的轉移，從而取得更好的模擬效果。URANS方程對NS方程進行平均，其中雷諾應力項需要用某種湍流模型進行假定，本文采用SST k-ω湍流模型，因為SST k-ω湍流模型采用Standard k-ω湍流模型求解近壁區域，采用k-ε湍流模型求解湍流區域，并且在兩模型之間進行平滑過渡，因此SST k-ω模型較其它模型能更好地模擬近壁區流場，得到的結果也更加準確。DES模擬中除了最基本的DES模型外，還有DDES（Delayed DES）和IDDES（Improved Delayed DES）兩種改進方法，DDES引入過渡函數，推遲RANS到LES的轉換，從而解決了由于RANS到LES提前轉換導致雷諾應力模化不足，產生網格誘導分離的問題；DDES是一種包含壁面LES的方法，它兼具DDES和LES壁面模型的特點，能夠解決邊界層附近“Log-Layer Mismatch”的問題，并且能夠加快分離區RANS到LES的轉換。

計算模型為直徑D=19 mm的圓柱，其展向長度為10倍直徑，即190 mm長。流場入口距圓心為5倍直徑，平行流線方向兩個面距圓心為5倍直徑，出口距圓心為25倍直徑。計算域如圖1所示。邊界條件：face 2、4、6為速度入口邊界條件；face 5為壓力出口邊界條件；face 3為壁面邊界條件，垂直于圓柱的兩個面face 1和face 7設為對稱面。

圖1　幾何模型及邊界條件

入口條件采用速度入口，速度大小為69.19 m/s，流體介質為空氣，馬赫數為0.2，在該馬赫數下可以不考慮流體的可壓縮性。噪聲計算之前，先用CFD計算圓柱繞流流場信息。流場計算過程中，先進行定常計算，定常計算的目的是為非定常計算提供一個穩定的流場，以便流動更好地收斂；再進行非定常計算。控制方程使用基于單元中心的有限體積法離散，其中對流項采用2階迎風差分格式，擴散項采用中心差分格式；壓力與速度的耦合使用PISO算法；離散得到的代數方程使用Guass-Seidel迭代求解。

典型計算域網格如圖2所示，網格壁面進行了加密處理，在軸向進行了200層網格劃分，由于使用不同計算方法時，壁面所需的網格數不同，所以進行了六套網格劃分，劃分的網格第一層網格高度分別為0.006 mm、0.015 mm、0.05 mm、0.09 mm、0.15 mm 和0.3 mm，對應的網格數分別為，760萬、610萬、486萬、410萬、378萬和304萬網格。

圖2　圓柱網格示意圖

首先進行網格驗證計算，流場計算使用URANS、LES、IDDES三種計算方法，將這三種計算方法都應用于六套網格中。提取流場得到的升力系數頻率，該頻率即為周期性漩渦脫落的頻率，再通過公式（4）得到St數；通過總阻力Fd和式（5）計算得到阻力系數Cd，對其取平均值得到-；引入描述網格壁面距離的無量綱參數y+，y+的表達式如式（6）所示，其中Δy表示離壁面的距離，這里指第一層網格高度，uτ是壁面摩擦速度是壁面切應力，μ是流體動力黏性系數，v是流體運動黏性系數，ρ是流體密度。試驗［1］測得St數為0.2、為1.2，對比試驗結果，做出St數和-Cd隨圓柱壁面y+變化的規律分別如圖3和圖4所示。

圖3　St數隨y+變化曲線

圖4　隨y+變化曲線

從圖3和圖4可以看出，使用LES計算時，在y+值小于1時才能得到最佳的計算結果，這是因為LES是基于濾波函數的，無法模擬比自身網格尺度更小的渦，所以如果壁面網格不夠精細的話，LES就無法捕捉足夠的流動，計算結果就會很差。另一方面，使用URANS和DES時，壁面y+值并不是越小越好，而是有一個最佳范圍，由于兩種模型都是使用SST k-ω湍流模型，所以均在y+值等于10左右能夠得到最好的模擬效果。從文獻［9］可知，流場中的壁面區域可劃分為三個范圍，分別是黏性底層、過渡層和對數律層，y+代表與壁面的無量綱距離，當y+＜5時，所對應的區域是黏性底層，此時的速度沿壁面法線方向呈纖細分布，流場為層流狀態，如果在這些區域用湍流來模擬，預報結果將很差。所以使用DES或者URANS模型時，壁面y+值必須超過黏性底層，文獻［9］給出的y+=5作為黏性底層的最外層，本文的計算中發現y+在10左右能取得較好的計算結果。這也說明使用DES進行流場仿真時，所需的網格數量要小于LES模型，如果模擬有許多彎曲壁面的模型，DES所需網格數量將大大減小。

網格驗證計算后，分別采用URANS、DES、DDES、IDDES和LES結合最適應的網格進行流場仿真，將結果與試驗結果對比，如表2所示。

表2　不同流場模擬方法計算的阻力系數和St數的最優結果與試驗值對比（試驗值St=0.2，=1.2）

表2　不同流場模擬方法計算的阻力系數和St數的最優結果與試驗值對比（試驗值St=0.2，=1.2）

流場模擬方法URANS DES DDES IDDES LES -Cd -Cd誤差Sty+ 1.116 1.217 1.198 1.199 1.247 7.00% 1.42% 0.17% 0.08% 3.92% 0.240 0.218 0.221 0.212 0.201 St誤差20.06% 8.91% 10.74% 6.17% 0.70% 9.20 10.09 9.37 9.89 0.59網格數410萬410萬410萬410萬760萬

對比表2中不同湍流模擬方法計算阻力系數和St數，LES方法模擬的渦脫落頻率和試驗值誤差最小，而使用IDDES方法預報的阻力系數值和試驗結果誤差最小。綜合而言，使用LES結合動態Smagorinsky-Lilly模型的方法預報的結果最為理想，頻率和阻力的誤差均在4%以下。IDDES模型和DDES模型模擬的阻力系數和試驗誤差很小，IDDES方法要優于DDES方法，DDES方法又優于DES法。使用URANS模型方法與其他方法相比，誤差較大。

圖5描繪了一個周期內4個階段渦量云圖，對于圓柱繞流問題，當流場雷諾數Re處于65到2×105之間時，尾流中會形成穩定的渦列，渦從圓柱表面脫落并一直延續至下游，在Re=90 000的條件下，形成的尾渦也會很不均勻。對比發現：用LES法得到的流場渦量云圖非常清晰，不僅可以觀察到尾渦交替性的脫落，也可以發現渦核逐漸發展、分離為幾個渦，形成湍流的現象。而用URANS法得到的云圖比較均勻，與實際的尾渦發展成極不均勻的湍流的情況不符。DES在近壁區使用雷諾平均應力方程，在遠場使用大渦模擬，所以在近壁區渦量變化更接近URANS法，而遠離壁面區域接近LES法，從渦量變化上來看，DES在壁面處可以觀察到均勻的渦脫落，而遠離壁面處則更接近LES的渦量場。

圖5　圓柱徑向截面渦量一周期變化云圖

從流場預報來看，LES結合動態Smagorinsky-Lilly模型的方法能最真實、準確地模擬流場變化，但同時它所需的網格數量很大，圓柱表面y+值需要保證在1以下。而使用DES和URANS時，y+值取在10左右最為合適，其中IDDES方法計算的結果較好。所以進行流噪聲預報時，流場預報可以采用LES方法和IDDES方法。

3　聲場預報

流體流經圓柱時，會交替性地產生脫落的漩渦，同時在圓柱面上形成脈動的升力，這種脈動升力就是噪聲形成的主要原因；此外，脫落的漩渦逐漸形成湍流，也會產生四極子噪聲源。目前對流噪聲的預報主要有兩種途徑：一種是基于FW-H積分形式方程的解法，另一種是先從CFD計算結果中獲取聲源場，再使用邊界元法求解聲學Helmholtz方程。兩種方法均是基于Lighthill聲類比理論，在流場模擬準確的基礎上，本文對比兩種方法的差異。

由Focus Williams和Hawkings提出的FW-H方程是一種非齊次的波動方程［10］，它可以通過直接推導流體的連續性方程和Navier-Stokes方程獲得，FW-H方程有如下形式

式中ui代表流體在xi方向的速度分量；un代表流體在 f=0面上的法向速度分量；vi代物面速度在xi方向的速度分量；vn代物面速度在 f=0面上的法向分量；δ（f）代表狄拉克函數；H（f）代表Heaviside函數。如果假設流動是在自由空間中的并且聲源和接受者之間沒有障礙，則方程可以用積分法解析的求解出來，完整的解析解包括面積分和體積分。面積分代表單極子、偶極子和面上四極子聲源的作用，體積分則表示聲源面外部四極子源的作用。本例中流動是亞音速的，體積分可以忽略。時域FW-H積分解可用時域形式求解出來。

邊界元方法同樣是以Lighthill方程為控制方程，對公式（3）進行傅里葉變化，可以得到非齊次的Helmholtz方程，并代入基本解。本例的流場是用不可壓縮流體計算，認為壁面壓力脈動pCFD并沒有計算得到壁面聲壓pa，任一點聲壓的表達式為

其中G（r1，r2）是Helmholtz方程的基本解，pCFD（r2）是流場計算的表面脈動壓力，pa（r2）是表面聲壓。此時，要首先解決邊界上的聲壓問題，才能進一步進行邊界元求解。得到了聲壓的表達式后，用邊界元對邊界面進行網格離散，將型函數Nj代入式（8），最終得到邊界元離散形式的解為

通過式（9），就可以通過先求得壁面的聲壓值，再求解空間任一點的聲壓值。

對于CFD計算得到的流體脈動壓力數據，如果要采用FW-H積分法預報輻射噪聲，則直接打開Fluent的噪聲計算模塊計算即可，方法方便快捷。如果要采用邊界元求解輻射噪聲，還需將流體數據導入聲學邊界元計算程序，并且要劃分邊界元網格，但在用邊界元求解時可以考慮聲源與接受者之間存在障礙的情形，并且可以計算流激振動噪聲。如圖6所示取距圓柱正上方128 D處為監測點，計算值如表3所示，該點處的總聲壓級試驗結果［1］為100 dB，不同流場得到的聲壓曲線對比如圖7所示，積分法和邊界元方法的對比如圖8所示。

表3　監測點總聲壓級與試驗值對比/dB

圖6　監測點示意圖

在表3監測點（0，128 D）總聲壓級的對比中，LES和試驗結果最為接近，IDDES次之，URANS結果最差。從圖7中的LES、IDDES和URANS計算聲壓值以及試驗測量的聲壓值對比發現：基于LES的輻射噪聲預報結果與試驗結果最為接近；基于IDDES的輻射噪聲預報結果和和LES相差不大，只是峰值頻率與試驗結果稍有偏差。URANS法無法捕捉到峰值頻率，且計算得到的聲壓級頻譜曲線非常光順，這是因為URANS方法的局限性，脈動量都是通過湍流模型求得的原因。綜合而言，LES結合動態Smagorinsky-Lilly模型方法預報的流噪聲最為準確；IDDES方法也比較準確，但由于其網格要求并不嚴格，在計算量方面占有優勢；而URANS方法由于將流場的脈動量進行了均化，而代表脈動的雷諾應力項只能通過湍流模型捕捉，所以預報流噪聲結果較差。噪聲計算方法方面，從圖8中可以看出，積分法和邊界元法計算結果基本一致。

圖7　不同流場模型預報輻射噪聲對比

圖8　積分法和邊界元法對比

4　結語

本文分別用大渦模擬（LES）、脫體渦模擬（DES）和瞬態雷諾平均法（URANS）的方法模擬了三維圓柱繞流流場，并分別用FW-H積分法和邊界元法計算了遠場輻射噪聲，對比各算例的計算結果，得到如下結論：

（1）使用LES要保證壁面網格無量綱高度y+值小于1，因為LES是基于濾波函數的方法，y+值大于1會導致網格尺度太大，從而無法模擬足夠精細的流動。使用DES結合SST k-ω湍流模型時，y+取在10左右的范圍較合適，因為y+＜5處于壁面黏性底層，此時流動全為層流，湍流模型不起作用，所以壁面網格并非越精細越好。

（2）流場預報方法中，LES結合動態Smagorinsky-Lilly模型的方法計算的輻射噪聲和試驗結果最為接近，但網格量較大；DES法中的IDDES計算的輻射噪聲結果也較為準確，在保證計算精度的前提下可有效降低網格數量；URANS方法預報結果誤差較大。聲輻射預報方法中，FW-H積分法和邊界元法預報遠場直發聲結果基本相同。

［1］REVELL J D，PRYDZ R A，HAYS A P.Experimental study of aerodynamic noise vs drag relationships for circular cylinders［J］.Aiaa Journal，1977，16（9）：889-897.

［2］BRENTLNER K S，COX J S，RUMSEY C L，et al. Computation of sound generated by flow over a circular cylinderanacousticanalogyapproach［C］//Second ComputationalAeroacoustics（CAA） Workshopon Benchmark Problems.NASA.1997.

［3］COX J S，CHRISTOPHER L.Computation of sound generated by viscous flow over a circular cylinder［C］// Proceedings oftheASME/JSME/IMechE/CSME/IAHR4th International Symposium on Fluid-Structure Interactions，Aeroelasticity，Flow-Induced Vibration&Noise.1997.

［4］INOUE HATAKEYAMA.Sound generation by a twodimensional circularcylinderinauniformflow［J］. Journal of Fluid Mechanics，2002，471（1）：285-314.

［5］TANG K F.Numerical simulation of flow-induced noise by means of the hybrid method with LES and aeroacoustic analogy［D］.University of Siegen，Germany，2004.TANG K F.

［6］江文成，張懷新，孟堃宇.基于邊界元理論求解水下潛艇流噪聲的研究［J］.水動力學研究與進展，2013，28（4）.

［7］LIGHTHILL M J.On sound generated aerodynamically.I. general theory［J］.Proceedings of the Royal Society A Mathematical Physical&Engineering Sciences，1952，211（1107）：564-587.

［8］LIGHTHILL M J.On sound generated aerodynamically. II.turbulence as a source of sound［J］.Royal Society of London Proceedings，1954，222（1148）：1-32.

［9］FLUENT Inc.Fluent User's Guide，Fluent Inc，2003.

［10］WILLIAMS J E F，HAWKINGS D L.Sound generation byturbulenceandsurfacesinarbitrarymotion［J］. Philosophical Transactions of the Royal Society A Mathematical Physical&Engineering Sciences，1969，264（1151）：321-342.

ComparativeStudyontheFlowFieldandAcousticFieldSimulation forNoisePredictionInducedbytheFlowaroundaCylinder

ZHANG Han-qin，CHENMing，SUN Guo-cang
（Wuhan 2nd Ship Design and Research Institute，Wuhan 430205，China）

Lighthill acoustic analogy method is used to predict the noise induced by the flow around a three dimensional cylinder.Firstly，the flow field of sound source area is predicted by large-eddy simulation（LES），detached-eddy simulation（DES）and Transient Reynolds-average method（TRAM）respectively.The proper methods to simulate the flow field are selected via comparing the contours of vorticity and pressure.Secondly，the FW-H integral equation method and BEM method based on Lighthill acoustic analogy equation are used to predict the far-field noise，difference of these methods is analyzed by comparing their results with Revell's experimental data.It is shown that the flow field simulation result of LES has a best agreement with the experimental data，IDDES method can effectively reduce the number of grid of the flow field with the precision guaranteed.TRAM can yield large errors.FW-H integral equation method and BEM method can essentially yield the same results in predicting radiation noise.

acoustics；flow-induced noise；LES；DES；Lighthill acoustics analogy

O422.6

ADOI編碼：10.3969/j.issn.1006-1335.2016.03.006

1006-1355（2016）03-0026-06

2015-11-17

國家自然科學基金資助項目（51409199）

張翰欽（1992-），男，江西省南昌市人，碩士生，主要研究方向為水下結構物流噪聲。E-mail：hackinzhq@sina.com