陳武 , ，韓斐 , ，周毅 ,

（1.中國空氣動力研究與發(fā)展中心 氣動噪聲控制重點實驗室，四川 綿陽 621000；2.南京理工大學 能源與動力工程學院，南京 210094）

流體流經(jīng)單個或多個障礙物的現(xiàn)象在航空航天領(lǐng)域中處處可見（如飛機起落架及艙體[1]、直升機旋翼[2]等）。柱體繞流是典型的繞流問題之一，由于存在尾流間的相互作用及不同尾流與柱體間的相互作用，相比較單個柱體產(chǎn)生的尾流場，多柱體后形成的尾流場的演化特性更加復雜，研究內(nèi)涵會更加豐富[3]。

氣動噪聲是在湍流運動中，伴隨流體的壓力脈動，通過流體的彈性和慣性作用向遠離流動區(qū)域空間傳播而形成的流動噪聲[4]。目前多柱體繞流流場特性和氣動噪聲特性[5]的研究在航空航天領(lǐng)域中具有重大的研究意義。國內(nèi)外學者對多柱體繞流及其氣動噪聲問題進行了一系列的研究[6-16]。Ishigai 等[6]對3 種不同排列的多圓柱組合進行了實驗研究，結(jié)果表明渦形成區(qū)域和康達效應(yīng)（Coanda effect）對流場特性具有重要影響。Auger 等[7]通過實驗研究發(fā)現(xiàn)，當相鄰圓柱之間的間距不同時，流體流經(jīng)并列圓柱后產(chǎn)生的旋渦脫落有明顯的變化。Zheng 等[8]對3 個并排方柱后的流動進行了數(shù)值仿真研究，確定了流場中存在5 種不同的流態(tài)，并詳細研究了每種流態(tài)對應(yīng)的湍流特性。劉歡[9]采用有限體積方法對低雷諾數(shù)下等邊布置三方柱繞流進行了數(shù)值研究，預測了臨界間距比隨雷諾數(shù)的增加而增大。Xu 等[10]對3 ～ 4 個并排的柔性圓柱體進行了實驗研究，探討了流致振動問題。Inoue 等[11-12]采用二維非定?？蓧嚎sNavier-Stokes 方程組的直接求解方法，研究了在低馬赫數(shù)下的均勻來流場中并列雙方柱和并列三方柱所產(chǎn)生的氣動噪聲特性。杜炳鑫等[13]提出了一種基于流聲分解法的繞流噪聲預測方法，對低雷諾數(shù)（Re=200）下的串列和并列雙圓柱繞流噪聲問題進行了研究。高威等[14]采用改進型延遲分離渦（Improved delayed detached eddy simulation, IDDES）模型對串列雙圓柱繞流進行數(shù)值仿真，分析了不同來流速度、圓柱間距比和圓柱直徑下串列雙圓柱繞流的氣動噪聲特性。余雷等[15]采用基于非線性k-ε模型的限制數(shù)值尺度（Limited numerical scales, LNS）方法對串列雙圓柱繞流算例進行了研究。Zhou 等[16]采用直接數(shù)值模擬方法對高間距比的雙方柱繞流進行研究，并將近場區(qū)域非 Kolmogorov 理論的-5/3能譜冪律函數(shù)與極端減速事件相聯(lián)系。

上述的研究大多集中在相同尺度柱體后的尾流及其氣動噪聲問題。近年來，具有多尺度/分形輪廓的物體（分形格柵[17-21]、不規(guī)則端板[22-23]和多尺度邊緣的葉輪[24-25]等）產(chǎn)生的尾流引起學者廣泛的研究。部分學者認為，分形的多尺度邊界在某種程度上可以調(diào)控湍流特性[26]。

