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

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

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

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

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

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

1 理論基礎(chǔ)

1.1 控制方程

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

(1)

(2)

1.2 湍流模型

由于在時(shí)間和空間上對計(jì)算成本的需求較高，目前主要的CFD求解器均采用求解平均量或?yàn)V波量的方法來模擬小尺度流動(dòng)以節(jié)省計(jì)算成本.本文就多個(gè)湍流模型開展?jié)撏У乃畡?dòng)力性能及流場特征計(jì)算，對比各種湍流模型的計(jì)算結(jié)果，進(jìn)而篩選適合潛艇水動(dòng)力性能及流場特征仿真計(jì)算的湍流模型.

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

2 仿真基礎(chǔ)

2.1 坐標(biāo)系和計(jì)算域

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

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

圖1 幾何模型及計(jì)算域

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

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

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

3 收斂性分析

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

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

表3 網(wǎng)格收斂性計(jì)算結(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.理論上，當(dāng)P=2時(shí)，網(wǎng)格正交性或湍流模型非線性對計(jì)算結(jié)果準(zhǔn)確性的影響程度最小[14,18]，選用中等密度網(wǎng)格進(jìn)行數(shù)值計(jì)算具有較高可靠性.可知方案2～4具有更好的收斂率且離散精度更接近理論值，故后續(xù)的計(jì)算中采用網(wǎng)格方案3進(jìn)行計(jì)算.

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

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

4 結(jié)果與討論

4.1 水動(dòng)力性能

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

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

表6 不同湍流模型計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果對比

4.2 流場特征

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

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

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

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

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

圖5 流場提取截面

圖6 流場時(shí)均軸向速度的周向分布

圖7 Cp的周向分布

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

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

綜上可知，就模型尺度下的潛艇而言，RANS類湍流模型不能有效地捕捉潛艇流場的流動(dòng)細(xì)節(jié)，對湍流應(yīng)力有明顯的抑制；DES類湍流模型普遍會(huì)抑制潛艇尾流的流動(dòng)渦結(jié)構(gòu)，同時(shí)預(yù)報(bào)的宏觀力亦存在明顯的偏差；LES類湍流模型在捕捉潛艇流動(dòng)渦結(jié)構(gòu)、湍流應(yīng)力方面表現(xiàn)優(yōu)異，但是除LES-S模型外，針對宏觀力的預(yù)測均出現(xiàn)明顯的誤差.因此在開展?jié)撏畡?dòng)力及流場特征分析的綜合研究中，LES-S模型是一個(gè)較為優(yōu)異的選擇.

4.3 潛艇周圍渦場特征

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

圖9 Sub-off周圍渦系

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

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

5 結(jié)語

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

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

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