蔣靖偉，馬 騁，錢正芳，2，陳 科，蔡昊鵬，張 赫

（1海軍裝備研究院，北京 100161；2海軍工程大學(xué)，武漢 430033；3中國艦船研究院，北京100101）

各向異性湍流積分長度數(shù)值預(yù)報研究

蔣靖偉1，馬 騁1，錢正芳1，2，陳 科1，蔡昊鵬1，張 赫3

（1海軍裝備研究院，北京 100161；2海軍工程大學(xué)，武漢 430033；3中國艦船研究院，北京100101）

文章針對目前模型試驗(yàn)只能獲得各向同性湍流積分長度這一現(xiàn)狀，建立了各向異性湍流積分長度數(shù)值預(yù)報方法：（1）RANS預(yù)報初始定常流場；（2）大渦模擬計(jì)算脈動速度場；（3）脈動速度相關(guān)函數(shù)計(jì)算各向同性湍流積分長度；（4）計(jì)算各向異性湍流積分長度。通過二維水翼模型試驗(yàn)驗(yàn)證了水中各向同性湍流積分長度的數(shù)值預(yù)報精度，平均誤差6.81%；通過風(fēng)洞螺旋槳驗(yàn)證了空氣中各向異性湍流積分長度計(jì)算方法。對SUBOFF槳盤面處各向異性湍流積分長度的預(yù)報研究發(fā)現(xiàn)，各方向湍流積分長度外半徑大于內(nèi)半徑；潛艇上部渦團(tuán)擴(kuò)散較下部明顯，內(nèi)半徑湍流積分長度波峰隨著角度增大而內(nèi)移周向分布受指揮臺圍殼和穩(wěn)定翼馬蹄渦影響，波峰出現(xiàn)在20°和90°左右的頻率峰值受指揮臺圍殼和穩(wěn)定翼馬蹄渦脫落頻率（約15 Hz）影響，峰值集中在30 Hz以下。該文提出的各向異性湍流積分長度預(yù)報方法，為下一步螺旋槳低頻寬帶噪聲預(yù)報中構(gòu)造新的湍流波數(shù)譜以提高預(yù)報精度提出了思路。

各向異性湍流積分長度；數(shù)值預(yù)報；CFD；大渦模擬

0 引 言

螺旋槳低頻寬帶噪聲的預(yù)報是目前潛艇螺旋槳噪聲研究的重點(diǎn)和熱點(diǎn)[1-8]，是潛艇進(jìn)一步減震降噪需要突破的關(guān)鍵技術(shù)。譜方法是目前預(yù)報低頻寬帶噪聲的主流方法。譜方法由Blake（1986）[1]提出完整的理論，隨后Kirschner（1993）[2]提出條帶理論，建立了螺旋槳低頻寬帶噪聲的數(shù)值預(yù)報方法，其計(jì)算流程如圖1。

圖1 條帶理論計(jì)算流程Fig.1 Procedure diagram of strip theory

其中作為表征湍流信息的湍流波數(shù)譜對螺旋槳低頻寬帶噪聲的精確預(yù)報起著重要作用，湍流波數(shù)譜由湍流積分長度和湍流度構(gòu)成，其表達(dá)式為（以三維Liepmann譜[3]為例）

目前湍流積分長度的獲取主要是通過模型試驗(yàn)[6-8]，模型試驗(yàn)主要有以下不足：（1）無論是風(fēng)洞熱線單點(diǎn)測量法還是空泡水筒中PIV試驗(yàn)方法，都是基于各向同性假設(shè)，測量的是各向同性湍流積分長度；（2）試驗(yàn)費(fèi)時費(fèi)力，加工模型造成額外成本；（3）不同測試系統(tǒng)（風(fēng)洞、水池、循環(huán)水槽等）的試驗(yàn)測量結(jié)果對于潛艇艉部各向異性湍流場的適用性并沒有進(jìn)行過系統(tǒng)的研究分析。為彌補(bǔ)試驗(yàn)及目前理論研究中的不足，本文建立了各向異性湍流積分長度的數(shù)值預(yù)報方法，通過數(shù)值計(jì)算結(jié)果對比二維水翼驗(yàn)證了水槽各向同性湍流積分長度，對比空氣中螺旋槳試驗(yàn)數(shù)據(jù)驗(yàn)證了風(fēng)洞各向異性湍流積分長度數(shù)值預(yù)報方法；其次，用建立的數(shù)值方法預(yù)報了SUBOFF潛艇槳盤面處各向異性湍流積分長度，對不同半徑、不同頻率的各向異性湍流積分長度進(jìn)行了研究分析；文章的最后對全文進(jìn)行了總結(jié)展望。

