各向異性湍流積分長度數值預報研究

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

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

各向異性湍流積分長度數值預報研究

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

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

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

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

0 引 言

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

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

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

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

1 各向異性湍流積分長度數值預報方法

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

數值預報方法流程如圖2所示。

圖2 湍流積分長度數值預報方法流程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的比率）控制方程中引入了計入阻尼影響的導數擴散項，通過引入湍流切變應力修正了渦粘性系數，并且修改了部分封閉方程的常系數。SST k-ω模式實際上混合了k-ω模式和k-ε模式，使得模式同時具有了k-ω模式計算近壁面區域粘性流動的可靠性和k-ε模式計算遠場自由流動的精確性。湍動能k和耗散率ε的控制方程為：

1.2 LES及DKET 亞格子模式

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

方程（4）右端有不封閉項

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

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

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

亞格子應力表達式為:

亞格子湍動能通過求解湍動能輸運方程來封閉：

式中：νt為粘性系數，

1.3 湍流積分長度計算式

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

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

對于平穩過程中系綜平均等于隨機過程的時間平均，則湍流的系綜平均可以表示為

則流場中任一點的脈動速度為

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

定義式中delta函數為

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

式中：r為兩點的間距。

2 算例分析

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

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

圖3 試驗模型示意圖Fig.3 Schematic diagram of experimental model

圖4二維水翼計算域網格Fig.4 Mesh 2-D hydrofoil

圖5比較了試驗和數值計算的結果。

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

表1 各向同性湍流積分長度計算試驗結果比較Tab.1 Comparison between numerical and experimental results of isotropy turbulence integral length

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

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

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

Blake（2002）在文獻[9]中給出了各向異性湍流積分長度的試驗及數值預報方法，計算公式在1.3中給出，試驗相關參數見表2，圖6給出了試驗結果、Blake計算結果和本文計算結果。

表2 試驗相關參數Tab.2 Test parameters

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

圖6 計算及試驗結果Fig.6 Numerical prediction and experimental data of

圖7 SUBOFFFig.7 SUBOFF

2.3 SUBOFF艉部湍流積分長度預報

2.3.1 模型尺寸及網格劃分

本節采用美國國防部先進研究項目局（DARPA）提出的SUBOFF AFF-8標準模型為研究對象，研究潛艇艉部槳盤面處湍流積分長度的數值預報，并且對相關特性進行了研究。

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

2.3.2 流場計算結果驗證

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

圖8 指揮臺圍殼不同高度壓力系數比較：（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處的計算及試驗值比較Fig.9 Comparison between calculation and experiment at r/R=0.25

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

2.3.3 槳盤面處湍流積分長度

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

圖10 計算結果Fig.10 Calculation results of

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

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

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

3 結 語

本文對潛艇艉部湍流積分長度進行了數值預報及研究，首先通過與二維標準水翼和風洞三維螺旋槳試驗數據比較，分別驗證了二維水中各向同性和三維風洞各向異性湍流積分長度的數值預報方法；然后對SUBOFF艉部各向異性湍流積分長度進行了預報研究，得出以下結論：

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

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

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

[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].水動力學研究與進展,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]朱錫清,李 亞,孫紅星.船舶螺旋槳葉片與艉部湍流場相互作用噪聲的預報研究[J].聲學技術,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]熊紫英,孫紅星,朱錫清.入射湍流與螺旋槳相互作用的低頻寬帶噪聲預報研究[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]張國平.空泡水筒中湍流參數PIV試驗測量和分析方法研究報告[R].無錫:中國船舶科學研究中心,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聯合測試分析方法研究[C]//第十三屆船舶水下噪聲學術討論會論文集.中國鷹潭,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]劉建華,謝 華,田于逵.湍流積分長度新型聯合測試方法研究[C]//第二十三屆全國水動力學研討會暨第十屆全國水動力學學術會議論文集.中國西安,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].北京:清華大學出版社,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

國家自然科學基金資助項目（No.51079158）

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

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