潛艇阻力及流場數(shù)值仿真策略優(yōu)化分析

李 鵬， 王 超， 孫華偉， 郭春雨

(哈爾濱工程大學 船舶工程學院， 哈爾濱 150001)

隨著計算機的發(fā)展，計算流體力學(Computational Fluid Dynamics, CFD)被越來越多地應用到海洋結(jié)構(gòu)物的性能研究中.潛艇是重要的海洋結(jié)構(gòu)物，準確了解其不同航速下的水動力性能、精確掌握其流場特征，對提高潛艇自身性能具有重要意義[1].

國內(nèi)外學者針對潛艇性能開展過大量的數(shù)值研究.Pan等[2]采用RANSk-ω(k為湍動能，ω為單位湍動能耗散率)和雷諾應力湍流模型(Reynolds Stress Model, RSM)對潛艇的阻力進行預報，取得了較好的結(jié)果.Posa等[3-5]基于大渦模擬湍流模型(Large Eddy Simulation, LES)對潛艇的尾流場進行仿真計算，結(jié)果表明LES能很好地捕捉流動特征，同時為流噪聲的預報提供精確輸入.Fu等[6-8]就分離渦湍流模型(Detach Eddy Simulation, DES)下的潛艇水動力性能及流場特征進行仿真，結(jié)果證實前者對垂直壁面方向的無量綱化距離y+的敏感性微弱，后者對空間離散精度的高需求性.此外，李士強等[8]還就k-ω、k-ε(ε為湍流耗散率)模型對潛艇附體變動的影響進行仿真，并證實關(guān)于舵失速的問題，在失速前k-ε獲取的渦流脫落頻率較k-ω更高.

綜上，國內(nèi)外學者就直航潛艇的性能進行較多分析，但未闡明不同湍流模型在潛艇水動力性能及流場特征計算中的優(yōu)勢或短板.本文針對這一問題，首先展多套網(wǎng)格及時間步長的收斂性分析.接著，針對不同的湍流模型開展特定工況的數(shù)值仿真計算，比較預報結(jié)果的差異.最后，基于LES-Smagorinsky模型的預報結(jié)果對潛艇渦系特征進行分析.

1 理論基礎(chǔ)

1.1 控制方程

時均化后不可壓縮流體的連續(xù)性方程與動量方程為

(1)

(2)

1.2 湍流模型

由于在時間和空間上對計算成本的需求較高，目前主要的CFD求解器均采用求解平均量或濾波量的方法來模擬小尺度流動以節(jié)省計算成本.本文就多個湍流模型開展?jié)撏У乃畡恿π阅芗傲鲌鎏卣饔嬎悖瑢Ρ雀鞣N湍流模型的計算結(jié)果，進而篩選適合潛艇水動力性能及流場特征仿真計算的湍流模型.

雷諾平均納維爾-斯托克斯(Reynolds-Averaged Navier-Stokes, RANS)是一類湍流模型集合的總稱，涉及多個子模型，這些模型的特征都是在N-S(Navier-Stokes)方程的基礎(chǔ)上添加了一個額外的黏度項且都包含一個通用項湍動能k，即單位質(zhì)量湍流脈動的動能，而這些模型的區(qū)別在于對額外的黏度項處理方式不同：SA(Spalart-Allmaras)模型是一方程湍流模型，無壁面函數(shù)要求但增加一個新的變量SA黏度，其優(yōu)點是對計算內(nèi)存要求較低，有良好的收斂性，缺點是不適合求解剪切流與分離流；k-ε模型基于k和ε建立，是具有良好收斂性的兩方程模型，自提出后逐漸演化出realizablek-ε，RNGk-ε等，局限性體現(xiàn)在對無滑移壁面、逆壓梯度等的流動模擬不精確，耗散不易計算；k-ω模型計算結(jié)果對初始條件較為敏感且收斂相對困難，但相對于ε，ω更容易求解.LES是基于自相似理論，通過過濾函數(shù)將大渦和小渦分離開，大渦直接解析，小渦是各項同性，因此可對其進行?；嬎悖诮馕龃鬁u和模化小渦的過程中，產(chǎn)生的亞格子尺度應力項通過亞格子(Subgrid Scale, SGS) 模型計算.DES模型結(jié)合RANS/LES模型各自的優(yōu)點，在近壁面采用RANS進行計算而利用LES模型的算法對湍流核心區(qū)域計算，經(jīng)過多年發(fā)展逐漸演變出延遲分離渦湍流 (Delayed Detach Eddy Simulation, DDES)、改進的延遲分離渦湍流 (Improved Delayed Detached Eddy Simulation, IDDES) 等模型變體，對非定常湍流流動的模擬有較大的優(yōu)勢.

