基于轉捩模型及聲比擬方法的高精度圓柱分離渦/渦致噪聲模擬?

王光學 王圣業(yè) 葛明明 鄧小剛

1)(國防科技大學空天科學學院,長沙 410073)

2)(中山大學物理學院,廣州 510275)

(2017年12月17日收到；2018年8月10日收到修改稿)

1 引 言

圓柱尾跡渦及渦致噪聲問題在航空航天、風工程等實際工程中非常具有代表性.如何準確預測渦街特性、渦街產聲大小一直是計算流體力學(CFD)及計算氣動聲學(CAA)的熱點.尤其是在亞臨界雷諾數(shù)范圍內(1×103—2×105),圓柱邊界層為層流而尾跡已轉變?yōu)橥牧鳒u,對CFD中的湍流模型是個巨大的挑戰(zhàn)[1].近年來,隨著計算機技術的飛速發(fā)展,一些簡單的學術問題通過采用大渦模擬(LES)方法已經(jīng)能得到很好的結果[2?5].然而,對于如大型客機起落架、大型通風管道等雷諾數(shù)較高并且特征尺度較大的工程外形,邊界層內最大湍流尺度相對于幾何尺度也會變得很小,使計算花費巨增.因此,基于雷諾平均Navier-Stokes(RANS)方程的湍流模型仍然是工業(yè)應用CFD/CAA的中堅力量.

目前工程應用中,轉捩主要依靠經(jīng)驗或半經(jīng)驗公式確定.其中,Menter和Langtry[6]提出的基于當?shù)仃P聯(lián)的γ-Reθ轉捩模型,憑借其與現(xiàn)代CFD方法良好的兼容性,在工程領域得到了廣泛的應用.然而,當轉捩后出現(xiàn)大規(guī)模非定常分離時,γ-Reθ模型仍受到傳統(tǒng)RANS方法的局限,在分離區(qū)會高估渦黏性而無法準確預測氣動力.另一方面,傳統(tǒng)的尺度模擬方法,如分離渦模擬(DES)[7]、尺度自適應模擬(SAS)等[8],基于全湍流模型,對于高雷諾數(shù)的大分離問題能很好地求解,而對于中等雷諾數(shù)(1×104—3×105)時邊界層附近存在轉捩的問題卻無法準確預測.

2006年,Langtry等[9]便針對該問題指出:“眾所周知RANS湍流模型在大規(guī)模分離區(qū)無法準確模擬；而混合RANS/LES方法,如DES或SAS,在準確捕捉分離流動的物理特性方面可能是必要的.”其后,諸多學者為實現(xiàn)該目的,在結合γ-Reθ模型和尺度求解方法方面做了許多工作.2011年,S?rensen等[10]提出了一種結合DES和γ-Reθ模型的分離渦模擬方法,在雷諾數(shù)10—1×106大范圍內的圓柱繞流問題中均能得到與實驗符合的平均阻力結果.2013年,You和Kwon[11]發(fā)展出了基于γ-Reθ模型的SAS方法,也能在圓柱繞流問題中很好地預測平均氣動力.2014年,Qiao等[12]通過轉捩平板算例重點關注了DES方法產生的?；瘧p耗(MSD)問題,表明該問題會阻礙邊界層轉捩的發(fā)生,而采用延遲分離函數(shù)可使該現(xiàn)象得到改善.2017年,Hodara和Smith[1]又通過結合γ-Reθ模型, 將Sánchez-Rocha和Menon[13]的Hybrid LDKM方法(一種加權函數(shù)型混合RANS/LES方法)擴展至轉捩分離流動的模擬中,在NACA 63-415機翼和圓柱繞流的氣動力預測方面取得了良好的效果.然而,上述學者對模擬結果的討論主要集中于宏觀氣動力,對速度梯度、壓力脈動以及聲壓級等更為精細的流場變化關注較少.同時,梯度和脈動量的準確模擬更依賴于低耗散的高精度數(shù)值格式,上述工作也未能涉及.

本文采用CFD和聲比擬相結合的方法,在近場采用Tran-DDES方法結合七階WCNS-E8T7格式對非定常湍流流場進行高保真求解；在遠場使用基于Ffowcs Williams-Hawkings(FW-H)方程的積分法模擬觀測區(qū)噪聲.算例考慮亞臨界雷諾數(shù)(約4×104)的圓柱分離/噪聲問題,及圓柱尾跡中放置翼型后的流動干擾問題.

