瞿 堯 劉滟鈺 王 寧 柴應(yīng)彬 楊 平 張?jiān)侜t*

(武漢理工大學(xué)船海與能源動(dòng)力工程學(xué)院1) 武漢 430061) (星云智熵科技有限公司2) 北京 100020)

0 引 言

在流體力學(xué)和流體聲學(xué)的研究過(guò)程中，鈍體繞流問(wèn)題一直是重要的研究課題之一.隨著大型客機(jī)的迅速發(fā)展，越來(lái)越多學(xué)者對(duì)飛機(jī)噪聲進(jìn)行了研究.Chow等[1]對(duì)大型客機(jī)A340進(jìn)行了實(shí)驗(yàn)研究.實(shí)驗(yàn)結(jié)果表明:起落架的噪聲比襟翼噪聲要高6 dB，圓柱繞流噪聲對(duì)飛機(jī)具有很大的影響.相比于單圓柱繞流，雙圓柱繞流更符合實(shí)際情況，現(xiàn)實(shí)中更多的是流體在多個(gè)物體之間的相互作用下流動(dòng).研究雙圓柱繞流產(chǎn)生的噪聲中，上游產(chǎn)生的渦會(huì)作用在下游圓柱面上產(chǎn)生更小的渦，下游產(chǎn)生的渦脫落具有更為復(fù)雜的流動(dòng)特性，這對(duì)復(fù)雜的氣動(dòng)噪聲研究具有重要價(jià)值.

美國(guó)航天局(NASA)研究中心[2]在空氣動(dòng)力學(xué)基礎(chǔ)研究風(fēng)洞(BART)和靜流設(shè)備(QFF)風(fēng)洞中進(jìn)行了大量的空氣動(dòng)力學(xué)和氣動(dòng)噪聲試驗(yàn).國(guó)內(nèi)外學(xué)者參照這些試驗(yàn)的模型做了很多研究，Lockard等[3]利用三維求解器 CFL3D模擬計(jì)算了串列雙圓柱繞流，利用 FW-H 方程計(jì)算得到遠(yuǎn)場(chǎng)噪聲，得到了與試驗(yàn)結(jié)果吻合的計(jì)算結(jié)果.Papaioannou等[4]比較了雙圓柱的三維模擬同二維模擬的區(qū)別，結(jié)果表明，三維模擬能更好地模擬尾跡的渦量場(chǎng)、圓柱表面受力及渦脫落頻率，得到了與實(shí)驗(yàn)更接近的臨界間距比.趙超等[5]基于不同RANS模型的DES方法研究了雙圓柱繞流表面壓力系數(shù)均方根值的對(duì)比，得出SST-Kω模型對(duì)比試驗(yàn)有較好的吻合.馬瑞軒等[6]基于SST-SAS自適應(yīng)尺度模擬結(jié)合FW-H方法計(jì)算了雙圓柱繞流流動(dòng)特性和聲場(chǎng)特性.趙威等[7]建立與試驗(yàn)?zāi)Ｐ拖喈?dāng)?shù)恼归L(zhǎng)并利用Lighthill聲類比方法求出遠(yuǎn)場(chǎng)聲壓.劉國(guó)慶[8]通過(guò)Lighthill聲類比方法得到很好的仿真結(jié)果，并對(duì)圓柱的間距、排列雷諾數(shù)對(duì)噪聲的影響作了分析.

基于此，文中通過(guò)LES方法計(jì)算標(biāo)準(zhǔn)雙圓柱繞流模型，利用基于有限元求解變分形式的聲類比方程和基于面積分求解聲類比方程計(jì)算遠(yuǎn)場(chǎng)噪聲，分別采用聲學(xué)計(jì)算軟件ACTRAN和流場(chǎng)計(jì)算軟件Fluent中的聲學(xué)模塊計(jì)算，并分析兩種方法存在的差異性.

1 數(shù)值求解理論

1.1 大渦模擬

大渦模擬(large eddy simulation，LES)是最近幾十年來(lái)興起的重要數(shù)值模擬方法，它主要通過(guò)濾波函數(shù)處理Navier-Stokes(N-S)方程，過(guò)濾掉較小尺度的渦，然后精準(zhǔn)計(jì)算該尺度以上所有湍流的運(yùn)動(dòng).

濾波器公式定義為

(1)

式中：D為流體域；G為決定大渦尺度的濾波函數(shù).