在航空航天工業(yè)領(lǐng)域，分形流動具有相當廣闊的應(yīng)用前景和潛力[27-29]。航空航天面臨的主要環(huán)境問題之一是機體噪聲，而擾流板是產(chǎn)生機體噪聲的主要部件之一。Nedi?等[27]研究發(fā)現(xiàn)機翼裝載大孔隙的分形擾流板，能夠在保持與常規(guī)擾流板類似空氣動力學性能的同時將其誘導的氣動噪聲聲壓級降低4 dB。Nedi?等[28]對NACA0012 機翼常規(guī)后緣和非平坦鋸齒/分形后緣產(chǎn)生的渦脫落特性及其氣動性能進行了實驗對比研究。實驗結(jié)果表明多尺度/分形機翼后緣能夠有效改善空氣動力學性能。Gehlert 等[29]研究了分形端板對NACA0012 半展翼翼尖渦的影響。結(jié)果發(fā)現(xiàn)，分形端板的引入對渦的幾何結(jié)構(gòu)和湍流特性都有很大的影響。Laizet 等[26]證明了多尺度/分形湍流場中存在所謂的空間尺度展開（Space-scale unfolding, SSU）機制，能夠顯著提高標量混合效率。Laizet 等[19]首次對二維多尺度方柱陣流場進行了數(shù)值研究，揭示了分形排列的影響機制。Baj 等[30]用粒子圖像測速技術(shù)和三重分解技術(shù)，對多尺度方柱陣的湍流近場進行了實驗研究。Tao 等[31]采用直接數(shù)值模擬的方法研究了兩種具有相同阻塞率的方柱布置（分形布置和規(guī)則布置），并使用快照本征正交分解方法定量的評估了兩種湍流中相干結(jié)構(gòu)特征的差異。

綜上所述，前人的研究集中在對分形結(jié)構(gòu)產(chǎn)生的旋渦脫落特性上，對于分形結(jié)構(gòu)產(chǎn)生的尾流氣動噪聲特性的研究稍顯不足。本文采用大渦模擬結(jié)合K-FWH 聲比擬方法對分形布置方柱尾流場和聲場進行數(shù)值模擬研究。此外，對相同阻塞率的規(guī)則布置方柱尾流場進行了對比仿真模擬，以比較其與分形布置方柱尾流場特性的異同。

1 數(shù)值方法及計算幾何模型

1.1 數(shù)值方法

本文利用壁面自適應(yīng)局部渦粘 （Wall adapting local eddy viscosity, WALE） 亞格子模型開展大渦模擬，該模型表達式[32]為

式中Ls=min（kd,CωΔ1/3），k為馮卡門常數(shù)，d為距離壁面最近的網(wǎng)格尺度，Cω為WALE 常數(shù)。

WALE 模型是在湍流結(jié)構(gòu)的運動和動力學性質(zhì)的基礎(chǔ)上，將轉(zhuǎn)動張量包含在模型中構(gòu)造而成。該模型的優(yōu)點是能夠正確的反映出近壁區(qū)域渦粘系數(shù)與垂直壁面距離的三次方成正比的特性，且模型中不含任何有關(guān)邊界幾何尺寸的參數(shù)，在復雜的湍流模擬中得到廣泛的應(yīng)用。

本研究對于聲場的求解采用K-FWH 聲壓時域解公式[33]，并根據(jù)實際情況簡化公式得到瞬時聲壓p'的時域解，其公式為

本文所采用的K-FWH 方程聲壓時域解公式由FWH 方程聲壓時域解結(jié)合Kirchhoff 方法得到，該公式不需要求解物體表面壓力的法向?qū)?shù)，簡化了聲壓的計算。

1.2 分形及規(guī)則布置

本文采用兩種不同排列方式的并列多方柱布置（規(guī)則布置方柱和分形布置方柱），研究在多方柱后形成的尾流特性及其氣動噪聲特性。在規(guī)則布置中，4 個相同的方柱體并排布置，與來流方向垂直（從下至上依次為Bar1、Bar2、Bar3和Bar4）。在分形布置方柱中，3 種分形尺度的方柱體并排布置，與來流方向垂直（從下至上依次為Bar1、Bar2、Bar3、···、Bar7）。

圖1 為兩種布置方柱在笛卡爾坐標系上的示意圖，其中左邊為前視圖，右邊為側(cè)視圖，虛線表示計算域邊界，具體的幾何細節(jié)和計算參數(shù)見表1。笛卡爾坐標系的原點設(shè)置在Y-Z平面的中心，所有柱體均垂直于X-Y平面。在兩種布置中，并排的方柱體放置在XBar/TBar=8 處，XBar表示計算域入口到柱體中心的流向距離，TBar表示柱體流向厚度?；诜街飨蚝穸萒Bar和入口來流速度為3.75 m/s。長度LX、LY、LZ分別表示沿流向、法向和展向的計算域長度，兩種布置的計算域大小相同。對于規(guī)則和分形布置，兩個相鄰柱體表面之間的垂直間距Lg/TBar分別為2.33 和1.33。在規(guī)則布置中，每個方柱法向厚度為TBar。在分形布置中，采用三次分形迭代的方柱，大尺寸、中尺寸和小尺寸柱體的法向厚度迭代比為T0∶T1∶T2= 4∶2∶1（見圖1b））。