1 各向異性湍流積分長度數(shù)值預(yù)報方法

各向湍流積分長度的數(shù)值預(yù)報，分為四個步驟：定常流場預(yù)報、脈動流場預(yù)報和積分長度計(jì)算。（1）在定常流場求解中，本文使用Shear-Stress Transport（SST）k-ω湍流模式來封閉雷諾應(yīng)力；（2）在脈動流場求解中，考慮到直接數(shù)值模擬（Direct Numerical Simulation,DNS）的計(jì)算量是目前難以承受的，本文采用大渦模擬（Large Eddy Simulation，LES）方法求解脈動速度場，亞格子模式采用動態(tài)湍動能輸運(yùn)（Dynamic Kinetic Energy Transport，DKET）；（3）在各向同性積分長度的計(jì)算中，通過脈動速度相關(guān)函數(shù)來計(jì)算；（4）在各向異性積分長度的計(jì)算中，采用Blake（2002）[9]提出的各向異性湍流積分長度計(jì)算公式。

數(shù)值預(yù)報方法流程如圖2所示。

圖2 湍流積分長度數(shù)值預(yù)報方法流程Fig.2 Procedure diagram of numerical prediction on anisotropic turbulence integral length

1.1 k-ω模式

Wilcox（2007）[10]提出了k-ω模式的修正模式Shear-Stress Transport（SST）k-ω模式，在ω（ε和k的比率）控制方程中引入了計(jì)入阻尼影響的導(dǎo)數(shù)擴(kuò)散項(xiàng)，通過引入湍流切變應(yīng)力修正了渦粘性系數(shù)，并且修改了部分封閉方程的常系數(shù)。SST k-ω模式實(shí)際上混合了k-ω模式和k-ε模式，使得模式同時具有了k-ω模式計(jì)算近壁面區(qū)域粘性流動的可靠性和k-ε模式計(jì)算遠(yuǎn)場自由流動的精確性。湍動能k和耗散率ε的控制方程為：

1.2 LES及DKET 亞格子模式

在獲得定常流場后，通過LES計(jì)算脈動的速度場，LES的基本思想是直接計(jì)算大尺度脈動，而對小尺度脈動做模式，因此首先要把小尺度脈動過濾掉，本文假定過濾過程和求導(dǎo)過程可以交換，對Navier-Stokes方程作過濾運(yùn)算，可得到LES控制方程為[11]

方程（4）右端有不封閉項(xiàng)

式中：τij為亞格子應(yīng)力，代表過濾掉的小尺度脈動和可解尺度湍流間的動量輸運(yùn)，需要?；?，本文采用DKET亞格子模式[12]進(jìn)行封閉。

該模式通過亞格子湍動能來定義粘性系數(shù)

式中：Δ=V1/3是亞格子尺度是亞格子湍動能。

亞格子應(yīng)力表達(dá)式為:

亞格子湍動能通過求解湍動能輸運(yùn)方程來封閉：

式中：νt為粘性系數(shù)，

1.3 湍流積分長度計(jì)算式

湍流積分長度是指流場中最大尺度渦結(jié)構(gòu)的量度，即大渦的尺度，湍流積分長度主要有下面三種定義[9]：（1）基于流場中兩點(diǎn)速度的相關(guān)函數(shù)；（2）基于湍流中的某些物理特征長度（例如流場中渦的半寬長度或者固壁面邊界層中的動量厚度）；（3）基于湍流波譜最佳擬合函數(shù)，通過調(diào)整湍流積分長度和湍流度來使得假設(shè)的曲線能夠以最小誤差吻合湍流波譜。上述三種積分長度的定義都是針對基于各向同性假設(shè)和Taylor冷凍假設(shè)的傳統(tǒng)積分長度，對于普遍存在的各向異性湍流積分長度并不能很好地進(jìn)行函數(shù)表達(dá)。因此，本文采用Blake提出的基于互功率譜卷積定義的各向異性湍流積分長度[9]，具體定義和計(jì)算公式如下。