將空間區(qū)域離散后得到有限控制體積.

(2)

式中：V為控制體積

濾波器函數(shù)G(x,x′)取為

(3)

將N-S方程經(jīng)過(guò)濾波函數(shù)過(guò)濾過(guò)后，得到如下方程.

(4)

式中：τij為亞格子應(yīng)力張量，定義如下.

(5)

亞格子應(yīng)力張量τij的計(jì)算需要建立亞格子應(yīng)力模型計(jì)算，基于Boussinesq假定，引入了大尺度應(yīng)變率張量Sij和亞格子黏性系數(shù)υsgs，計(jì)算公式為

(6)

式中：υsgs為小渦的黏性系數(shù)，而大尺度應(yīng)變張量Sij的計(jì)算公式為

(7)

目前常用的亞格子應(yīng)力模型有Smagorinsky模型、WALE模型、動(dòng)態(tài)Smagorinsky模型.本文使用更容易收斂的Smagorinsky亞格子應(yīng)力模型.

1.2 聲類比方程

面積分聲類比方程是通過(guò)引入廣義函數(shù)，通過(guò)N-S方程整理得到的，也是聲比擬方程中最常用的形式，具體形式為

(8)

式中：

ρ0為遠(yuǎn)場(chǎng)流體密度；p′為遠(yuǎn)場(chǎng)聲壓；c0為遠(yuǎn)場(chǎng)聲速；f為積分表面；ui為xi方向上的速度分量；un為f=0面上的法向速度；vn是物面速度的法向分量；δ(f)為狄拉克函數(shù)；H(f)為赫維賽德階躍函數(shù)；Tij為L(zhǎng)ighthill聲類比應(yīng)力張量；pij為表面載荷.式(8)右邊第一項(xiàng)為單極子聲源，第二項(xiàng)為偶極子聲源，第三項(xiàng)為四極子聲源.

2 計(jì)算模型介紹

2.1 流場(chǎng)計(jì)算模型

采用NASA在消聲室做的QFF模型試驗(yàn)，該試驗(yàn)下風(fēng)速Uo=43.4 m/s.試驗(yàn)中，兩個(gè)相同直徑的圓柱沿著風(fēng)口依次串聯(lián)排列，其中圓柱直徑D=0.057 15 m，雷諾數(shù)計(jì)算為Re=1.66×105，兩個(gè)圓柱之間的中心距離為L(zhǎng)z=16D.選取3D作為模型展長(zhǎng),然后利用文獻(xiàn)[8]的修正方案修正結(jié)果.為了充分捕捉繞流過(guò)后的聲源，使結(jié)果更加的準(zhǔn)確.取上游段長(zhǎng)度為L(zhǎng)1=5D，下游段長(zhǎng)度L2=16.3D，整個(gè)外流域尺寸24D×11D×3D，建模時(shí)圓心取在上游圓柱底部圓心處.

由于計(jì)算的模型為高雷諾數(shù)下圓柱繞流，圓柱尾部會(huì)產(chǎn)生大量不同尺度的渦，為了充分捕捉不同尺度的渦來(lái)提高計(jì)算精度，對(duì)圓柱周圍及圓柱尾流部分進(jìn)行加密處理.采取結(jié)構(gòu)化網(wǎng)格劃分，網(wǎng)格總數(shù)為350萬(wàn)，圖1為俯視平面的網(wǎng)格圖.

圖1 網(wǎng)格示意圖

流場(chǎng)入口邊界條件為速度入口，出口邊界條件為壓力出口，上下壁面為周期性邊界，圓柱壁面為無(wú)滑移壁面.流場(chǎng)先進(jìn)行定常計(jì)算，采用的是標(biāo)準(zhǔn)k-ε模型，然后取時(shí)間步長(zhǎng)Δt=0.001 s，計(jì)算一定步數(shù)后，待模型收斂效果很好時(shí)，進(jìn)行LES模型計(jì)算，Δt=2×10-5s.待LES模型計(jì)算收斂時(shí)，提取非定常流場(chǎng)計(jì)算結(jié)果，時(shí)間步數(shù)設(shè)為25 000步，計(jì)算時(shí)間即0.5 s，每隔25步保存A個(gè)結(jié)果，得到1 000個(gè)采樣數(shù)據(jù).

2.2 模型介紹

基于變分形式聲類比方程的ACTRAN聲學(xué)模型中，整個(gè)計(jì)算網(wǎng)格域分為聲源域、傳播域、無(wú)限域三個(gè)部分，見(jiàn)圖2.