表1 規(guī)則與分形模型具體幾何參數(shù)Tab.1 Computational and geometrical parameters of regular and multiscale models

圖1 兩種布置方柱的示意圖Fig.1 Two kinds of square column layout diagram

兩種布置具有相同的阻塞率（即σ= 0.30），阻塞率表達式為

式中：Af為方柱橫截面面積之和；At為計算域橫截面面積（即LY×LZ）。例如，規(guī)則布置方柱的阻塞率σ= 4TBar/LY。

在劃分網(wǎng)格時采用指數(shù)函數(shù)對柱體壁面附近網(wǎng)格進行加密，使壁面網(wǎng)格第一層網(wǎng)格Y+≈1。遠離壁面的部分采用均勻加密，最后得到規(guī)則布置的網(wǎng)格數(shù)為6.00 × 106，分形布置為6.96 × 106。展向設(shè)置為周期性邊界，方柱表面為無滑移壁面邊界，其余邊界為遠場邊界?？臻g離散為二階迎風格式，時間離散為二階隱式格式。時間步長設(shè)置為Δt=5 × 10-6s，在流場充分發(fā)展之后，對流場結(jié)果和噪聲結(jié)果進行統(tǒng)計，取樣統(tǒng)計時間為7 × 104個時間步長（對應(yīng)總時間為0.35 s）。

2 數(shù)值方法驗證

在研究規(guī)則和分形布置方柱引起的湍流特性和氣動噪聲特性之前，需要驗證本文所使用的數(shù)值模擬方法的準確性。因此，本文對單方柱繞流進行數(shù)值模擬計算用以驗證流場計算方法的可行性。與此同時，由于單方柱繞流氣動噪聲的實驗結(jié)果較少，國內(nèi)外學者對于串列雙圓柱繞流氣動噪聲的研究較多，本文采用串列雙圓柱模型驗證噪聲計算方法的正確性。

2.1 流場仿真方法驗證

本文研究對象為規(guī)則與分形布置多柱體繞流，由于多柱體繞流缺乏相對應(yīng)的實驗或仿真數(shù)據(jù)，無法直接將數(shù)值結(jié)果與前人結(jié)果進行直接比較。單方柱是多方柱構(gòu)成的基本元素，多方柱尾流場是多個單方柱尾流相互作用形成的。因此，本小節(jié)通過單方柱流場特性與前人結(jié)果進行對比，進而驗證所采用流場仿真方法的正確性。

本文研究單方柱的邊長為L，入口來流速度為Uin，基于方柱邊長和來流速度的雷諾數(shù)為Rein=UinL/v= 104。

選取方柱柱體中心作為坐標系原點，計算域為：-10L≤LX≤20L，-9L≤LY≤9L，0≤LZ≤4L。對方柱周圍進行加密處理，壁面第一層網(wǎng)格Y+≈2，X-Y平面上的計算域網(wǎng)格如圖2 所示，總網(wǎng)格數(shù)為2.36 ×106，單方柱網(wǎng)格加密方式與分形及規(guī)則布置方柱的網(wǎng)格加密方式相同。單方柱繞流的邊界條件設(shè)置與規(guī)則及分形布置方柱繞流相同。時間步長為Δt=5 ×10-6s，取樣統(tǒng)計時間為70 000 個時間步長（總時間為0.35 s）。

圖2 X-Y 平面上的單方柱計算網(wǎng)格Fig.2 Computational grids in the X-Y plane for a single square cylinder

圖3 為無量綱平均流向速度u1/Uin在中心線上的流向演化，文獻[31,35]為數(shù)值結(jié)果，文獻[36-39]為實驗結(jié)果。從圖3 中可以看到流向速度沿著中心線緩慢下降，靠近壁面時驟然下降直至壁面流向速度為0，在壁面后方產(chǎn)生回流區(qū)，隨后速度突然上升至約0.7 倍的入口流向速度并趨于平緩。在柱體后方上游區(qū)域，流向速度u1/Uin與前人實驗及模擬的結(jié)果吻合較好，在方柱后方下游區(qū)域，本文結(jié)果與Tao 等[31]（Rein=1 200）以及Portela 等[35]（Rein=3 900）的數(shù)值模擬結(jié)果類似，結(jié)果落在前人結(jié)果組成的合理區(qū)間范圍內(nèi)。

圖3 無量綱平均流向速度在中心線的流向演化Fig.3 Streamwise evolution of the normalized mean streamwise velocity u1/Uin along the centerline

