李 鵬, 王 超, 孫華偉, 郭春雨
(哈爾濱工程大學(xué) 船舶工程學(xué)院, 哈爾濱 150001)
隨著計(jì)算機(jī)的發(fā)展,計(jì)算流體力學(xué)(Computational Fluid Dynamics, CFD)被越來(lái)越多地應(yīng)用到海洋結(jié)構(gòu)物的性能研究中.潛艇是重要的海洋結(jié)構(gòu)物,準(zhǔn)確了解其不同航速下的水動(dòng)力性能、精確掌握其流場(chǎng)特征,對(duì)提高潛艇自身性能具有重要意義[1].
國(guó)內(nèi)外學(xué)者針對(duì)潛艇性能開(kāi)展過(guò)大量的數(shù)值研究.Pan等[2]采用RANSk-ω(k為湍動(dòng)能,ω為單位湍動(dòng)能耗散率)和雷諾應(yīng)力湍流模型(Reynolds Stress Model, RSM)對(duì)潛艇的阻力進(jìn)行預(yù)報(bào),取得了較好的結(jié)果.Posa等[3-5]基于大渦模擬湍流模型(Large Eddy Simulation, LES)對(duì)潛艇的尾流場(chǎng)進(jìn)行仿真計(jì)算,結(jié)果表明LES能很好地捕捉流動(dòng)特征,同時(shí)為流噪聲的預(yù)報(bào)提供精確輸入.Fu等[6-8]就分離渦湍流模型(Detach Eddy Simulation, DES)下的潛艇水動(dòng)力性能及流場(chǎng)特征進(jìn)行仿真,結(jié)果證實(shí)前者對(duì)垂直壁面方向的無(wú)量綱化距離y+的敏感性微弱,后者對(duì)空間離散精度的高需求性.此外,李士強(qiáng)等[8]還就k-ω、k-ε(ε為湍流耗散率)模型對(duì)潛艇附體變動(dòng)的影響進(jìn)行仿真,并證實(shí)關(guān)于舵失速的問(wèn)題,在失速前k-ε獲取的渦流脫落頻率較k-ω更高.
綜上,國(guó)內(nèi)外學(xué)者就直航潛艇的性能進(jìn)行較多分析,但未闡明不同湍流模型在潛艇水動(dòng)力性能及流場(chǎng)特征計(jì)算中的優(yōu)勢(shì)或短板.本文針對(duì)這一問(wèn)題,首先展多套網(wǎng)格及時(shí)間步長(zhǎng)的收斂性分析.接著,針對(duì)不同的湍流模型開(kāi)展特定工況的數(shù)值仿真計(jì)算,比較預(yù)報(bào)結(jié)果的差異.最后,基于LES-Smagorinsky模型的預(yù)報(bào)結(jié)果對(duì)潛艇渦系特征進(jìn)行分析.
時(shí)均化后不可壓縮流體的連續(xù)性方程與動(dòng)量方程為
(1)
(2)
由于在時(shí)間和空間上對(duì)計(jì)算成本的需求較高,目前主要的CFD求解器均采用求解平均量或?yàn)V波量的方法來(lái)模擬小尺度流動(dòng)以節(jié)省計(jì)算成本.本文就多個(gè)湍流模型開(kāi)展?jié)撏У乃畡?dòng)力性能及流場(chǎng)特征計(jì)算,對(duì)比各種湍流模型的計(jì)算結(jié)果,進(jìn)而篩選適合潛艇水動(dòng)力性能及流場(chǎ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ū)別在于對(duì)額外的黏度項(xiàng)處理方式不同:SA(Spalart-Allmaras)模型是一方程湍流模型,無(wú)壁面函數(shù)要求但增加一個(gè)新的變量SA黏度,其優(yōu)點(diǎn)是對(duì)計(jì)算內(nèi)存要求較低,有良好的收斂性,缺點(diǎn)是不適合求解剪切流與分離流;k-ε模型基于k和ε建立,是具有良好收斂性的兩方程模型,自提出后逐漸演化出realizablek-ε,RNGk-ε等,局限性體現(xiàn)在對(duì)無(wú)滑移壁面、逆壓梯度等的流動(dòng)模擬不精確,耗散不易計(jì)算;k-ω模型計(jì)算結(jié)果對(duì)初始條件較為敏感且收斂相對(duì)困難,但相對(duì)于ε,ω更容易求解.LES是基于自相似理論,通過(guò)過(guò)濾函數(shù)將大渦和小渦分離開(kāi),大渦直接解析,小渦是各項(xiàng)同性,因此可對(duì)其進(jìn)行模化計(jì)算,在解析大渦和?;u的過(guò)程中,產(chǎn)生的亞格子尺度應(yīng)力項(xiàng)通過(guò)亞格子(Subgrid Scale, SGS) 模型計(jì)算.DES模型結(jié)合RANS/LES模型各自的優(yōu)點(diǎn),在近壁面采用RANS進(jìn)行計(jì)算而利用LES模型的算法對(duì)湍流核心區(qū)域計(jì)算,經(jīng)過(guò)多年發(fā)展逐漸演變出延遲分離渦湍流 (Delayed Detach Eddy Simulation, DDES)、改進(jìn)的延遲分離渦湍流 (Improved Delayed Detached Eddy Simulation, IDDES) 等模型變體,對(duì)非定常湍流流動(dòng)的模擬有較大的優(yōu)勢(shì).
本文基于植入到STAR-CCM+中的湍流模型:RANSk-ω、k-ε、SA、LES-Samgorinsky(為行文簡(jiǎn)潔,后續(xù)以LES-S替代)、WALE(LES-W)、Dynamic Smagorinsky(LES-DS)以及DES-SA、DES-EBKE、DDES、IDDES模型開(kāi)展Sub-off水動(dòng)力性能及流場(chǎng)特征的計(jì)算,進(jìn)而甄選合適的湍流模型開(kāi)展?jié)撏畡?dòng)力性能及流場(chǎng)特征的計(jì)算.
本文采用的幾何模型是Sub-off模型,詳細(xì)的模型資料可參考文獻(xiàn)[9-11],后續(xù)的研究中,在原始幾何模型尺度的基礎(chǔ)上進(jìn)行縮放,縮放比λ=0.688 7,模型總長(zhǎng)L=3 m,最大直徑D=0.35 m,若無(wú)特殊說(shuō)明本文的計(jì)算結(jié)果均是基于這一參數(shù)下的幾何模型.
本文設(shè)定笛卡爾坐標(biāo)系Oxyz,其原點(diǎn)位于艇首,如圖1(a)所示,圖中Uref為自由流場(chǎng)的自由來(lái)流速度.本文設(shè)定計(jì)算域?yàn)閳A柱,具體尺寸及邊界條件的設(shè)置見(jiàn)圖1(b),參考文獻(xiàn)[12]的研究,當(dāng)目標(biāo)幾何模型在計(jì)算域流向投影面積小于1%時(shí),數(shù)值計(jì)算結(jié)果可忽略阻塞效應(yīng)的影響,通過(guò)簡(jiǎn)單的數(shù)據(jù)計(jì)算可知,Sub-off模型在計(jì)算域上下底面的投影面積占比為0.7%(小于1%),故可忽略阻塞效應(yīng)的影響.
圖1 幾何模型及計(jì)算域
圖2 網(wǎng)格方案3簡(jiǎn)圖
為統(tǒng)一和方便作者進(jìn)行數(shù)值和試驗(yàn)研究,結(jié)合計(jì)算和實(shí)際條件需求,對(duì)航速為2.53 m/s的潛艇模型的水動(dòng)力性能及流場(chǎng)特征進(jìn)行收斂性計(jì)算,以對(duì)比不同湍流模型獲得的仿真結(jié)果.
表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δ 評(píng)估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)格正交性或湍流模型非線性對(duì)計(jì)算結(jié)果準(zhǔn)確性的影響程度最小[14,18],選用中等密度網(wǎng)格進(jìn)行數(shù)值計(jì)算具有較高可靠性.可知方案2~4具有更好的收斂率且離散精度更接近理論值,故后續(xù)的計(jì)算中采用網(wǎng)格方案3進(jìn)行計(jì)算. 表4是不同時(shí)間步長(zhǎng)計(jì)算得到的判定系數(shù)結(jié)果,以0.001 s的計(jì)算結(jié)果為基準(zhǔn)進(jìn)行對(duì)比.可以看出, 不同時(shí)間步長(zhǎng)對(duì)計(jì)算結(jié)果的影響較小,為提高計(jì)算效率,后續(xù)的計(jì)算中均采用0.001 s為時(shí)間步長(zhǎng). 表4 不同時(shí)間步長(zhǎng)計(jì)算結(jié)果 不同湍流模型對(duì)近壁面流動(dòng)、黏性應(yīng)力等的解析不同,造成仿真得到的阻力有所差異,表5為不同湍流模型計(jì)算得到的潛艇阻力. 表5 不同湍流模型的艇體阻力仿真結(jié)果 表6 不同湍流模型計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果對(duì)比 采用Q標(biāo)準(zhǔn)來(lái)判定流體為渦旋或應(yīng)變主導(dǎo)[19],圖3所示為不同湍流模型仿真計(jì)算得到的Sub-off周圍的渦場(chǎng).圖3清晰地表明LES-W、LES-DS對(duì)潛艇重要渦特征——馬蹄渦的捕捉不足,不同的是除此外的8個(gè)湍流模型均很好地捕捉到了圍殼附近馬蹄渦的流動(dòng)特征,另外這兩種湍流模型預(yù)報(bào)的潛艇尾肩部的渦結(jié)構(gòu)亦沒(méi)有其他模型強(qiáng)烈.選用的10種湍流模型中,RANS類湍流模型對(duì)潛艇尾部渦強(qiáng)度的預(yù)報(bào)結(jié)果顯然低于其他模型,這可能是RANS速度平均后的結(jié)果,是RANS模型的算法特點(diǎn)所決定.值得注意的是,10個(gè)模型均沒(méi)能很好地捕捉或預(yù)報(bào)圍殼稍渦的演化,但對(duì)于尾舵稍渦的演化均能較好地捕捉.推測(cè)原因:一是圍殼稍渦的強(qiáng)度沒(méi)有尾舵稍渦強(qiáng)度大,二是空間分辨率或時(shí)間分辨率不足以捕捉其演化過(guò)程.限于計(jì)算資源以及本文研究重點(diǎn),在此不做深入探討. 圖3 不同湍流模型計(jì)算艇體周圍Q采樣體積(Q=5 s-2) 伴流是考核潛艇綜合性能的重要指標(biāo),對(duì)伴流的準(zhǔn)確預(yù)報(bào)十分重要.圖4所示為不同湍流模型計(jì)算得到的Sub-off尾流伴流分?jǐn)?shù)δ隨距離變化的曲線,圖中,0為潛艇中軸線的最尾端,D、2D、3D、4D為沿著x正方向與潛艇中軸線最尾端的距離.可見(jiàn),LES-DS、LES-W的預(yù)報(bào)結(jié)果與其他模型有明顯偏差,但同屬大渦模擬類中LES-S的預(yù)報(bào)結(jié)果有良好的“隨群”性,結(jié)合對(duì)于宏觀力的預(yù)報(bào)結(jié)果可知,LES-S模型在預(yù)報(bào)潛艇宏觀變量方面有著較好的適用性. 圖4 不同湍流模型預(yù)報(bào)尾流伴流(0位于艇體尾端) 為了深入探究不同湍流模型對(duì)Sub-off流場(chǎng)的預(yù)報(bào)結(jié)果, 本文提取了如圖5虛線所示位置的流場(chǎng)時(shí)均變量并加以分析,時(shí)均流場(chǎng)的采樣總長(zhǎng)是流經(jīng)一倍艇長(zhǎng)時(shí)間. 圖6所示為時(shí)均軸向速度的周向變化曲線,圖中θ為圖5中提取截面沿著白色虛線環(huán)向位置的角度,u為流場(chǎng)局部軸向速度,通過(guò)u/Uref對(duì)軸向速度進(jìn)行無(wú)因次化對(duì)比, 分析其變化規(guī)律,為了有效說(shuō)明流場(chǎng),將Sub-off周圍流場(chǎng)分為內(nèi)流場(chǎng)(r/D<0.501)和外流場(chǎng)(r/D>0.501),其中r為對(duì)應(yīng)曲線的半徑值.可知,10種湍流模型對(duì)Sub-off的外流場(chǎng)及其平行中體內(nèi)流場(chǎng)的預(yù)報(bào)結(jié)果吻合程度非常高,推測(cè)由于遠(yuǎn)離壁面或壁面流動(dòng)簡(jiǎn)單,各湍流模型對(duì)其處理無(wú)明顯差異,但即使對(duì)平行中體內(nèi)流場(chǎng)的流動(dòng)預(yù)報(bào)差異微弱,Sub-off尾內(nèi)流場(chǎng)的預(yù)報(bào)差異卻十分明顯.這表明,湍流模型對(duì)因?yàn)閹缀翁卣髯兓l(fā)流動(dòng)變化更為敏感,對(duì)流場(chǎng)“傳遞過(guò)程”引發(fā)的變化敏感性較弱.圖7所示為流場(chǎng)時(shí)均壓力的周向分布曲線,可以看出,Sub-off周圍流場(chǎng)的時(shí)均壓力系數(shù)Cp(Cp=p/(0.5ρUref))分布曲線對(duì)湍流模型地變化不敏感,各個(gè)模型均有良好的“隨群”性,由于尾部附體以及尾部幾何收縮造成的負(fù)壓力梯度,在x/L=0.9處,流場(chǎng)壓力波動(dòng)仍較為劇烈.圖6和7的結(jié)果表明:時(shí)均速度較時(shí)均壓力對(duì)湍流模型的改變更敏感,尤其是近壁面流動(dòng),如圖6(a)、7(a)所示,這是不同湍流模型調(diào)用不同壁面解析函數(shù)造成的.圖6(b)、7(b)中不同湍流模型得到的流場(chǎng)流速和壓力的周向分布基本相同,說(shuō)明在遠(yuǎn)流場(chǎng)各個(gè)湍流模型預(yù)報(bào)的數(shù)值仿真結(jié)果基本相同,但圖7(b)中x/L=0.90處的流速和壓力分布仍有相對(duì)明顯的波動(dòng),這是由于艇體尾部逆壓力梯度造成尾部流速降低,邊界層增厚,從而改變流場(chǎng)流速和壓力基本對(duì)應(yīng)的變化趨勢(shì). 圖5 流場(chǎng)提取截面 圖6 流場(chǎng)時(shí)均軸向速度的周向分布 圖7 Cp的周向分布 據(jù)上推斷,10種湍流模型在Sub-off流場(chǎng)時(shí)均變量的捕捉上存在差異,且主要應(yīng)該體現(xiàn)在流場(chǎng)黏性上,進(jìn)而導(dǎo)致其他流場(chǎng)變量的變化,在本文采用的時(shí)間及空間分辨率的情況下,LES-DE、LES-W同比效果更差,原因可能是空間或時(shí)間分辨率不足導(dǎo)致的湍流模型對(duì)流場(chǎng)黏性特征的預(yù)報(bào)存在差異. 不同湍流模型對(duì)潛艇周圍流場(chǎng)參數(shù)的平均量預(yù)測(cè)結(jié)果存在差異,同時(shí)對(duì)流場(chǎng)的二次量如雷諾應(yīng)力的預(yù)測(cè)結(jié)果亦有不同.圖8所示為不同湍流模型仿真得到的潛艇周圍流場(chǎng)的雷諾正應(yīng)力RNS,其余湍流模型的結(jié)果由于數(shù)值過(guò)小,出于文章簡(jiǎn)潔的考慮未在文中列出.如同文獻(xiàn)[17],LES模型對(duì)流場(chǎng)二次特征的捕捉較好,DES模型次之,RANS模型最差,故研究流場(chǎng)二次量時(shí)模型的選擇會(huì)引起數(shù)值仿真結(jié)果的巨大差異. 綜上可知,就模型尺度下的潛艇而言,RANS類湍流模型不能有效地捕捉潛艇流場(chǎng)的流動(dòng)細(xì)節(jié),對(duì)湍流應(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模型外,針對(duì)宏觀力的預(yù)測(cè)均出現(xiàn)明顯的誤差.因此在開(kāi)展?jié)撏畡?dòng)力及流場(chǎng)特征分析的綜合研究中,LES-S模型是一個(gè)較為優(yōu)異的選擇. 基于上述分析得到的數(shù)值仿真策略對(duì)Sub-off的流場(chǎng)特征進(jìn)行計(jì)算.由于艇體表面附體的存在,使得其周圍流場(chǎng)較為復(fù)雜.圖9所示為潛艇周圍的渦系簡(jiǎn)圖,渦系通過(guò)Q標(biāo)準(zhǔn)的重新采樣體積進(jìn)行表達(dá),Q標(biāo)準(zhǔn)是Hunt等[19]提出的一項(xiàng)渦的判斷方法.其中附著在圍殼或尾舵附近的稍渦、馬蹄渦(項(xiàng)鏈渦)會(huì)改變螺旋槳的來(lái)流均勻性而增加螺旋槳的載荷波動(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所示為圍殼附近多截面流場(chǎng)變量,圖中:wi,wj,wk為渦量在3個(gè)方向的分量,圖11對(duì)比圍殼和尾舵表面的壓力變化.從圖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中均對(duì)變量進(jìn)行無(wú)因次化處理以方便對(duì)比分析.文獻(xiàn)[17]的研究表明馬蹄渦的存在與圍殼(尾舵)、艇體邊界層的三維流動(dòng)有關(guān),本文不對(duì)其進(jìn)行深入研究. 圖9 Sub-off周圍渦系 圖10 圍殼附近多截面流場(chǎng)變量 圖11 圍殼和尾舵表面Cp分布 本文以Sub-off為幾何模型,針對(duì)特定工況,就10種湍流模型計(jì)算獲得的結(jié)果進(jìn)行分析,得出如下結(jié)論: 通過(guò)空間及時(shí)間收斂性分析,基于理查森外推法判定本文選用的網(wǎng)格方案具有良好的網(wǎng)格收斂性,最終采用了合適的網(wǎng)格開(kāi)展后續(xù)的數(shù)值計(jì)算.本文選用的湍流模型中,LES-S模型在潛艇的宏觀量如阻力、時(shí)均流場(chǎng)、伴流方面具有良好的性能或“隨群性”;同時(shí),在流場(chǎng)的二次量預(yù)報(bào)方面,LES類模型獲得的結(jié)果較其他模型更加豐富. 潛艇渦系因其表面附體的存在較為復(fù)雜,但由于圍殼和尾舵的幾何相似性,其渦系具有內(nèi)部形成機(jī)理的相似性,但因?yàn)楦襟w裝配位置及參數(shù)的差異,尾舵的渦系耦合和演化有較大差異.4 結(jié)果與討論
4.1 水動(dòng)力性能
4.2 流場(chǎng)特征
4.3 潛艇周圍渦場(chǎng)特征
5 結(jié)語(yǔ)