2 湍流及噪聲模型

2.1 SST-γ-Reθ轉捩模型

γ-Reθ轉捩模型通過經(jīng)驗關聯(lián)函數(shù)控制邊界層內間歇因子的生成,再通過間歇因子控制湍流模型中湍流的生成.其無量綱守恒形式的輸運方程為:

其中源項及系數(shù)具體參見文獻[6].該模型包含三個關鍵的經(jīng)驗關聯(lián)函數(shù):轉捩動量厚度雷諾數(shù)Reθt、轉捩區(qū)長度Flength和臨界動量厚度雷諾數(shù)Reθc,并對使用的CFD軟件平臺的差異較為敏感.本文作者在前期工作中已依據(jù)T3系列低速平板實驗對經(jīng)驗關聯(lián)函數(shù)進行了標定,并在低速問題中開展了應用[14].

γ-Reθ轉捩模型通常需要耦合k-ω剪切應力輸運(SST)模型構成四方程模型進行求解.無量綱SST模型輸運方程的守恒形式為:

其中源項及系數(shù)具體參見文獻[6].γ-Reθ模型和SST模型的結合主要是通過間歇函數(shù)γ來修正k方程的生成項和破壞項:

其中F1,orig是原SST模型中的混合函數(shù)；γfi是γ-Reθ模型經(jīng)過分離誘導轉捩修正后的間歇函數(shù)；γD,min和F3參見文獻[6],ω方程不做修正.

2.2 Tran-DDES方法

分離渦模擬方法由Spalart等[15]提出,將LES方法和SA模型結合,并通過比較當?shù)鼐W(wǎng)格尺度與壁面距離實現(xiàn)兩種方法的自動切換.其后,Strelets[16]借鑒Spalart的思想,通過比較當?shù)鼐W(wǎng)格尺度與湍流長度尺度,將DES方法引入SST模型.然而,早期DES方法的限制器在復雜網(wǎng)格上處理RANS和LES的轉換過程中過于生硬,會造成模化應力損耗等問題.2006年,Spalart等[17]借鑒Menter的SST模型構造思想,采用“延遲LES函數(shù)”改善了原版本的MSD問題,發(fā)展出了DDES方法.由于γ-Reθ模型是耦合SST模型求解的,而DDES方法僅是針對k方程耗散項進行的修正,因此轉捩-分離渦模擬(Tran-DDES)方法是直接對(5)式的耗散項修正為:

這里DDES限制器為

轉捩限制器和分離渦模擬限制器均作用于k方程耗散項上,兩者之間可能會產生干擾,需要引起關注.要避免DDES判斷在轉捩區(qū)開啟,保證在層流-轉捩區(qū)中經(jīng)驗關聯(lián)函數(shù)不被DDES修正干擾,本文采用了如下保護函數(shù):

而fd是原DDES方法中的延遲函數(shù),

這里,ν和νt分別為運動黏性系數(shù)和渦黏系數(shù),κ為馮.卡門常數(shù)0.41；dw為壁面距離.

2.3 Ffowcs Williams-Hawkings聲比擬方法

遠場噪聲采用FW-H方程對近場計算得到的非定常流場參數(shù)進行積分得到.聲源面可采用固體邊界面或可穿透聲源面兩種.本文采用固體邊界面,并且假設聲源面包含了所有聲源信息時的情況下,時域求解的FW-H方程[18]為

其中聲源面用f(xi,t)=0表示,r表示觀測點距聲源點的距離,Li和H分別代表載荷噪聲(偶極子聲源)和厚度噪聲(單級子聲源).tadv代表推遲時間,即觀測點接收到聲波的時刻.

由于計算網(wǎng)格量限制,計算模型的展向長度往往與實驗模型存在差別.因此需要根據(jù)聲源信息的相干性對結果進行校正.較為常用的校正公式為Kato公式[19]:

其中,(Spp(f))exp為校正后聲壓；(Spp(f))sim為數(shù)值積分得到的待校正聲壓；Lsim,LC和Lexp分別代表計算模型展向長度、流動展向相關長度和實驗模型展向長度.