為了進一步驗證本文的數(shù)值結(jié)果是否可以準確模擬單方柱后方的尾流，本文還與進一步比較驗證了流場中的單點二階統(tǒng)計量。圖4 給出了中心線上3 個方向的無量綱均方根速度其中仿真結(jié)果與Portela 等[35]的結(jié)果吻合。本文的數(shù)值結(jié)果相比于Tao 等[31]（Rein=1 200）的直接數(shù)值模擬結(jié)果，兩者結(jié)果在柱體前方區(qū)域類似。在方柱后方下游區(qū)域，本文的在峰值處略高于Tao 等[31]的直接數(shù)值模擬結(jié)果，但在峰值過后結(jié)果吻合良好，在峰值過后略低于Tao 等[31]的直接數(shù)值模擬結(jié)果，在峰值過后低于Tao等[31]的直接數(shù)值模擬結(jié)果，但與Portela 等[35]的結(jié)果較為接近，在合理區(qū)間范圍內(nèi)。由于本文使用的是大渦模擬方法，而Tao 等[31]使用的是直接數(shù)值模擬方法，直接數(shù)值模擬與大渦模擬相比對于空間分辨率需求很高，因而計算量大、耗時多，對計算資源消耗較高，對計算機內(nèi)存依賴性強，精度會更為準確，但本文的數(shù)值結(jié)果與Tao 等[31]也較為接近，可見本文的計算結(jié)果在保證數(shù)值模擬結(jié)果準確性的同時可節(jié)約一定的計算資源。

圖4 均方根速度沿中心線的流向演化Fig.4 Streamwise evolution of all three rms velocities along the centerline

方柱后部渦脫落會在方柱表面引起周期性壓力波動變化，進而誘導遠場氣動噪聲。本文計算了表征流動周期性的相似準則：斯特勞哈爾數(shù)St=fL/Uin= 0.122。在當前入口雷諾數(shù)下Rein= 104，斯特勞哈爾數(shù)計算結(jié)果接近文獻[36]研究結(jié)果 （St≈0.13）。

圖5 為不同雷諾數(shù)下的平均阻力系數(shù)CD，文獻[31, 35, 40-41]表示數(shù)值結(jié)果，文獻[42-45]表示實驗結(jié)果。圖5 表明，雷諾數(shù)Rein=104時，平均阻力系數(shù)CD落在前人結(jié)果的合理區(qū)間范圍內(nèi)，進一步驗證了本文計算結(jié)果的正確性。

圖5 不同雷諾數(shù)Rein 下的平均阻力系數(shù)CDFig.5 Mean drag coefficient CD at different Rein

此外，本文計算了方柱表面的平均壓力系數(shù)CP分布，即

式中：P為當?shù)貕毫Γ籔in為入口壓力； ρ為流體密度；Uin為入口速度。圖6 為單方柱A-B、B-C、C-D 和D-A 這4 個面的平均壓力系數(shù)CP分布，A-B、C-D面位垂直于流向的面，B-C、D-A 面位平行于流向的面。文獻[35, 41]表示數(shù)值結(jié)果，文獻[46-48]表示實驗結(jié)果。圖6 表明平均阻力系數(shù)CP結(jié)果與前人的研究結(jié)果一致，進一步證實了本文所使用的流場計算方法可以很好的預測固體表面的渦脫落現(xiàn)象。

圖6 單方柱4 個面上的平均壓力系數(shù)CP 分布Fig.6 Distribution of the mean pressure coefficient CP over 4 sides of a single square cylinder

圖7 為不同雷諾數(shù)Rein下升/阻力系數(shù)均方根值之比，其中文獻[35, 40-41]表示數(shù)值結(jié)果，文獻[43]表示實驗結(jié)果。從圖7 中可以看出，本文結(jié)果十分接近Portela 等[35]的預測結(jié)果，說明達到捕捉表面壓力脈動的精度要求。

圖7 不同雷諾數(shù)Rein 下的平均阻力系數(shù)與平均升力系數(shù)均方根值之比Fig.7 The ratio of the mean drag coefficient to the mean lift coefficient rms at different Rein

本小節(jié)計算了單方柱繞流的各項速度及壓力指標以驗證所采用數(shù)值方法的正確性，首先計算了流向平均速度及3 個方向的均方根速度在中心線上的流向演化，其次計算了斯特勞哈爾數(shù)、平均阻力系數(shù)、升/阻力系數(shù)均方根值之比以及單方柱表面的壓力系數(shù)分布，驗證結(jié)果顯示仿真結(jié)果與前人的實驗及數(shù)值結(jié)果相吻合，證實了本文所使用流場計算方法的可靠性。