（10）式定義了湍流的系綜平均：

對于平穩(wěn)過程中系綜平均等于隨機(jī)過程的時間平均，則湍流的系綜平均可以表示為

則流場中任一點(diǎn)的脈動速度為

定義“速度的j分量在ri方向上的湍流積分長度”，

定義式中delta函數(shù)為

式中：αc是r2方向上的湍流渦團(tuán)拉伸參數(shù)，若要計(jì)算r2方向上的積分長度，只需將拉伸參數(shù)αc改成βc即可是各向同性湍流積分長度，可以由下式求得：

式中：r為兩點(diǎn)的間距。

2 算例分析

2.1 二維水翼各向同性湍流積分長度

本節(jié)通過與國際上已被廣泛應(yīng)用于隨邊湍流和輻射噪聲研究的二維水翼的PIV和TR-PIV（高速PIV）試驗(yàn)數(shù)據(jù)進(jìn)行對比分析，驗(yàn)證了各向同性湍流積分長度的數(shù)值預(yù)報方法。二維水翼弦長為364 mm，厚度為20.32 mm，展長為450 mm。水翼頭部為5:1的橢圓，中間為平板，下表面為平面，上表面尾緣部分為與下表面夾角為45°的圓弧，水翼模型的剖面示意圖及坐標(biāo)如圖3（a），圖3（b）是PIV測量流場區(qū)域示意圖，圖3（c）給出了湍流積分長度計(jì)算取點(diǎn)位置，試驗(yàn)在中國船舶科學(xué)研究中心空泡水筒中進(jìn)行[6]。計(jì)算域網(wǎng)格如圖4，網(wǎng)格總量27萬。

圖3 試驗(yàn)?zāi)Ｐ褪疽鈭DFig.3 Schematic diagram of experimental model

圖4二維水翼計(jì)算域網(wǎng)格Fig.4 Mesh 2-D hydrofoil

圖5比較了試驗(yàn)和數(shù)值計(jì)算的結(jié)果。

圖5 不同位置剖面湍流積分長度預(yù)報結(jié)果與試驗(yàn)比較Fig.5 Comparison between numerical and experimental data of turbulence integral length at different locations

表1 各向同性湍流積分長度計(jì)算試驗(yàn)結(jié)果比較Tab.1 Comparison between numerical and experimental results of isotropy turbulence integral length

從圖5中可以看出，DKET亞格子模式基本反應(yīng)出湍流積分長度的結(jié)構(gòu)特征：（1）對于近隨邊（x= 6.3,15 mm）區(qū)域的高湍流度部分（-20 mm≤y≤15 mm），CFD計(jì)算出了顯著的雙峰結(jié)構(gòu)，并且在x=15 mm處雙峰間距大于x=6.3 mm處，反應(yīng)出漩渦的擴(kuò)散；但是計(jì)算的兩峰間隔較試驗(yàn)結(jié)果偏大，即CFD流場計(jì)算時漩渦渦核間距偏大。（2）對于稍遠(yuǎn)離隨邊（x=25.7,36.2 mm）區(qū)域，在（y≤0 mm）部分，計(jì)算結(jié)果與試驗(yàn)結(jié)果較吻合，波峰稍偏左，這是由于之前的渦核間距計(jì)算偏大；在（y＞0 mm）部分，計(jì)算峰值比試驗(yàn)結(jié)果小1 mm左右。（3）對于y=0.8,-9.7 mm兩個剖面，y=0.8計(jì)算偏差較大，而y=-9.7 mm計(jì)算較為精確。總體情況來看，在y≤0 mm區(qū)域計(jì)算的較為準(zhǔn)確，而在y＞0 mm區(qū)域偏差較大。表1給出了不同剖面上試驗(yàn)和計(jì)算結(jié)果的均值。