本文基于植入到STAR-CCM+中的湍流模型：RANSk-ω、k-ε、SA、LES-Samgorinsky(為行文簡潔，后續(xù)以LES-S替代)、WALE(LES-W)、Dynamic Smagorinsky(LES-DS)以及DES-SA、DES-EBKE、DDES、IDDES模型開展Sub-off水動力性能及流場特征的計算，進而甄選合適的湍流模型開展?jié)撏畡恿π阅芗傲鲌鎏卣鞯挠嬎?

2 仿真基礎(chǔ)

2.1 坐標系和計算域

本文采用的幾何模型是Sub-off模型，詳細的模型資料可參考文獻[9-11]，后續(xù)的研究中，在原始幾何模型尺度的基礎(chǔ)上進行縮放，縮放比λ=0.688 7，模型總長L=3 m，最大直徑D=0.35 m，若無特殊說明本文的計算結(jié)果均是基于這一參數(shù)下的幾何模型.

本文設定笛卡爾坐標系Oxyz，其原點位于艇首，如圖1(a)所示，圖中Uref為自由流場的自由來流速度.本文設定計算域為圓柱，具體尺寸及邊界條件的設置見圖1(b)，參考文獻[12]的研究，當目標幾何模型在計算域流向投影面積小于1%時，數(shù)值計算結(jié)果可忽略阻塞效應的影響，通過簡單的數(shù)據(jù)計算可知，Sub-off模型在計算域上下底面的投影面積占比為0.7%(小于1%)，故可忽略阻塞效應的影響.

圖1 幾何模型及計算域

2.2 網(wǎng)格劃分和工況設定

圖2 網(wǎng)格方案3簡圖

為統(tǒng)一和方便作者進行數(shù)值和試驗研究，結(jié)合計算和實際條件需求，對航速為2.53 m/s的潛艇模型的水動力性能及流場特征進行收斂性計算，以對比不同湍流模型獲得的仿真結(jié)果.

3 收斂性分析

表2 4套網(wǎng)格方案的計算結(jié)果

網(wǎng)格收斂率K=εi, i+1/εi, i-1,εi, i+1是較密兩方案結(jié)果差，εi, i-1是較稀兩方案結(jié)果差.K的計算結(jié)果可能是K>1，0

表3 網(wǎng)格收斂性計算結(jié)果

P=ln (εi, i+1/εi, i-1)lnδ

評估3套網(wǎng)格的收斂指數(shù)：

Ii, i+1=DS|ei, i+1|/(δP-1)

式中：DS=1.25為網(wǎng)格收斂安全系數(shù)[18]；ei, i+1=εi, i+1.理論上，當P=2時，網(wǎng)格正交性或湍流模型非線性對計算結(jié)果準確性的影響程度最小[14,18]，選用中等密度網(wǎng)格進行數(shù)值計算具有較高可靠性.可知方案2～4具有更好的收斂率且離散精度更接近理論值，故后續(xù)的計算中采用網(wǎng)格方案3進行計算.

表4是不同時間步長計算得到的判定系數(shù)結(jié)果，以0.001 s的計算結(jié)果為基準進行對比.可以看出, 不同時間步長對計算結(jié)果的影響較小，為提高計算效率，后續(xù)的計算中均采用0.001 s為時間步長.

表4 不同時間步長計算結(jié)果

4 結(jié)果與討論

4.1 水動力性能

不同湍流模型對近壁面流動、黏性應力等的解析不同，造成仿真得到的阻力有所差異，表5為不同湍流模型計算得到的潛艇阻力.

表5 不同湍流模型的艇體阻力仿真結(jié)果

表6 不同湍流模型計算結(jié)果與實驗結(jié)果對比

4.2 流場特征