2.2 噪聲仿真方法驗證

本文采用簡化后的時域解公式（式（2））對聲場進行計算。在本課題組前期研究[49]中對串列雙圓柱繞流的噪聲進行了數(shù)值仿真，使用了3 套不同精細程度的網(wǎng)格計算了串列雙圓柱繞流在某一觀測點處的噪聲頻譜分布，3 個算例所對應(yīng)的網(wǎng)格數(shù)分別為238 萬、665 萬和478 萬，計算結(jié)果與Langley 研究中心[50]在靜流設(shè)施（Quiet flow facility, QFF）以及基礎(chǔ)空氣動力學研究（Basic aerodynamic research,BART）風洞中的實驗結(jié)果相吻合，如圖8 所示，其中算例3 是在算例1、2 基礎(chǔ)上對網(wǎng)格所做相應(yīng)改進，吻合度較高。本文采用的噪聲仿真方法與前期串列雙圓柱繞流氣動噪聲研究相同，保證了本文所使用聲比擬方法的可靠性。

圖8 觀測點處的聲壓級功率譜密度Fig.8 Power spectral density of sound pressure level at observation points

3 計算結(jié)果與分析

3.1 湍流統(tǒng)計特性

在本節(jié)中將討論規(guī)則和分形布置方柱繞流的湍流場特性，包括二維瞬時渦量、速度梯度張量第二不變量Q值瞬時等值面以及X-Y平面的平均流向速度、均方根速度和平均壓力。

圖9 為在任意選定展向位置Z處的無量綱瞬時渦量場，其中紅色表示高渦量值≥5，藍色表示低渦量值≤0.1。在圖9 中可以清楚的看到對于兩種情況的尾流，在方柱后方上游區(qū)域均呈現(xiàn)出高度間歇性。流體空間間歇性是指流體在空間中的狀態(tài)呈現(xiàn)非均勻性，主要包括速度、壓力、渦量等的非均勻性，而隨著尾流向下游區(qū)域發(fā)展，由于各尾流間相互作用，渦量場逐漸變得均勻，空間間歇性減小。

圖9 X-Y 平面上的瞬時渦量等值面云圖Fig.9 Contour plots of the magnitude of instantaneous vorticity in the X-Y plane

圖9 的瞬時渦量場為二維瞬時云圖，為了更加清晰地說明規(guī)則和分形布置的尾流渦結(jié)構(gòu)，在此利用速度梯度張量第二不變量Q準則對流場進行可視化處理。圖10 為與圖9 相同時刻下的Q值瞬時等值面圖，其中=1，顏色代表無量綱瞬時流向速度u1/Uin。圖10a）表明在規(guī)則布置方柱流場中，相鄰柱體后方的渦呈反相位脫落；圖10b）表明在分形布置方柱流場中，渦脫落行為較為不規(guī)則。

圖10 Q 值瞬時等值面圖（ =1）Fig.10 Instantaneous isosurfaces of Q-value (=1)

除了尾流的瞬時渦結(jié)構(gòu)外，本文研究了平均流向速度場和均方根速度場。圖11 為X-Y平面上的無量綱平均流向速度場u1/Uin云圖。隨著向下游發(fā)展，規(guī)則布置方柱流場趨于均勻分布，而分形布置方柱流場中仍能觀察到最大尺寸柱體后方的大尺度尾流的影響。

圖11 無量綱平均流向速度場u1/Uin云圖Fig.11 Contour plots of the normalized mean streamwise velocity field (u1/Uin)

圖12 為X-Y平面上的無量綱均方根速度場云圖。從圖12a）中可以看出規(guī)則布置每個方柱后產(chǎn)生的均方根速度分布相似，圖12b）顯示分形布置3 種尺度的方柱后產(chǎn)生不同的均方根速度分布，最大尺度方柱后部的湍流強度最大。

圖12 無量綱均方根速度云圖Fig.12 Contour plots of the normalized RMS velocity ()

圖13 為X-Y平面上的無量綱化平均壓力場（P-Pin）/（）云圖。圖13a）顯示規(guī)則布置中各方柱后部的低壓區(qū)分布相似，隨著向下游發(fā)展，各方柱后方壓力恢復速度相當。圖13b）表明不同尺度方柱后的壓力場分布不同，最大尺寸柱體后的壓力恢復速度最慢。上述發(fā)現(xiàn)與平均流向速度場分布的結(jié)果相呼應(yīng)。