從圖5及表1可以看出，數(shù)值計(jì)算的湍流積分長度能滿足預(yù)報精度要求。

2.2 空氣螺旋槳各向異性湍流積分長度

Blake（2002）在文獻(xiàn)[9]中給出了各向異性湍流積分長度的試驗(yàn)及數(shù)值預(yù)報方法，計(jì)算公式在1.3中給出，試驗(yàn)相關(guān)參數(shù)見表2，圖6給出了試驗(yàn)結(jié)果、Blake計(jì)算結(jié)果和本文計(jì)算結(jié)果。

表2 試驗(yàn)相關(guān)參數(shù)Tab.2 Test parameters

圖中實(shí)線是本文計(jì)算結(jié)果，虛線是Blake[9]的計(jì)算結(jié)果，從圖中可以看出，本文計(jì)算的無論是與Blake的計(jì)算結(jié)果還是試驗(yàn)結(jié)果都較為接近，驗(yàn)證了各向異性積分長度計(jì)算方法的計(jì)算精度。

圖6 計(jì)算及試驗(yàn)結(jié)果Fig.6 Numerical prediction and experimental data of

圖7 SUBOFFFig.7 SUBOFF

2.3 SUBOFF艉部湍流積分長度預(yù)報

2.3.1 模型尺寸及網(wǎng)格劃分

本節(jié)采用美國國防部先進(jìn)研究項(xiàng)目局（DARPA）提出的SUBOFF AFF-8標(biāo)準(zhǔn)模型為研究對象，研究潛艇艉部槳盤面處湍流積分長度的數(shù)值預(yù)報，并且對相關(guān)特性進(jìn)行了研究。

SUBOFF AFF-8帶有指揮臺圍殼和十字型穩(wěn)定翼，總長L=4.356 m，進(jìn)流段長1.016 m，平行舯體長2.229 m，最大直徑D=0.508 m，去流段長1.111 m（后體端部長 0.095 m），穩(wěn)定翼后緣位于4.007 m處。指揮臺圍殼前緣位于0.924 m處，后緣位于1.923 m處，長0.368 m，高0.460 m，其截面為橢圓，長短軸比為2:1，頂部為有外凸的橢圓蓋。槳盤面位于x/L=0.978處，直徑d=D/2。具體外形參見文獻(xiàn)[13]。坐標(biāo)系定義及網(wǎng)格如圖7，網(wǎng)格總量1 200萬。

2.3.2 流場計(jì)算結(jié)果驗(yàn)證

圖8對比了指揮臺圍殼0.1H、0.5H和0.9H高度處（指揮臺圍殼總高度為H）壓力系數(shù)的試驗(yàn)和計(jì)算結(jié)果，圖9對比了槳盤面上r/R=0.25處三向平均速度和湍動能的周向分布，在圖10中對比了軸向速度周向平均值沿徑向分布。

圖8 指揮臺圍殼不同高度壓力系數(shù)比較：（a）h=0.1H，（b）h=0.5H，（c）h=0.9HFig.8 Pressure coefficient comparison at different height of the fairwater

圖9槳盤面上r/R=0.25處的計(jì)算及試驗(yàn)值比較Fig.9 Comparison between calculation and experiment at r/R=0.25

圖8和圖9的計(jì)算結(jié)果與試驗(yàn)及與Bull[13]的比較分析可以看出，本文的CFD計(jì)算結(jié)果無論是在艇前部指揮臺圍殼處，還是艇艉槳盤面處的流場計(jì)算，均有較高的精度，可以作為計(jì)算湍流積分長度的流場輸入。

2.3.3 槳盤面處湍流積分長度

本節(jié)數(shù)值預(yù)報了SUBOFF槳盤面處各向異性湍流積分長度，圖10給出基于速度相關(guān)函數(shù)的湍流積分長度不同半徑處周向分布和周向平均值的徑向分布。

圖10 計(jì)算結(jié)果Fig.10 Calculation results of