3 高精度數(shù)值格式

加權緊致非線性格式(WCNS)首先由Deng和Zhang[20]提出.之后,一系列WCNS格式被發(fā)展[21,22]并廣泛用于流動問題中[23,24].本文采用的WCNS-E8T7格式是由Liu等[22]通過在經(jīng)典的5階顯式加權插值模板基礎上添加緊致項構造而來.該格式在保持較短的插值模板長度的基礎上,實現(xiàn)了精度及色散、耗散特性的顯著提高.同時,在遇到激波或數(shù)值不穩(wěn)定情況時,通過設計加權系數(shù)使緊致項關閉,7階緊致插值格式又可退回至5階顯式插值,從而提高格式在較復雜外形求解時的魯棒性.

具體的,7階加權緊致插值公式為

權ωk由下式給出,

這里,

其中β0,β1和β2是光滑因子,和5階顯式加權插值相同.β5被設計用來控制緊致項,具體參考文獻[22].

Wang等[25]在渦輸運、三角翼等典型渦流動算例中,對WCNS-E8T7格式進了驗證.結果表明該格式能夠在較粗糙的網(wǎng)格上達到甚至超過傳統(tǒng)二階格式在極密網(wǎng)格下對渦結構的分辨率,即計算效率更高.Seo等[2],Boudet等[3]以及Jiang等[5]在模擬圓柱或圓柱-翼型算例時均采用了高階精度的空間離散格式,以更好地避免格式色散、耗散誤差對流場/聲場的影響.因此為更方便地對比湍流模型的差異,本文空間離散同樣基于高階精度格式.為本文研究的轉捩后渦街結構的保持、渦致噪聲的準確模擬奠定了基礎.另外需要指出的是,在計算網(wǎng)格導數(shù)及雅克比時,本文采用了對稱守恒網(wǎng)格導數(shù)方法(SCMM)[26],以保證WCNS格式在曲線網(wǎng)格求解時的數(shù)值精度.

4 計算結果

4.1 單圓柱轉捩分離

單圓柱算例考慮ReD=4.3×104雷諾數(shù)情況,此時圓柱壁面附近的層流邊界層開始分離,剪切層隨后發(fā)生轉捩.流場結果參考Szepessy和Bearman[27]在ReD=4.3×104,圓柱展長Lz=2.7D的實驗結果和Seo等[2]在ReD=4.6×104,圓柱展長Lz=3.0D的不可壓縮LES結果.

本算例采用7階緊致格式WCNS-E8T7,同時添加準3階格式MUSCL(κ= ?1/3)[28]作為參考.湍流模型選擇Tran-DDES并和傳統(tǒng)SSTDDES進行對比.計算網(wǎng)格采用O型拓撲,包含180×180×30網(wǎng)格單元,與文獻[2]中網(wǎng)格設置保持一致.時間推進采用雙時間步方法,100個左右對流周期(T=D/U)后,擬周期的渦街流動建立.此時,統(tǒng)計平均過程開始并維持200T.

圖1給出了120T后升阻力系數(shù)隨時間的變化過程,對比圖2中Seo等的LES結果,Tran-DDES方法得到升阻力系數(shù)波動幅值與LES更為接近.表1給出了氣動力的統(tǒng)計結果,升力系數(shù)隨卡門渦街脫落周期振蕩,實驗測得的斯特勞哈爾數(shù)為St=0.19,WCNS結合Tran-DDES與之符合很好,而結合傳統(tǒng)DES則表現(xiàn)稍差.另一方面,觀察MUSCL結合Tran-DDES結果,其明顯偏離了實驗值.這是由低階格式較大的色散誤差造成的,并且會對后續(xù)聲學結果的處理造成影響.CL,rms和CD,rms分別為升、阻力系數(shù)的均方根值,反映擬周期振蕩的振幅大小.MUSCL和WCNS結合Tran-DDES均能與實驗符合,比SST-DDES表現(xiàn)得更好.

圖1 單圓柱算例升阻力系數(shù)隨時間變化Fig.1.Time variations of the drag and lift cofiicients for the single cylinder.

圖2 Seo等的單圓柱LES結果[2]Fig.2.The LES results obtained by Seo et al.using LES[2].