圖13 無量綱平均壓力場（P-Pin）/（）云圖Fig.13 Contour plots of the normalized mean pressure field(P-Pin)/()

3.2 法向分布特性

圖10 表明規(guī)則布置和分形布置方柱的平均速度流場有較大區(qū)別，為了進一步研究兩者速度場之間的異同，本文給出了無量綱平均流向速度和法向速度在選定的5 個流向位置（XBar/TBar=2, 5, 10, 15,20）處的法向分布，見圖14 和圖15。由圖14 可得，兩種布置方柱在上游區(qū)域同時存在著類似射流（u1/Uin＞1）和類似尾流（u1/Uin＜1）的流場特性。值得關(guān)注的是，在下游遠場區(qū)域，規(guī)則布置方柱的流向速度分布有趨于均勻分布的趨勢，而分形布置方柱在下游區(qū)域仍能顯示出明顯的尾流分布特征。

圖14 無量綱化平均流向速度u1/Uin 的法向分布Fig.14 Vertical distribution of the normalized mean stream- wise velocity u1/Uin

從圖15 中可以看出，在尾流的上游區(qū)域（XBar/TBar=2），兩種布置方柱都表現(xiàn)出強烈的法向運動。在上游區(qū)域（XBar/TBar=2），柱體后方的尾流存在著正負交替的法向速度分布，促使流向速度的法向分布趨于平滑（見圖13）。在XBar/TBar≥5 之后，由于流速在Y方向的平均梯度（即法向運動驅(qū)動力）變?。ㄒ妶D14），導致兩種布置方柱后方尾流的法向運動變得不明顯。

圖16 無量綱速度均方根值（）的法向分布Fig.16 Vertical distribution of the normalized mean squares velocity ()

綜上所述，在上游區(qū)域（XBar/TBar=2），兩類流場均呈現(xiàn)出較強的間歇非均勻性。在下游區(qū)域（XBar/TBar≥15），分形布置流場中仍然能觀察到大尺寸柱體對流向速度的影響，而該區(qū)域的法向速度和均方根速度在法向上趨于均勻分布。

3.3 旋渦脫落統(tǒng)計特性

多柱體繞流場中的重要特征是存在渦脫落以及渦之間的相互作用。本小節(jié)將對相關(guān)的重要問題（例如柱體的渦脫落頻率、渦脫落相位關(guān)系、相位鎖定的渦脫落行為以及分形布置對于規(guī)則布置能否改變相干渦結(jié)構(gòu)等）進行討論。

圖17 為規(guī)則布置方柱的升/阻力系數(shù)時間演化曲線。圖中升力系數(shù)CL和阻力系數(shù)CD的初始時刻（t=0）相同且為任意選定。圖17a）表明在規(guī)則布置方柱中4 個方柱的渦脫落周期相同，各柱體之間存在旋渦脫落鎖定，即相鄰柱體的渦呈反相位脫落，其峰值大小同樣保持一致。也就是說，對于任意的一個方柱，當它的升力系數(shù)CL為最大值時，其相鄰方柱的升力系數(shù)為最小值。從圖17b）中可以看出在規(guī)則布置方柱中，相鄰方柱的阻力系數(shù)隨時間演化似乎是鎖定的，即阻力系數(shù)呈同相位變化。

圖17 規(guī)則布置4 個方柱的升力系數(shù)CL 和阻力系數(shù)CD時程曲線Fig.17 Time trace curves of the lift coefficient CL and the drag coefficient CD behind four cylinders for the regular array

為了定量的評估渦脫落行為，本文對7 × 104個時間步（對應(yīng)規(guī)則布置方柱約48 個渦脫落周期）的升/阻力系數(shù)時間演化數(shù)據(jù)進行分析，通過快速傅里葉變換，得到相應(yīng)的功率譜密度。由于在規(guī)則布置方柱中，4 個方柱的升/阻力系數(shù)變化相似，渦脫落周期相同。因此，在圖18 中只給出了一個任意選定方柱的升/阻力系數(shù)的功率譜密度，兩條垂直虛線為fTBar/Uin=0.176 和0.352。圖18a）為規(guī)則布置方柱的升力系數(shù)CL頻譜，主峰對應(yīng)的斯特勞哈爾數(shù)為St=0.176。相比較于單方柱繞流（St=0.122），規(guī)則布置中的尾流相互作用提升了旋渦脫落頻率。圖18b）為規(guī)則布置方柱的阻力系數(shù)CD頻譜，主峰對應(yīng)的斯特勞哈爾數(shù)為St=0.352，對應(yīng)升力系數(shù)斯特勞哈爾數(shù)的兩倍。