從圖10（a），（b），（c）可以看出，受到指揮臺圍殼和穩(wěn)定翼馬蹄渦的影響，湍流積分長度在20°和90°附近出現(xiàn)峰值；受渦團(tuán)擴(kuò)散影響，外半徑湍流積分長度要大于內(nèi)半徑。圖10（d）可以看出三向速度的湍流積分長度的周向平均隨徑向變化趨勢大致一致；在內(nèi)半徑處三個湍流積分長度數(shù)值上大小較為一致，在外半徑數(shù)值上差別較小，而與差別較大，這是由于內(nèi)半徑湍流的各向異性較弱，而外半徑隨著渦團(tuán)的擴(kuò)散，湍流的各向異性逐漸增強(qiáng)。

圖11 不同半徑處的頻率分布Fig.11 Frequency distribution ofat different radius

圖12 不同頻率在不同周向位置的徑向分布Fig.12 Radial distribution ofat topic frequency at different circumferential locations

3 結(jié) 語

本文對潛艇艉部湍流積分長度進(jìn)行了數(shù)值預(yù)報及研究，首先通過與二維標(biāo)準(zhǔn)水翼和風(fēng)洞三維螺旋槳試驗(yàn)數(shù)據(jù)比較，分別驗(yàn)證了二維水中各向同性和三維風(fēng)洞各向異性湍流積分長度的數(shù)值預(yù)報方法；然后對SUBOFF艉部各向異性湍流積分長度進(jìn)行了預(yù)報研究，得出以下結(jié)論：

（1）指揮臺圍殼和穩(wěn)定翼馬蹄渦對艉部渦團(tuán)位置有明顯影響，周向波峰出現(xiàn)在20°和90°左右；受渦團(tuán)擴(kuò)散影響，徑向分布外半徑積分長度要大于內(nèi)半徑；

（2）指揮臺圍殼和穩(wěn)定翼馬蹄渦脫落頻率（約15 Hz）對的頻率峰值有明顯影響，峰值集中在30 Hz以下；徑向分布整體呈現(xiàn)外半徑大于內(nèi)半徑趨勢；潛艇上部渦團(tuán)擴(kuò)散較下部明顯，導(dǎo)致內(nèi)半徑波峰隨著角度增大而內(nèi)移。

下一步工作將集中在對各向異性湍流積分長度的機(jī)理研究和通過各向異性湍流積分長度構(gòu)造新的湍流波數(shù)譜并評估對低頻寬帶噪聲的影響這兩個方面。

[1]Blake W K.Mechanics of flow-induced sound and vibration[M].New York:Academic Press,1986:735-756.

[2]Kirschner I N,Corriceau P J,Muench J D.Validation of propeller turbulence ingestion acoustic radiation model using wind tunnel measurements[C].Louisiana:Flow Noise Modeling,Measurement and Control,ASME,1993:175-186.

[3]朱錫清,唐登海,孫紅星.船舶螺旋槳低頻噪聲研究[J].水動力學(xué)研究與進(jìn)展,2000,15(1):74-81.Zhu Xiqing,Tang Denghai,Sun Hongxing.Study of low-frequency noise induced by marine propeller[J].Chinese Journal of Hydrodynamics,2000,15(1):74-81.

[4]朱錫清,李 亞,孫紅星.船舶螺旋槳葉片與艉部湍流場相互作用噪聲的預(yù)報研究[J].聲學(xué)技術(shù),2006,25(4):361-364. Zhu Xiqing,Li Ya,Sun Hongxing.Prediction of noise induced by interaction between turbulence flow and propeller blades [J].Technical Acoustics,2006,25(4):361-364.

[5]熊紫英,孫紅星,朱錫清.入射湍流與螺旋槳相互作用的低頻寬帶噪聲預(yù)報研究[J].中國造船,2014,55(3):3-11. Xiong Ziying,Sun Hongxing,Zhu Xiqing.Prediction of low-frequency broadband noise induced by the interaction between injected turbulence and propeller[J].Shipbuilding of China,2014,55(3):3-11.