本文中采用的網(wǎng)格較密,在平均流場特性方面,WCNS與MUSCL結果差異不大.下面主要對比湍流模型之間的差異.圖3對比了WCNS-E8T7格式得到的表面壓力分布,右圖給出了Seo等的LES結果作為參考.Tran-DDES在背壓區(qū)包括吸力峰部分均與實驗值符合得很好,而SST-DDES的在背風區(qū)負壓值則偏小.觀察圖4中平均流場的渦強分布,采用全湍流模型計算時,圓柱兩側的湍流剪切層抵抗失穩(wěn)能力更強.即SST-DDES預測的回流區(qū)長度更長,在背風面附近流動的動能也更大,因此靜壓比Tran-DDES得到的要小.雖然在圓柱表面,湍流邊界層的摩擦阻力大于層流狀態(tài),但湍流尾跡引起的壓差阻力減小占主導,因而造成總阻力偏小.

表1 單圓柱算例氣動力統(tǒng)計結果Table 1.Statistical results of aerodynamic cofiicients for the single cylinder.

圖3 單圓柱表面壓力分布,右圖為Seo等的LES結果[2]Fig.3.Mean pressure cofiicient around the cylinder,right:results by Seo et al.using LES[2].

圖4 單圓柱平均流場的渦強分布Fig.4.Distributions of mean vorticity intensity for the single cylinder.

然后分析遠場觀測點噪聲結果,參考Jacob等[29]在ReD=4.8×104,圓柱展長Lz=30D的實驗數(shù)據(jù).采用FW-H聲比擬方法模擬遠場觀測點的噪聲大小,觀測點位于圓柱垂直于來流的軸線上R=185D處.由圓柱繞流聲場的輻射指向特性可知,垂直于來流的軸線為偶極子軸,并且聲壓級最大.由于實驗模型的展長為Lz=30D,而計算模型的展長僅為Lz=3D.因此本文采用了Kato公式,根據(jù)聲源信息的相干性對結果進行校正.圖5給出了該觀測點處功率譜密度計算結果與實驗的對比.隨著卡門渦街的周期性脫落,聲壓級在能譜圖上隨頻率(St數(shù))形成了主峰、二次峰和三次峰等.主峰為圓柱壁面交替的正負壓力脈動產生的偶極子聲源,Jacob實驗測得的主峰處St數(shù)接近0.2,相較Seo實驗偏大,這與兩次實驗的來流雷諾數(shù)及模型差異有關.另一方面,由于FW-H積分公式中時間導數(shù)通常僅采用二階差分格式,使觀測點處聲壓級譜上的St數(shù)均略低于氣動力統(tǒng)計得到的St數(shù)(表1),也部分導致了WCNS結合SST-DDES的主頻更接近實驗值.然而,對比主峰處聲壓幅值,WCNS結合Tran-DDES的結果優(yōu)于SST-DDES.MUSCL結合Tran-DDES得到的主頻和幅值均偏離實驗稍遠,這與表1中St數(shù)趨勢相符合.

圖5 圓柱上方R=185D處聲壓功率譜密度結果對比Fig.5.Comparison of PSD of acoustic pressure at R=185D directly above the cylinder.

4.2 圓柱-翼型干擾流動

單圓柱算例中,主要研究了轉捩模型對圓柱壁面附近非定常流場的影響.可以看到,是否考慮轉捩對預測的渦街形態(tài)有很大影響.下面進一步研究圓柱繞流與下游翼型的相互干擾.如圖6,圓柱的直徑為0.01 m,基于圓柱直徑的雷諾數(shù)和單圓柱算例近似,為ReD=4.8×104.下游翼型為NACA0012,弦長c=0.1 m.翼型與圓柱的距離為0.1 m.以翼型前緣為原點,遠場噪聲觀測點Pfar的坐標為(0.05 m,1.85 m).

圖6 圓柱-翼型實驗設置及觀測點位置Fig.6.Setup of the rod-airfoil experiment and location of measurements of interest.