圖18 規(guī)則布置方柱中選定方柱的升力系數(shù)CL 和阻力系數(shù)CD 功率譜密度Fig.18 Power spectral density of the lift coefficient CL and drag coefficient CD of the chosen square cylinders for the regular array

Zheng 等[8]研究了3 個并列方柱的升力系數(shù)頻譜，其3 個柱體的旋渦脫落頻率分別為St=0.172,0.183, 0.172。由于存在旋渦脫落的相位滯后，導致其頻譜中出現(xiàn)一個表征更大尺度運動的次級峰（St=0.011）。然而，本研究中規(guī)則布置各相鄰方柱之間的升力系數(shù)相位鎖定，因此在頻譜上未觀察到顯著的次級峰現(xiàn)象。

圖18 為多尺度分形布置方柱的升力系數(shù)CL和阻力系數(shù)CD的時間演化曲線。由于分形布置方柱Y軸中心對稱，故選取了對應(yīng)于圖1b）中分形布置方柱的下半?yún)^(qū)域柱體（即Bar1、Bar2、Bar3和Bar4）的統(tǒng)計結(jié)果進行展示。

從圖19 中可以看出，分形布置方柱后方的旋渦脫落相關(guān)性較低，不同尺度柱體的渦脫落周期有明顯差異。同樣對分形布置方柱的渦脫落行為進行定量評估，將7 × 104個時間步的升/阻力系數(shù)時間演化數(shù)據(jù)進行快速傅里葉變化，得到分形布置方柱的功率譜密度。

圖19 分形布置方柱升力系數(shù)CL 和阻力系數(shù)CD 時程曲線Fig.19 Time trace curves of the lift coefficient CL and the drag coefficient CD of square cylinders for the multiscale array

圖20 為分形布置方柱中選定方柱的升/阻力系數(shù)頻譜，圖20a）中3 條虛線為相應(yīng)柱體的渦脫落頻率，分別為fTBar/Uin=0.141, 0.192, 0.793，圖20b）中3 條虛線表示兩倍的柱體渦脫落頻率，分別為fTBar/Uin=0.282, 0.384, 1.586。從圖20a）中可以看出，多尺度的分形布置方柱產(chǎn)生的渦脫落頻率不同，大中小3 個尺度柱體的渦脫落頻率分別為fTBar/Uin=0.141, 0.192, 0.793。在升力系數(shù)頻譜中，由于相鄰柱體尾流間的影響，各柱體的渦脫落過程會導致在相鄰柱體的頻譜中出現(xiàn)對應(yīng)的峰值。圖20b）顯示阻力系數(shù)頻譜在相應(yīng)柱體的兩倍渦脫落頻率處未存在明顯峰值。上述發(fā)現(xiàn)表明，分形布置方柱產(chǎn)生的尾流間相互作用會對柱體阻力產(chǎn)生顯著的影響，尤其是最小尺度柱體。

圖20 分形布置方柱中選定方柱的升力系數(shù)CL 和阻力系數(shù)CD 功率譜密度Fig.20 Power spectra density of the lift coefficient CL and the drag coefficient CD of the chosen square cylinder in the multiscale case

圖21 為規(guī)則布置和分形布置各方柱體的平均阻力系數(shù)CD。根據(jù)面積權(quán)重計算可得規(guī)則布置和分形布置方柱的平均阻力系數(shù)分別為3.475 和3.677。因此，在雷諾數(shù)Rein=104下，兩者的平均阻力系數(shù)大致相等，分形結(jié)構(gòu)并不會明顯改變平均阻力系數(shù)。

圖21 規(guī)則布置和分形布置各方柱體的平均阻力系數(shù)CDFig.21 Mean drag coefficient CD of different square cylinders for the regular array and the multiscale array

3.4 氣動噪聲特性

圖22 給出了規(guī)則布置方柱和分形布置方柱位于r/TBar=32 處的遠場噪聲聲壓級（Sound pressure level, SPL）指向分布圖，其中r是與觀測點相同的X-Y平面上離坐標原點的距離。聲壓級指向性可由各點瞬時聲壓的均方根計算得到。