采用Q標準來判定流體為渦旋或應變主導[19]，圖3所示為不同湍流模型仿真計算得到的Sub-off周圍的渦場.圖3清晰地表明LES-W、LES-DS對潛艇重要渦特征——馬蹄渦的捕捉不足，不同的是除此外的8個湍流模型均很好地捕捉到了圍殼附近馬蹄渦的流動特征，另外這兩種湍流模型預報的潛艇尾肩部的渦結(jié)構(gòu)亦沒有其他模型強烈.選用的10種湍流模型中，RANS類湍流模型對潛艇尾部渦強度的預報結(jié)果顯然低于其他模型，這可能是RANS速度平均后的結(jié)果，是RANS模型的算法特點所決定.值得注意的是，10個模型均沒能很好地捕捉或預報圍殼稍渦的演化，但對于尾舵稍渦的演化均能較好地捕捉.推測原因：一是圍殼稍渦的強度沒有尾舵稍渦強度大，二是空間分辨率或時間分辨率不足以捕捉其演化過程.限于計算資源以及本文研究重點，在此不做深入探討.

圖3 不同湍流模型計算艇體周圍Q采樣體積(Q=5 s-2)

伴流是考核潛艇綜合性能的重要指標，對伴流的準確預報十分重要.圖4所示為不同湍流模型計算得到的Sub-off尾流伴流分數(shù)δ隨距離變化的曲線，圖中，0為潛艇中軸線的最尾端，D、2D、3D、4D為沿著x正方向與潛艇中軸線最尾端的距離.可見，LES-DS、LES-W的預報結(jié)果與其他模型有明顯偏差，但同屬大渦模擬類中LES-S的預報結(jié)果有良好的“隨群”性，結(jié)合對于宏觀力的預報結(jié)果可知，LES-S模型在預報潛艇宏觀變量方面有著較好的適用性.

圖4 不同湍流模型預報尾流伴流(0位于艇體尾端)

為了深入探究不同湍流模型對Sub-off流場的預報結(jié)果， 本文提取了如圖5虛線所示位置的流場時均變量并加以分析，時均流場的采樣總長是流經(jīng)一倍艇長時間. 圖6所示為時均軸向速度的周向變化曲線，圖中θ為圖5中提取截面沿著白色虛線環(huán)向位置的角度，u為流場局部軸向速度，通過u/Uref對軸向速度進行無因次化對比, 分析其變化規(guī)律，為了有效說明流場，將Sub-off周圍流場分為內(nèi)流場(r/D<0.501)和外流場(r/D>0.501)，其中r為對應曲線的半徑值.可知，10種湍流模型對Sub-off的外流場及其平行中體內(nèi)流場的預報結(jié)果吻合程度非常高，推測由于遠離壁面或壁面流動簡單，各湍流模型對其處理無明顯差異，但即使對平行中體內(nèi)流場的流動預報差異微弱，Sub-off尾內(nèi)流場的預報差異卻十分明顯.這表明，湍流模型對因為幾何特征變化引發(fā)流動變化更為敏感，對流場“傳遞過程”引發(fā)的變化敏感性較弱.圖7所示為流場時均壓力的周向分布曲線，可以看出，Sub-off周圍流場的時均壓力系數(shù)Cp(Cp=p/(0.5ρUref))分布曲線對湍流模型地變化不敏感，各個模型均有良好的“隨群”性，由于尾部附體以及尾部幾何收縮造成的負壓力梯度，在x/L=0.9處，流場壓力波動仍較為劇烈.圖6和7的結(jié)果表明：時均速度較時均壓力對湍流模型的改變更敏感，尤其是近壁面流動，如圖6(a)、7(a)所示，這是不同湍流模型調(diào)用不同壁面解析函數(shù)造成的.圖6(b)、7(b)中不同湍流模型得到的流場流速和壓力的周向分布基本相同，說明在遠流場各個湍流模型預報的數(shù)值仿真結(jié)果基本相同，但圖7(b)中x/L=0.90處的流速和壓力分布仍有相對明顯的波動，這是由于艇體尾部逆壓力梯度造成尾部流速降低，邊界層增厚，從而改變流場流速和壓力基本對應的變化趨勢.

圖5 流場提取截面

圖6 流場時均軸向速度的周向分布

圖7 Cp的周向分布

據(jù)上推斷，10種湍流模型在Sub-off流場時均變量的捕捉上存在差異，且主要應該體現(xiàn)在流場黏性上，進而導致其他流場變量的變化，在本文采用的時間及空間分辨率的情況下，LES-DE、LES-W同比效果更差，原因可能是空間或時間分辨率不足導致的湍流模型對流場黏性特征的預報存在差異.