本算例采用7階緊致格式WCNS-E8T7,湍流模型選擇Tran-DDES和SST-DDES進行對比.計算網(wǎng)格采用H-O型拓撲,圓柱表面仍為180個網(wǎng)格單元,翼型表面為320個網(wǎng)格單元.展向拉伸0.03 m,均布30個網(wǎng)格單元,總網(wǎng)格單元數(shù)約300萬.文獻[3]采用傳統(tǒng)LES方法模擬,圓柱和翼型表面分別布置了200,350個網(wǎng)格單元,展向布置30個網(wǎng)格單元.但由于法向方向網(wǎng)格量稍高,當前網(wǎng)格的總單元量反而高于文獻[3].計算采用雙時間步方法,時間步長取1.5×10?6s,平均統(tǒng)計時長0.03 s.結果參考Jacob等[29]在里昂中央理工大學(Ecole Centrale de Lyon)的亞聲速無回聲風洞的實驗結果,展向長度0.3 m.

首先分析氣動結果. 圖7展示了0.03 s時,Tran-DDES模型和SST-DDES模型結合高精度WCNS-E8T7得到的Q等值面分布.圓柱兩端剪切層失穩(wěn)形成卡門渦街,尾跡渦打到機翼并破裂.特別是渦街破裂后的較小尺度湍流結構也能被清晰地分辨.

以翼型頂點為原點,觀察平均流場中四個典型剖面的速度分布(圖8).x/c= ?0.87和x/c=?0.255為圓柱尾跡中兩個典型站位,平均速度在圓柱流向軸線處達到最小值.在x/c=?0.87站位,兩種模型的結果相比實驗值均存在巨大偏差.Agrawal和Sharma[30],Jiang等[5]學者采用LES方法預測的速度值在中線處也明顯低于實驗值.Agrawal對該站位的實驗值提出了質疑,因為尾跡速度虧損隨著遠離圓柱而減弱,即在x/c=?0.87處的速度最小值應該小于x/c= ?0.255處.在x/c= ?0.255站位處,兩種模型得到的速度型沒有明顯差異,并且速度最小值均比實驗結果低.x/c=0.25位于翼型最大厚度處,圓柱尾跡流過翼型前緣到最大厚度處是加速的過程.觀察該站位實驗值,在y/c=0.3附近速度達到最大值.兩種模型均預測出了該趨勢,與實驗值符合較好.x/c=2.0為翼型尾跡,兩種方法也均能較好地預測出速度型,SST-DDES比Tran-DDES稍接近實驗值,但差異并不大.

圖7 圓柱-翼型干擾流動Q等值面分布,顏色由馬赫數(shù)標識Fig.7.Iso-surface of the Q-criterion=0 for the rod-airfoil,colored by Mach number.

圖8 圓柱-翼型平均流場中四個典型站位速度分布圖Fig.8.Comparisons of mean velocity at four slice stations in the flow-field of rod-airfoil.

圖9 圓柱-翼型算例圓柱附近平均流場的渦強分布Fig.9.Distributions of mean vorticity intensity for the single cylinder.

圖10 圓柱-翼型流動中四個典型站位速度均方根值分布圖Fig.10.Comparisons of root mean square value of velocity at four slice stations in the flow-field of rod-airfoil.

上述四個站位的平均速度分布,Tran-DDES與SST-DDES的差異并不明顯.進一步觀察兩種方法得到的平均渦量分布(圖9),并與單圓柱構型情況的數(shù)值結果(圖4)進行對比.可以發(fā)現(xiàn),全湍流模型無法模擬層流-轉捩過程而導致剪切層失穩(wěn)推遲的問題減弱甚至消除了.文獻[22]指出翼型的加入也會反過來影響圓柱流場,包括圓柱表面氣動力及渦街形態(tài).由于翼型的影響,SST-DDES方法模擬的圓柱尾跡回流區(qū)被壓縮,并且恰好更接近真實情況.這也是許多文獻中[31,32]采用全湍流方法仍能在圓柱-翼型干擾流動中取得較好結果的原因.