圖2 雙圓柱繞流聲學(xué)計(jì)算模型

ACTRAN聲學(xué)軟件當(dāng)中，為了更直觀的觀察聲學(xué)信號(hào)，通過(guò)傅里葉變換將時(shí)域下的瞬態(tài)流場(chǎng)信息轉(zhuǎn)化為頻域下的流場(chǎng)信息.其中，轉(zhuǎn)化后的分辨率Δf=1/(N·Δt)，聲學(xué)計(jì)算得最高頻率fmax=1/(2n·Δt).所以，計(jì)算所得的聲學(xué)信號(hào)分辨率為2 Hz，最高聲學(xué)頻率為1 000 Hz.

而基于面積分聲類比方程直接計(jì)算寬頻帶噪聲，由于計(jì)算的低馬赫數(shù)下圓柱繞流問(wèn)題，忽略四極子聲源，只考慮壁面上的偶極子聲源.因此在聲學(xué)模塊中選擇圓柱壁面為噪聲源，數(shù)據(jù)保存和ACTRAN一致.

為了驗(yàn)證流場(chǎng)計(jì)算結(jié)果的準(zhǔn)確性，見(jiàn)圖3.

圖3 壓力系數(shù)測(cè)點(diǎn)位置示意圖

取圓柱中間截面，以圓柱中心為坐標(biāo)原點(diǎn)，圓柱駐點(diǎn)為起始點(diǎn)，順時(shí)針?lè)较蛎扛?0°取一個(gè)監(jiān)測(cè)點(diǎn).計(jì)算上、下游圓柱表面的壓力系數(shù)和壓力均方根系數(shù)，與試驗(yàn)結(jié)果進(jìn)行對(duì)比.其中壓力系數(shù)計(jì)算公式為

(9)

式中:p0為參考聲壓；ρ0為流體密度；u0為流速.

三維建模原點(diǎn)位置在上游圓柱底部中心處，遠(yuǎn)場(chǎng)噪聲監(jiān)測(cè)點(diǎn)以(9.11D，-2.4D，1.5D)為圓心，布置A(-8.33D，27.715D，1.5D)、B(9.11D，32.49D，1.5D)、C(26.55D，27.815D，1.5D)三個(gè)監(jiān)測(cè)點(diǎn)，將仿真模擬計(jì)算得到的聲學(xué)結(jié)果與試驗(yàn)比較.

3 流場(chǎng)結(jié)果分析

3.1 壓力系數(shù)

圖4為上、下游圓柱表面各測(cè)點(diǎn)的壓力系數(shù).整體具有很好的對(duì)稱性，壓力系數(shù)從駐點(diǎn)角度開(kāi)始向兩邊先下降到最低點(diǎn)，然后上升到流動(dòng)分離點(diǎn)時(shí)壓力系數(shù)趨于平緩，平緩的這一部分圓柱表面的剪切力趨于零.

圖4 上、下游圓柱表面各點(diǎn)壓力系數(shù)

表1為分離點(diǎn)角度對(duì)比表，由表1可知，上、下游的壓力系數(shù)與試驗(yàn)值吻合的很好，渦脫落分離點(diǎn)都與試驗(yàn)接近，最大分離點(diǎn)誤差為3.7%.說(shuō)明了流場(chǎng)計(jì)算結(jié)果可靠，可以作為后續(xù)計(jì)算聲場(chǎng)的流場(chǎng)結(jié)果.

表1 分離點(diǎn)角度對(duì)比表

3.2 流場(chǎng)特性分析

渦脫落是圓柱繞流的重要特征之一，而雙圓柱繞流的渦脫落比單圓柱繞流的渦脫落要復(fù)雜的多，圖5為雙圓柱繞流X-Y截面和X-Z截面速度云圖.

圖5 截面云圖

由圖5可知，兩個(gè)圓柱下游部分存在著不同尺度的渦，由渦聲理論可知，漩渦的產(chǎn)生、拉伸、變形、和破碎都會(huì)產(chǎn)生噪聲，這也是繞流的主要噪聲來(lái)源，結(jié)合漩渦結(jié)構(gòu)能夠很好的知道圓柱繞流的發(fā)聲機(jī)理.為了更加形象的看出繞流產(chǎn)生的脫落渦，圖6為雙圓柱繞流的速度渦量云圖.