不同湍流模型對潛艇周圍流場參數(shù)的平均量預測結(jié)果存在差異，同時對流場的二次量如雷諾應力的預測結(jié)果亦有不同.圖8所示為不同湍流模型仿真得到的潛艇周圍流場的雷諾正應力RNS，其余湍流模型的結(jié)果由于數(shù)值過小，出于文章簡潔的考慮未在文中列出.如同文獻[17]，LES模型對流場二次特征的捕捉較好，DES模型次之，RANS模型最差，故研究流場二次量時模型的選擇會引起數(shù)值仿真結(jié)果的巨大差異.

綜上可知，就模型尺度下的潛艇而言，RANS類湍流模型不能有效地捕捉潛艇流場的流動細節(jié)，對湍流應力有明顯的抑制；DES類湍流模型普遍會抑制潛艇尾流的流動渦結(jié)構(gòu)，同時預報的宏觀力亦存在明顯的偏差；LES類湍流模型在捕捉潛艇流動渦結(jié)構(gòu)、湍流應力方面表現(xiàn)優(yōu)異，但是除LES-S模型外，針對宏觀力的預測均出現(xiàn)明顯的誤差.因此在開展?jié)撏畡恿傲鲌鎏卣鞣治龅木C合研究中，LES-S模型是一個較為優(yōu)異的選擇.

4.3 潛艇周圍渦場特征

基于上述分析得到的數(shù)值仿真策略對Sub-off的流場特征進行計算.由于艇體表面附體的存在，使得其周圍流場較為復雜.圖9所示為潛艇周圍的渦系簡圖，渦系通過Q標準的重新采樣體積進行表達，Q標準是Hunt等[19]提出的一項渦的判斷方法.其中附著在圍殼或尾舵附近的稍渦、馬蹄渦(項鏈渦)會改變螺旋槳的來流均勻性而增加螺旋槳的載荷波動、提高空化風險及增加水動力噪聲，湍流度巨大的渦結(jié)構(gòu)會提高艇體表面激勵力，引發(fā)艇身結(jié)構(gòu)振動，削弱潛艇的隱身性能.艇體附著渦是由邊界層流動引起的速度梯度變化形成，在邊界層較為穩(wěn)定的平行中體段強度微弱，在幾何特征變化明顯的首尾段有較強的艇體附著渦.圖10所示為圍殼附近多截面流場變量，圖中：wi,wj,wk為渦量在3個方向的分量，圖11對比圍殼和尾舵表面的壓力變化.從圖10可以看出，圍殼在軸向首先出現(xiàn)正壓力梯度，在尾段出現(xiàn)負壓力梯度，而其附近的渦量表明稍渦的脫落和壓力梯度的變化存在密切關(guān)系，同時值得注意的是wk可能是稍渦變化的重要變量，存在明顯脫落和分離的wi、wj是圍殼隨邊渦變化產(chǎn)生的重要變量.由于尾舵和圍殼相似的幾何特征，本文認為尾舵稍渦和圍殼稍渦的形成及影響因素相似，但尾舵稍渦的強度明顯大于圍殼稍渦.從圖11可以看出，圍殼和尾舵頂面和側(cè)面的壓力梯度差異可能是導致強度不同的關(guān)鍵因素， 需要注意的是圖10和11中均對變量進行無因次化處理以方便對比分析.文獻[17]的研究表明馬蹄渦的存在與圍殼(尾舵)、艇體邊界層的三維流動有關(guān)，本文不對其進行深入研究.

圖9 Sub-off周圍渦系

圖10 圍殼附近多截面流場變量

圖11 圍殼和尾舵表面Cp分布

5 結(jié)語

本文以Sub-off為幾何模型，針對特定工況，就10種湍流模型計算獲得的結(jié)果進行分析，得出如下結(jié)論：

通過空間及時間收斂性分析，基于理查森外推法判定本文選用的網(wǎng)格方案具有良好的網(wǎng)格收斂性，最終采用了合適的網(wǎng)格開展后續(xù)的數(shù)值計算.本文選用的湍流模型中，LES-S模型在潛艇的宏觀量如阻力、時均流場、伴流方面具有良好的性能或“隨群性”；同時，在流場的二次量預報方面，LES類模型獲得的結(jié)果較其他模型更加豐富.

潛艇渦系因其表面附體的存在較為復雜，但由于圍殼和尾舵的幾何相似性，其渦系具有內(nèi)部形成機理的相似性，但因為附體裝配位置及參數(shù)的差異，尾舵的渦系耦合和演化有較大差異.