圖22 遠場聲壓級指向分布圖（r/TBar=32）Fig.22 Directivity of SPL at r/TBar=32

圖22 顯示了規(guī)則布置方柱呈現(xiàn)偶極子特性，且與來流方向平行的位置處聲壓級最大，與來流垂直方向的位置處聲壓級最小。而對于分形布置方柱，其指向性分布大致與規(guī)則布置方柱相同，也呈現(xiàn)偶極子分布特性，但聲壓級最小值略偏于來流垂直方向?？傮w而言，分形結(jié)構(gòu)對遠場噪聲的指向性分布影響不大。

本文選取坐標點（0, 32TBar, 0）作為觀測點計算遠場氣動噪聲特性。圖23 為規(guī)則布置和分形布置方柱在選定觀測點處不同柱體引起的瞬時聲壓隨時間演化曲線。從圖23a）中可以看出，在規(guī)則布置方柱中，Bar1～ Bar4的時程曲線變化與升力系數(shù)變化規(guī)律相似；圖23b）表明分形布置方柱中Bar1～ Bar7產(chǎn)生的噪聲則混亂無序。

圖23 規(guī)則與分形布置方柱在選定觀測點（0, 32TBar, 0）處不同柱體引起的瞬時聲壓隨時間變化曲線Fig.23 Time variation curves of the instantaneous acoustic pressure of different square cylinders at the observation point(0, 32TBar, 0) for the regular array and the multiscale array

圖24 為規(guī)則布置和分形布置方柱在選定觀測點（0, 32TBar, 0）處的噪聲功率譜密度（Power spectral density, PSD），其中垂直虛線分別表示fTBar/Uin=0.141, 0.176, 0.192 和0.793。從圖24 中可以看出，規(guī)則布置方柱的噪聲功率頻譜在低頻段只有一個峰值（St=0.176），而分形布置方柱分別在高、中、低區(qū)段產(chǎn)生3 個峰值（St=0.141, 0.192 和0.793）。低頻峰值為大尺寸柱體產(chǎn)生，中頻峰值為中尺寸的柱體產(chǎn)生，高頻峰值為最小尺寸柱體產(chǎn)生（見圖20）。這是由于規(guī)則布置方柱只存在單一尺度的旋渦脫落，而分形布置存在不同尺度的旋渦及其相互作用。

圖24 規(guī)則布置和分形布置方柱在選定觀測點（0, 32TBar, 0）處的聲壓級功率譜密度Fig.24 Power spectral density of the sound pressure level at the receiver point (0, 32TBar, 0) for the regular array and the multiscale array

上述頻譜分析表明，規(guī)則及分形布置方柱的聲壓級功率譜密度與對應(yīng)的升力系數(shù)頻譜密切相關(guān)。圖24 表明分形布置方柱的高頻能量大于規(guī)則布置，說明分形布置能夠?qū)⒌皖l聲壓轉(zhuǎn)化為高頻聲壓。

4 結(jié)論

1） 通過數(shù)值模擬得到單方柱繞流流場特性和串列雙圓柱繞流氣動噪聲特性，并與前人的實驗及數(shù)值結(jié)果進行對比，驗證了本文采用的大渦模擬結(jié)合K-FWH 聲類比方法的準確性。結(jié)果表明本文所采用的數(shù)值仿真方法能夠較好的預測柱體繞流問題的流場特性與聲場特性。

2） 規(guī)則布置方柱尾流場旋渦脫落呈現(xiàn)明顯的“相位鎖定”現(xiàn)象，而分形布置尾流場旋渦脫落混亂無序，兩種流場的平均阻力系數(shù)大致相等。由于分形布置不同尺度柱體產(chǎn)生的旋渦脫落頻率不同，相鄰柱體尾流間存在互相影響，從而會在升/阻力系數(shù)頻譜上產(chǎn)生次級峰。

3） 在選定觀測點處，規(guī)則布置不同柱體的瞬時聲壓值呈現(xiàn)“相位鎖定”現(xiàn)象，而分形布置不同柱體的瞬時聲壓值相對混亂。分形布置方柱由于存在不同尺度的旋渦脫落及其相互作用，在聲壓級功率頻譜上產(chǎn)生了3 個不同的主頻峰值。

4） 分形布置和規(guī)則布置方柱的聲壓級指向性分布均呈現(xiàn)偶極子分布特性，兩者在各方向的指向分布大致相同。分形布置能夠改變聲壓級功率譜密度分布，將低頻聲壓向高頻聲壓轉(zhuǎn)移。