圖6 速度渦量云圖

由圖6可知，上游圓柱存在著周期性脫落的漩渦，上游脫落渦抨擊下游圓柱表面產(chǎn)生了大量的破碎渦，同時(shí)自身也產(chǎn)生不同尺度的脫落渦，這部分是整個(gè)流場(chǎng)發(fā)聲的主要部分.

4 聲場(chǎng)結(jié)果分析

由于本文計(jì)算的展長(zhǎng)為3D，而試驗(yàn)計(jì)算的展長(zhǎng)為16D，所以需要引入一個(gè)修正系數(shù)來(lái)減小由于展長(zhǎng)差異造成的影響.采用文獻(xiàn)[8]的修正方法，具體形式為

當(dāng)LC≤LS

SPL=SPLS+10 lg(L/LS)

(10)

當(dāng)LC≤LS≤L

SPL=SPLS+20 lg(LC/LS)+10lg(L/LC)

(11)

當(dāng)L≤LC

SPL=SPLS+20lg (L/LS)

(12)

式中：SPLS為數(shù)值計(jì)算得到的聲壓級(jí)；L為實(shí)驗(yàn)的展向長(zhǎng)度；LS為數(shù)值計(jì)算得展向長(zhǎng)度；LC為展向相關(guān)尺度.

修正中所用到的展向相關(guān)尺寸是不同頻率和不同位置的函數(shù)，因此不同頻率下的修正值都應(yīng)該是不相同的.鑒于目前展向修正方法并不太完善，文中采用相同的相關(guān)尺度進(jìn)行修正，認(rèn)為所有頻率下LC=4D，則計(jì)算得結(jié)果修正值為9.77 dB，得到的聲壓級(jí)頻譜圖見(jiàn)圖7.

圖7 三個(gè)監(jiān)測(cè)點(diǎn)頻域聲壓級(jí)

由圖7可知，在100～500 Hz以內(nèi)結(jié)果與試驗(yàn)結(jié)果吻合度比較好，500 Hz以上和100 Hz以下的差異是由于試驗(yàn)風(fēng)洞上下是由兩塊平板固定圓柱，能量會(huì)被剛性平板反射出去，而流場(chǎng)計(jì)算時(shí)上下壁面采用的是周期性邊界.

表2～3為三個(gè)監(jiān)測(cè)點(diǎn)仿真的峰值頻率和峰值聲壓級(jí)，與實(shí)驗(yàn)結(jié)果進(jìn)行了比較.

表2 ACTRAN計(jì)算的峰值聲壓對(duì)比結(jié)果

表3 Fluent聲學(xué)模塊計(jì)算的峰值聲壓對(duì)比結(jié)果

由表2～3可知，峰值聲壓級(jí)和頻率與試驗(yàn)量級(jí)相差很小，二者之間峰值頻率的差異主要是由于數(shù)值模擬中無(wú)法模擬實(shí)驗(yàn)中轉(zhuǎn)捩帶讓渦脫落提前這一現(xiàn)象造成的.

5 結(jié) 束 語(yǔ)

文中分別采用聲學(xué)軟件ACTRAN和計(jì)算流體動(dòng)力學(xué)軟件Fluent中的聲學(xué)模塊，即采用基于有限元求解變分形式聲類比方程和基于面積分求解聲類比方程兩種計(jì)算方法，計(jì)算了雙圓柱繞流的流動(dòng)特征和聲場(chǎng)特征，并且和NASA的QFF試驗(yàn)進(jìn)行了比較.對(duì)于峰值頻率和聲壓級(jí)，經(jīng)過(guò)修正的ACTRAN計(jì)算更接近于實(shí)驗(yàn)值，峰值聲壓級(jí)誤差為1.1%，峰值頻率誤差為4.5%.而Fluent聲學(xué)模塊計(jì)算的聲壓級(jí)誤差為2.2%，峰值頻率誤差為7.4%.兩種計(jì)算方法采用了相同的流場(chǎng)結(jié)果，只是噪聲計(jì)算方法不同，但都能比較好的與試驗(yàn)結(jié)果進(jìn)行吻合，這表明了計(jì)算方法的合理性.對(duì)軍事領(lǐng)域最大聲壓級(jí)要求精度較高的可以采用ACTRAN軟件，而對(duì)于工況眾多的工程問(wèn)題，可以采用更為方便的Fluent聲學(xué)模塊.