然后,對比四個典型剖面的速度均方根值分布(圖10).在x/c=?0.87站位,計算結果仍與實驗值存在較大偏差.在x/c=?0.255站位,Tran-DDES和SST-DDES均與實驗符合很好,主要差異在y/c=0處速度均方根值是單個峰值還是兩個峰值.來流流過圓柱后,分離渦在圓柱軸線兩側交替產生,因此在均方根圖上會出現(xiàn)兩個峰值的現(xiàn)象.然而實驗也未捕捉到兩個峰值的現(xiàn)象,這與該站位遠離圓柱更接近翼型,并且受翼型流場干擾有關.在翼型最大厚度處的x/c=0.25站位,計算值均與實驗符合得很好,僅在近壁面處Tran-DDES預測的速度均方根值優(yōu)于SST-DDES.在翼型尾跡的x/c=2.0站位,實驗預測出了明顯的兩個峰值,并且由于圓柱和翼型存在高低差,均方根分布并不對稱.Tran-DDES成功預測出了兩個峰值與實驗趨勢符合.風洞中來流通過實驗模型湍流度會迅速升高,并且小尺度渦的速度脈動也會反映在均方根圖上.但對于湍流模型而言,湍流度無法直觀反映在流場中,并且小尺度渦會被平均(過濾)掉,因此在遠離軸線后脈動值會趨近于0.

最后對比聲學模擬結果,遠場噪聲采用FW-H方法模擬,并且同樣通過Kato公式,根據(jù)聲源信息的相干性對結果進行校正.圖11給出了采用不同模型所得到的Pfar點聲場頻譜形狀,并與實驗結果進行了對比.Tran-DDES和SST-DDES得到的主頻能較好地符合卡門渦街脫落的主頻,并在整體上的形狀也與實驗符合較好.但功率譜密度的量級存在一定的偏差,相比實驗值均偏大,部分原因可能是展向修正造成的.數(shù)值計算所得頻譜高頻振蕩較為嚴重,這是由于數(shù)值計算采樣時間不足采樣信號不穩(wěn)定造成的.后處理時需要對隨機信號進行多次平均才能得到較為穩(wěn)定的信號.這就需要延長采樣時間,并在進行快速傅里葉變換時取多次平均可以改善高頻振蕩.

圖12給出了以壁面為積分面在以翼型中點為圓心R=185D圓周上所得到的指向性結果與實驗的對比,圖中橫坐標OASPL(overall sound pressure level)表示總聲壓級.可以看出兩種計算模型預測的結果差異不大,并且均與實驗得到的指向性規(guī)律一致,即表現(xiàn)為偶極子指向性,且偶極子長軸垂直于來流方向.說明圓柱脫落渦對翼型表面的周期性撞擊所產生的非定常載荷是產生噪聲的主要原因.

圖11 圓柱-翼型算例Pfar點功率譜密度結果對比Fig.11.Comparison of PSD of acoustic pressure at Pfarfor the rod-airfoil.

圖12 圓柱-翼型算例聲壓級指向分布圖Fig.12.Directivity curves of OASPL in the flow-field of rod-airfoil.

5 結 論

本文基于七階WCNS-E8T7格式,結合DDES方法和FW-H聲比擬方法,對亞臨界雷諾數(shù)下單圓柱、圓柱-翼型的分離渦/渦致噪聲問題進行了數(shù)值模擬.針對亞臨界雷諾數(shù)下圓柱尾跡中的轉捩問題,發(fā)展了基于γ-Reθ模型的高精度Tran-DDES方法,并與傳統(tǒng)基于全湍流SST模型的高精度DDES方法進行了對比.

在單圓柱算例中,由于無法估計圓柱壁面附近的層流-轉捩過程,傳統(tǒng)DDES方法延遲了圓柱兩側剪切層的失穩(wěn),致使平均流場中回流區(qū)增大,壓差阻力偏小.同時,渦街產生的頻率變高,使聲壓頻譜整體向高頻偏移.而耦合了γ-Reθ轉捩模型的DDES方法能預測出與實驗符合的結果.

在圓柱-翼型算例中,由于翼型對圓柱附近流場的影響,傳統(tǒng)DDES方法對圓柱兩側剪切層失穩(wěn)的推遲減弱甚至消除了.因此在平均速度分布方面,SST-DDES方法與Tran-DDES方法并無明顯區(qū)別.但在脈動量預測以及脈動產生的噪聲預測方面,Tran-DDES方法與實驗符合得更好.

由于數(shù)值計算中無法承受加大展長所增加的巨大花費,本文采用了Kato公式擬合真實展長下的噪聲.但在聲壓主頻及功率譜密度量級方面,并不能完全與實驗符合,這也是下一步需要繼續(xù)改進的方面.

感謝中山大學國家超級計算廣州中心在計算資源方面提供的幫助.