(深海載人裝備國家重點實驗室,江蘇 無錫 214082)
新型水下潛器外形各異,在設計研發(fā)過程中水動力設計無母型參考,潛艇相關(guān)規(guī)范也不適用,只能依靠水池試驗和CFD仿真。在整個研發(fā)周期內(nèi),尤其是在初期多方案選型時,動穩(wěn)性評估都是必不可少的關(guān)鍵環(huán)節(jié)。盡管旋臂水池試驗能提供穩(wěn)定性評估需要的旋轉(zhuǎn)導數(shù),但是,試驗周期長,費用高,常常無法及時地用于初期的線型優(yōu)化設計中。國內(nèi)對旋臂水池試驗的數(shù)值仿真方法主要采用多參考系、動量源項法和滑移網(wǎng)格法。前兩種都是采用穩(wěn)態(tài)的方法來模擬旋臂水池中物體的運動過程,而滑移網(wǎng)格法采用瞬態(tài)模擬,理論上與物理實際的相近程度更高。雖然該方法計算量較大,但是,當前計算能力已經(jīng)達到了較高的水平,計算能力已經(jīng)不是突出的問題。因此,考慮采用某商業(yè)軟件使用滑移網(wǎng)格法對suboff標模進行數(shù)值仿真。對旋臂水池的仿真研究[1-3]未見具體的數(shù)值計算設置。因此,考慮對計算域大小、邊界層厚度、湍流模式、臨界雷諾數(shù)等因素對仿真的影響進行系統(tǒng)性分析,分析探討較為準確的旋臂水池數(shù)值模擬方法。
采用的計算模型為美國國防高級研究計劃署(DARPA)[4-5]提出的SUBOFF AFF-1和AFF-8。SUBOFF主體是水滴形的軸對稱體,全附體還包括指揮臺圍殼和尾翼,見圖1。
圖1 SUBOFF模型
采用潛艇操縱性規(guī)范中規(guī)定坐標系、量綱一的量化及表達方法。為方便與試驗結(jié)果對比,參照文獻中的規(guī)定,取特征長度L=4.261 m。坐標系原點位于浮心位置處。浮心位于主體軸線上,距艏2.013 m。計算中的工況參數(shù)也與試驗中的相同,水溫18℃,密度ρ=998.55 kg/m3,動力黏度取μ=0.001 053 kg/(m·s)。三向力和力矩系數(shù)分別以Cx、Cy、Cz、Cmx、Cmy和Cmz表達,均以特征長度L、旋轉(zhuǎn)線速度V量綱一的量化。無因次角速度r′=rL/V。
不可壓縮黏性流體的連續(xù)性方程和RANS方程可寫成如下形式。
(1)
(2)
式(1)、(2)是不封閉的,需要尋求補充關(guān)系(湍流模型)使問題封閉,采用realizablek-ε或SSTk-ω湍流模型進行湍流封閉[6]。
為減小網(wǎng)格數(shù)量,將整個環(huán)形計算區(qū)域分割成2個同心的圓環(huán),見圖2。大圓環(huán)為靜止區(qū)域,內(nèi)部小圓環(huán)為旋轉(zhuǎn)區(qū)域。艇體處于小圓環(huán)內(nèi),跟隨內(nèi)部圓環(huán)圍繞軸心旋轉(zhuǎn),實現(xiàn)艇體的旋轉(zhuǎn)模擬,動靜區(qū)域通過內(nèi)外圓環(huán)交界面?zhèn)鬟f數(shù)據(jù)。旋轉(zhuǎn)半徑分別取10、13、15、20和25 m。
圖2 計算域及計算域參數(shù)示意
采用非結(jié)構(gòu)化網(wǎng)格劃分方法對光體和全附體線型進行網(wǎng)格劃分,見圖3。在邊界層附近劃分棱柱層網(wǎng)格,并對翼、圍殼、艇體周圍等區(qū)域進行局部加密。全附體計算時,網(wǎng)格總數(shù)約200萬。
圖3 圓環(huán)計算域網(wǎng)格劃分
采用有限體積法(FVM)離散控制方程和湍流模型。對流項采用二階迎風格式,擴散項采用中心差分格式。對時間項采用一階離散格式,每個時間步旋轉(zhuǎn)0.5°,每步內(nèi)部迭代6次。對壓力速度方程采用分離式計算方法。
直航仿真中,為了減小遠場邊界的干擾,通常需要將艇體距離四周邊界2倍艇長以上。而旋臂水池仿真中艇體旋轉(zhuǎn)導致流動更加復雜。為了確定合適的邊界條件,分別取外部大圓環(huán)的半徑差Ro等于L、1.5L和2.0L進行計算,結(jié)果見表1。
表1 Ro對光體和全附體仿真結(jié)果的影響
從表1可以看到,增大Ro對光體和全附體的所有系數(shù)的影響非常微弱??梢奟o超過特征長度1倍時已經(jīng)滿足計算要求。
滑移網(wǎng)格法中,艇體周圍被小的圓環(huán)包裹。為探討旋轉(zhuǎn)內(nèi)部計算域?qū)Ψ抡娼Y(jié)果的影響,分別取內(nèi)部圓環(huán)半徑Ri為2H(H為艇高)和3H對光體和全附體進行計算,結(jié)果見圖4。
圖4 Ri對光體和全附體計算結(jié)果的影響
由圖4可見,Ri的增加對光體和全附體的Cy和Cmz幾乎沒有影響,2H已經(jīng)滿足計算需求。盡管Ri對Cx的計算有一定影響,但是程度較小,不足0.4%。此外,Cx不用于求線性旋轉(zhuǎn)導數(shù),是非重點關(guān)注量。增大Ri會大幅增加網(wǎng)格數(shù)量,計算效率大減。而當Ri從2H繼續(xù)減小時,會導致艇體與邊界距離太短,布置網(wǎng)格數(shù)量較少,容易降低網(wǎng)格質(zhì)量。因此,綜合認為旋臂水池仿真中Ri超過2H即可。
根據(jù)直航阻力計算時的經(jīng)驗,通常第一層邊界層厚度對湍流的計算有一定影響。為此,采用realizablek-ε湍流模式,對第一層網(wǎng)格厚度分別取1.2、2.4、3.6和4.8 mm,相對應的y+分別為40、80、120和160。并都取6層棱柱層網(wǎng)格,增長率默認為1.3,其余設置均相同。計算結(jié)果見圖5。
圖5 y+對光體和全附體計算的影響
由圖5可見,不同y+下,光體和全附體的力和力矩系數(shù)幾乎相同,表明其影響均較小??梢?,在旋臂水池的模擬中,采用realizablek-ε湍流模式下y+對力和力矩的影響較小,取值在40~160內(nèi)均可。
旋臂試驗中,通常需要先進行臨界雷諾數(shù)測量試驗,而臨界雷諾數(shù)通常是通過水動力系數(shù)的穩(wěn)定來判斷,不同標準可能會導致其數(shù)值有一定的變化。
由上節(jié)結(jié)論可知,采用realizablek-ε時,y+在合適的范圍內(nèi),對仿真結(jié)果影響較小。因此,在不同雷諾數(shù)計算時,僅需調(diào)整邊界層網(wǎng)格厚度保持相當?shù)膟+即可。對不同雷諾數(shù)下的工況進行仿真時,均采用realizablek-ε湍流模式,網(wǎng)格數(shù)量基本相同,僅對邊界層網(wǎng)格調(diào)整保持y+基本相同即可。從圖6計算結(jié)果可以看到,光體和全附體下Cx、Cy和Cmz都隨著線速度增加逐漸平穩(wěn)。Cx和Cy在超過4 m/s后變化幅度較小,初步可以認為4 m/s已經(jīng)到達臨界雷諾數(shù),此時雷諾數(shù)約17×106,比文獻[5]中的14×106稍大。
圖6 雷諾數(shù)對光體和全附體計算的影響
對臨界雷諾數(shù)附近不同速度進行線性擬合得到旋轉(zhuǎn)導數(shù)見表2。其中,光體的擬合結(jié)果見圖7。
圖7 光體線性擬合求旋轉(zhuǎn)導數(shù)
表2 不同線速度下旋轉(zhuǎn)導數(shù)計算結(jié)果對比
由表2可見,線速度對旋轉(zhuǎn)導數(shù)的計算有一定的影響,對光體的影響大于對全附體,影響在5%以內(nèi),對全附體的影響則在2%以內(nèi)。從數(shù)值上,隨著線速度的增加,光體旋轉(zhuǎn)導數(shù)則逐漸減小,全附體則逐漸增大。在超過4 m/s時,線速度對旋轉(zhuǎn)導數(shù)的影響可以忽略。
采用上文研究得到的計算設置,取文獻中的線速度6.5 kn,采用不同的旋轉(zhuǎn)半徑對光體和全附體進行水平面的旋轉(zhuǎn)仿真。采用2種湍流模式,以對比其影響,見表3。
表3 CFD仿真結(jié)果與試驗結(jié)果對比
從表3中可以看到,realizablek-ε和SSTk-ω模擬得到N′r結(jié)果幾乎完全相同,而Y′r則有一定的不同,前者與試驗結(jié)果符合更好。光體模型的計算結(jié)果較全附體差,這是由于光體外形光順,分離點較不固定的原因。分離點的變動對艇體流場和受力產(chǎn)生較大影響,對分離點的準確求取較為困難。
以Suboff為研究對象,通過仿真對比,發(fā)現(xiàn)計算內(nèi)部旋轉(zhuǎn)域縱向與艇體距離大于2倍艇高,外部計算域縱向與艇體距離了超過2倍艇長時,計算域大小對計算結(jié)果的影響可以忽略。y+在40~160時,對仿真結(jié)果的影響不顯著。雷諾數(shù)對旋轉(zhuǎn)導數(shù)的計算有一定影響,且對光體的影響較全附體大,必須超過臨界雷諾數(shù)以消除其影響。realizablek-ε和SSTk-ω兩種湍流模式下得到的N′r幾乎完全相同,而Y′r有一定的不同??傮w來看前者與試驗符合程度更高。全附體N′r模擬結(jié)果與試驗結(jié)果符合非常好,其余系數(shù)誤差基本在30%以內(nèi)。盡管以上因素會對計算結(jié)果產(chǎn)生一定的影響,但是,波動較小,顯現(xiàn)出較強的可靠性,可以采用該方法用于初期的方案優(yōu)選。