[6]張國平.空泡水筒中湍流參數(shù)PIV試驗(yàn)測量和分析方法研究報告[R].無錫:中國船舶科學(xué)研究中心,2010. Zhang Guoping.Comparison on turbulence measurement and analysis methods in cavitation tunnel by Particle Image Velocimeter[R].Wuxi:China Ship Scientific Research Center,2010.

[7]謝 華,劉建華,田于逵.湍流積分長度的HWA與LDA聯(lián)合測試分析方法研究[C]//第十三屆船舶水下噪聲學(xué)術(shù)討論會論文集.中國鷹潭,2011:345-350. Xie Hua,Liu Jianhua,Tian Yekui.Research and analysis on HWA and LDA joint test of turbulence integral length[C]// Proceedings of 13st Symposium on Underwater Noise.Yingtan,China,2011:345-350.

[8]劉建華,謝 華,田于逵.湍流積分長度新型聯(lián)合測試方法研究[C]//第二十三屆全國水動力學(xué)研討會暨第十屆全國水動力學(xué)學(xué)術(shù)會議論文集.中國西安,2011:283-289. Liu Jianhua,Xie Hua,Tian Yukui.Validation of a new method for turbulence integral length measurement based on a HWA and LDA combined system[C]//Proceedings of 23st Symposium on National Hydrodynamics.Xi’an,China,2011: 283-289.

[9]Lynch D A,Mueller T J,Blake W K.A correlation length scale for the prediction of aeroacoustic response[C]//AIAA. Colorado,2002.

[10]Wilcox D C.Formulation of the k-ω turbulence model revisited[C].Reno:45th AIAA Aerospace Sciences Meeting and Exhibit,2007.

[11]張兆順,崔桂香,許春曉.湍流理論與模擬[M].北京:清華大學(xué)出版社,2005.

[12]Won-Wook,Suresh Menon.Application of the localized dynamic subgrid-scale model to turbulent wall-bounded flows [R].Reno:Technical Report AIAA-97-0210,1997.

[13]Bull P.The validation of CFD prediction of nominal wake for the SUBOFF fully appended geometry[C]//Proceedings of 21st Symposium on Naval Hydrodynamics.Norway,1996.

Numerical prediction and research on anisotropic turbulence integral length

JIANG Jing-wei1,MA Cheng1,QIAN Zheng-fang1,2,CHEN Ke1,CAI Hao-peng1,ZHANG He3
(1 Naval Academy of Armament,Beijing 100161,China;2 Naval Univ.of Engineering,Wuhan 430033,China; 3 China Ship Research and Development Academy,Beijing 100101,China)

This paper proposes a numerical method to predict the anisotropic turbulence integral length in view of this situation that the classic(isotropic)turbulence integral length is mainly acquired from experiments.Where RANS is utilized to predict the initial steady flow field;LES is used to calculate fluctuating velocity;the classic turbulence integral length is calculated based on cross-correlation function and the anisotropic turbulence integral length is finally predicted.The numerical method to predict isotropic turbulence integral length is validated the experimental data from 2-D hydrofoil,and the average error is 6.81%; and anisotropic turbulence integral length prediction method is validated by wind tunnel test.Then the anisotropic turbulence integral length of SUBOFF is studied and understood:the integral lengths at outside radius are lager then those at inside radius in all directions;the vortexes in the upper part diffuse more intensely than that in the under part,leading to the peak bars moving inside along with the angle increasing; the peak bars oare all near 20°and 90°according to the exists of horse-shoe vortexes of conning tower and stabilizers;The peak frequency ofare both lower than 30 Hz as theshedding frequency of the horse-shoe vortexes is about 15 Hz.The current work lays the foundation of constructing new turbulence spectrum to improve the accuracy of prediction of low-frequency broadband noise.

anisotropic turbulence integral length;numerical prediction;CFD; Large Eddy Simulation,LES

U661.1

A

10.3969/j.issn.1007-7294.2015.10.001

1007-7294（2015）10-1161-12

2015-06-08

國家自然科學(xué)基金資助項(xiàng)目（No.51079158）

蔣靖偉（1991-），男，碩士研究生，E-mail:15365286622@163.com；

馬 騁（1963-），男，研究員，博士生導(dǎo)師，E-mail:Macheng50089@yahoo.com.cn。