魏可可,高霄鵬,馬 騁,董祖舜
(1. 湖北工程學(xué)院 機械工程學(xué)院, 湖北 孝感 432000; 2. 海軍工程大學(xué) 艦船與海洋學(xué)院, 湖北 武漢 430033; 3. 海軍研究院, 北京 100161)
作為現(xiàn)代國防力量的重要組成部分,潛艇在出海遠航及戰(zhàn)略威懾中均發(fā)揮著重要的作用,然而潛艇應(yīng)急上浮是一個極其危險的運動,潛艇從上浮到出水對潛艇穩(wěn)性和艇中人員安全有著很大的影響[1],準(zhǔn)確預(yù)報潛艇上浮運動是探索研究潛艇上浮運動控制規(guī)律以及潛艇安全上浮的前提。好多學(xué)者以SUBOFF標(biāo)準(zhǔn)模型為研究對象[2],對其上浮、操縱、直航等運動開展過研究[3];Carrica等[4-5]采用CFD Ship-Iowa V4對潛艇六自由度運動進行了的數(shù)值模擬,預(yù)報了不同工況下潛艇的操縱運動;基于真實槳和耦合虛擬槳的計算方法,Martin等[6]對潛艇的操縱運動進行了模擬,并對比分析了兩種方法的計算結(jié)果;Chase等[7]基于雷諾平均納維-斯托克斯(Reynolds average Navier-Stokes,RANS)方法,以SUBOFF為研究對象,網(wǎng)格劃分采用重疊網(wǎng)格,對潛艇的直航和自航運動開展了仿真模擬;周廣禮等[8-10]基于RANS方法,采用動網(wǎng)格,模擬了SUBOFF潛艇的上浮-出水運動,并監(jiān)測了潛艇姿態(tài)變化;孫斌等[11]和孫濟政等[12]對潛艇大型集體逃生艙做了上浮優(yōu)化方案,基于CFD對逃生艙的上浮開展了仿真模擬,對逃生艙的上浮速度進行了分析;錢恒等[13]依據(jù)客觀事實和相關(guān)實例,構(gòu)建潛艇極地破冰上浮自然環(huán)境風(fēng)險綜合評估指標(biāo)體系,并提出基于云變換-三角模糊數(shù)隨機模擬的模糊層次分析法風(fēng)險評估模型;呂幫俊等[14]基于潛艇操縱運動仿真數(shù)學(xué)模型,對高壓氣吹除主壓載水艙后潛艇應(yīng)急上浮過程的運動狀態(tài)進行了預(yù)報;閆朋[15]基于CFD軟件STARCCM+提出了一套有效的數(shù)值計算方法,通過對模型試驗結(jié)果與數(shù)值計算結(jié)果的分析,總結(jié)了初始航速、穩(wěn)心高、浮力載荷作用位置和艉舵角對潛艇應(yīng)急上浮運動的影響規(guī)律。
以上對潛艇上浮運動的研究主要局限于數(shù)值計算,對其相關(guān)試驗研究和數(shù)值方法的驗證較少。本文以全附體的SUBOFF(指揮臺含有舵)為研究對象,對其上浮運動的數(shù)值計算方法開展了研究,并與模型上浮試驗進行了對比驗證,證明了數(shù)值計算方法的可行性。
SUBOFF是由美國泰勒研究所提供的潛艇試驗?zāi)P?因為該模型是標(biāo)準(zhǔn)模型,所以進行過許多試驗,在原SUBOFF模型基礎(chǔ)上,在其指揮臺圍殼上增加了舵,修改的SUBOFF模型如圖1所示,表1則給出了模型相關(guān)參數(shù)。
表1 模型主參數(shù)
圖1 艇體三維模型Fig.1 Three dimensional model of hull
為了便于對潛艇上浮運動的分析和對相關(guān)結(jié)果參數(shù)的監(jiān)測,根據(jù)國際水池會議(ITTC)發(fā)布的權(quán)威標(biāo)準(zhǔn),采用E-ξηζ為大地坐標(biāo)系,原點E可選地球上某一定點,如海中的任一點;G-xyz為隨艇坐標(biāo)系,隨艇坐標(biāo)系的原點位于艇的重心G處。潛艇的上浮的運動可分解為水平面運動和垂直面運動,因此其坐標(biāo)系也可分為兩種,如圖2所示。
(a) 水平面運動坐標(biāo)系(a) Coordinate system of horizontal motion
漂角β為航速V與x軸之間的夾角;α為來流的攻角;V在隨艇坐標(biāo)系上的投影為(u,v,w);潛艇的角速度Ω在隨艇坐標(biāo)系上的投影為(p,q,r);潛艇所受的外力F在隨艇坐標(biāo)系上的投影為(Fx,Fy,Fz);潛艇所受的力矩Q可用隨艇坐標(biāo)系上(K,M,N)來表示;潛艇的位置坐標(biāo)可用隨艇坐標(biāo)系原點的地面坐標(biāo)值(ξG,ηG,ζG)來表示;潛艇的姿態(tài)角可用相對于大地坐標(biāo)系的三個姿態(tài)角(φ,θ,ψ)來表示,其中φ為橫傾角、θ為縱傾角、ψ為航向角。
潛艇的上浮到出水運動是復(fù)雜氣液耦合黏性流體運動,對其開展數(shù)值模擬時,主要以連續(xù)性方程和動量方程為基礎(chǔ),方程的形式為:
(1)
其中,ρ為流體密度,xi、xj為i、j方向的坐標(biāo),ui、uj為流體速度沿i、j方向的分量,p為靜壓,fi為單位質(zhì)量力。
對于潛艇上浮運動的模擬,基于重疊網(wǎng)格技術(shù),計算域主要由背景域和重疊網(wǎng)格域組成,重疊網(wǎng)格域位于背景域里面,重疊網(wǎng)格域是跟隨潛艇一起運動的,為避免重疊網(wǎng)格域隨艇體運動移出背景域所設(shè)定的加密區(qū),需將艇體在平面內(nèi)的運動投影至背景域,即通過自編程序?qū)撏г诖蟮刈鴺?biāo)系下的水平移動速度和垂直移動速度賦予背景域,背景域隨艇體移動。不同時刻潛艇上浮運動重疊網(wǎng)格應(yīng)用情況如圖3所示,潛艇運動過程中,在水平面內(nèi)重疊網(wǎng)格域與背景域不發(fā)生相對移動,僅存在相對轉(zhuǎn)動;而在垂直面內(nèi),重疊網(wǎng)格域相對整體域可垂向移動、首搖、縱搖及橫搖。
(a) 初始時刻(a) Initial moment
(b) 水面時刻(b) Moment of water surface
計算域采用方形域,尺寸大小為9L×4L×4.8L,潛艇上浮運動的計算域和艇體的網(wǎng)格劃分如圖4所示。
圖4 SUBOFF模型上浮運動的計算域及艇體的網(wǎng)格Fig.4 Grid of computational domain and hull of SUBOFF model buoyancy motion
為了更好地模擬潛艇上浮運動,需要湍流邊界層,而湍流模型以k-ω模型為基準(zhǔn),采用黏流方法開展?jié)撏细∵\動數(shù)值模擬時,不同的湍流模型對計算的精度是不同的,眾多學(xué)者已針對該問題開展了大量研究,但所得的結(jié)論存在較大差異。本文主要研究Standardk-ω和 剪切應(yīng)力輸運(shear stress transport,SST)(Menter)k-ω湍流模型對潛艇上浮運動參數(shù)結(jié)果的影響,對其開展優(yōu)化分析。
Standardk-ω湍流模型方程為:
(2)
SST(Menter)k-ω湍流模型方程為:
(3)
其中,V為體積,μ為動力黏度,σk和σω為模型系數(shù),fβ*為自由剪切修正因子,fβ為渦流延伸修正因子,Sk和Sω為源項,k表示湍動能,ω為特殊湍動能耗散,其余參數(shù)具體形式及參數(shù)值參見文獻[16]。
潛艇的上浮運動和水面直航運動都涉及湍流邊界層的運動,由于水面直航更加易于對潛艇拖曳阻力的測量,且湍流模式對潛艇水動力性能及運動形式并無影響,因此基于拖曳阻力試驗,對Standardk-ω湍流模型和SST(Menter)k-ω湍流模型進行優(yōu)選,為后續(xù)潛艇上浮運動的數(shù)值研究奠定基礎(chǔ)。為此,分別以SUBOFF全附體模型和裸艇體模型為研究對象,對其展開靜水中阻力的數(shù)值計算,并與試驗值進行對比,模型的三維圖如圖5所示,對比結(jié)果如圖6所示。
(a) 全附體模型(a) Model of full appendage
(a) 全附體模型結(jié)果(a) Results of the full appendage model
圖6給出了SUBOFF全附體模型和裸艇體模型在航速分別為5.93 kn、10.00 kn、11.85 kn、13.92 kn、16.00 kn、17.79 kn時Standardk-ω湍流模型和SST(Menter)k-ω湍流模型下的拖曳阻力數(shù)值計算結(jié)果,同時與試驗值進行了對比。從圖6中可知,無論是全附體模型還是裸艇體模型,在低航速階段,兩種湍流模型的數(shù)值計算結(jié)果與試驗結(jié)果基本一致,但在中高速階段,Standardk-ω湍流模型的數(shù)值計算結(jié)果精度更高。由此可知,Standardk-ω模型可作為優(yōu)選的湍流模型。
在網(wǎng)格上對離散的值建立差分方程的方法稱為差分格式,差分方程的本質(zhì)就是連續(xù)空間上的微分方程。當(dāng)界面的物理量需要考慮方向時,對網(wǎng)格的離散可采用迎風(fēng)格式,潛艇的上浮運動是一個多維多方向的運動,因此在數(shù)值模擬中采用迎風(fēng)格式對網(wǎng)格進行離散。
以全附體的SUBOFF為研究對象,對其開展了一階和二階迎風(fēng)格式的上浮運動數(shù)值計算,分別計算了全附體的SUBOFF模型上浮運動到出水過程中的縱向力、橫向力、垂向力、縱向力矩、橫向力矩以及垂向力矩,如圖7所示。從圖中可知,兩種迎風(fēng)格式下潛艇上浮運動的力和力矩的計算結(jié)果差別較小,即兩者的計算精度相當(dāng),但考慮實際的計算時間,二階迎風(fēng)格式要優(yōu)于一階迎風(fēng)格式,同時,基于一階迎風(fēng)格式,二階迎風(fēng)格式包含了物理量在節(jié)點間的分布曲線,因此也考慮了曲線曲率變化的影響。綜上,潛艇上浮運動的數(shù)值計算優(yōu)選二階迎風(fēng)格式。
(a) 縱向力(a) Longitudinal force
為了驗證數(shù)值計算方法,須開展?jié)撏У哪P蜕细≡囼灐DP蜕细≡囼炛饕揽可细≡囼炂脚_來進行,潛艇的上浮是從水下到水面的六自由度運動,因此該試驗平臺必須含有能固定模型及釋放模型的裝置。同時為了保證在試驗前,潛艇能下潛到指定的水深,因此該試驗裝置應(yīng)含有可以承載潛艇模型的平臺并能調(diào)節(jié)不同的潛深。該試驗平臺包含的基本裝置有:平臺底盤(內(nèi)置水密發(fā)動機)、支撐架、拖車平臺及模型、升降支座、釋放器、軌道、升降裝置、水下高清攝像機等設(shè)備,試驗平臺主要構(gòu)件的尺寸如表2所示。
表2 平臺主要結(jié)構(gòu)裝置的尺度及質(zhì)量
整個潛艇模型上浮系統(tǒng)主要由岸機系統(tǒng)和上浮試驗平臺組成,其中,岸機系統(tǒng)主要由監(jiān)測系統(tǒng)和控制系統(tǒng)組成,監(jiān)測系統(tǒng)通過零浮力光纖和數(shù)據(jù)線分別和艇體以及高清水下攝像機相連,負責(zé)實時監(jiān)測艇的傾角、加速度、水深等參數(shù)的時歷變化以及實時記錄艇體水下運動姿態(tài);潛艇模型上浮試驗系統(tǒng)如圖8所示。
圖 8 潛艇模型上浮系統(tǒng)示意圖Fig.8 Schematic diagram of submarine model uplift system
為了精確完成潛艇應(yīng)急上浮模型試驗,試驗時采用機載式數(shù)據(jù)采集記錄系統(tǒng)進行數(shù)據(jù)的采集與記錄,采用姿態(tài)儀對潛艇的橫傾、縱傾及艏搖等姿態(tài)角進行測量,采用深度傳感器對潛艇模型的水深進行實時測量,并在艇體的艏、舯、艉三個位置安裝了三軸加速度傳感器。
機載式數(shù)據(jù)采集記錄系統(tǒng)、姿態(tài)儀及加速度傳感器置于模型內(nèi)部密封艙體內(nèi)部,深度傳感器固定在模型密封艙蓋外端面,艇體上浮試驗儀器設(shè)備的具體布置位置及功能如下:
1)姿態(tài)儀固定于重心處,可以實時測量艇體上浮過程中的姿態(tài)角;加速度傳感器3個,均為3軸加速度傳感器,分別固定于艇體軸線上艏、舯、艉三個不同縱向位置處,以測量艇體上浮到出水全過程中的運動狀態(tài)。
2)水深傳感器固定在模型密封艙蓋外端面與水接觸,可以實時測量模型上浮過程中的潛深位置;機載式數(shù)據(jù)采集記錄系統(tǒng)置于艇體內(nèi)部,可以直接采集姿態(tài)儀、加速度、深度傳感器的測量數(shù)據(jù)。
3)岸上布置的上位機主要采集數(shù)據(jù)和實時監(jiān)測艇體的姿態(tài)角、加速度、水深等參數(shù)的時歷結(jié)果。
4)水下攝像機主要拍攝艇在不同工況下艇體上浮到出水的運動軌跡和運動姿態(tài),并實時傳回岸機。
在開展模型上浮試驗前,根據(jù)需開展的試驗工況對模型狀態(tài)進行調(diào)試,模型調(diào)試好在水下拖車平臺上進行安裝和固定,在岸機上對整個測試系統(tǒng)的可靠性進行檢查,確保系統(tǒng)在運行時無誤;對模型的艙蓋進行密封,對模型進行氣密性的檢查;對上浮試驗平臺通電,將模型下潛至所需的水深位置,啟動水下拖車平臺,釋放模型上浮;在上浮試驗過程中,實時對模型數(shù)據(jù)進行采集并對模型的軌跡及姿態(tài)錄像;完成上浮試驗后,對模型進行打撈回收,更換模型狀態(tài),進行下一個工況的模型試驗準(zhǔn)備。
根據(jù)不確定度的相關(guān)計算方法[17],對潛艇上浮運動的計算結(jié)果開展不確定性分析,研究所構(gòu)建的數(shù)值計算方法是否可行。驗證和確認是開展不確定性分析的重要的兩個流程:驗證就是對數(shù)值誤差的不確定度USN開展研究,迭代次數(shù)、網(wǎng)格劃分、時間步長等因素都會導(dǎo)致數(shù)值誤差的發(fā)生;確認就是對數(shù)值的不確定度USN開展研究,主要對誤差E和確認不確定度UV來開展對比分析[18-20]。
對于潛艇上浮數(shù)值模擬進行驗證,主要是針對網(wǎng)格收斂性和時間步長收斂性開展的。本文主要對不同正上浮力下潛艇零速自由上浮的數(shù)值計算結(jié)果進行網(wǎng)格收斂性和時間步長收斂性開展研究。
表3 不同網(wǎng)格類型
以不同正上浮力占比下的垂向速度為考量對象,表4和表5分別給出了網(wǎng)格收斂性結(jié)果和時間步長收斂性結(jié)果。
表4 網(wǎng)格收斂性驗證結(jié)果
表5 時間步長收斂性驗證結(jié)果
(4)
其中:ε32為粗、中網(wǎng)格數(shù)值結(jié)果之差,ε32=S3G-S2G;ε21為中、細網(wǎng)格數(shù)值結(jié)果之差,ε21=S2G-S1G;Pkest為數(shù)值計算中首項準(zhǔn)確度極限階數(shù)的估算值,參考文獻[19]取2。
表6 數(shù)值計算的確認結(jié)果
由表6可知,數(shù)值計算結(jié)果與試驗值的對比誤差|E|均小于確認不確定度UV,這說明三種工況下數(shù)值計算的結(jié)果都得到了確認,由此說明該數(shù)值計算方法可行和可靠。
對于潛艇的上浮運動問題,采用重疊網(wǎng)格,在一個時間步長內(nèi),依據(jù)所需的工況,輸入初始條件和邊界條件,初始化流場后可運行計算,基于控制方程求解艇體表面的速度場和壓力場,進而可求出艇體表面所受的力和力矩,將所求得的力和力矩帶入潛艇6DOF運動方程中就可求取所需艇體上浮運動參數(shù)。潛艇在靜水中上浮運動的數(shù)值計算方法如圖9所示。
圖9 數(shù)值計算的流程圖Fig.9 Flowchart of numerical calculation
試驗工況為:潛艇質(zhì)量為640.48 kg,正上浮力為621.56 N,潛深3.5 m,靜水,航速為零。在此工況下開展?jié)撏细∵\動的數(shù)值計算。在進行數(shù)值模擬時,選取分離式(segregated flow)求解控制方程,對于壓力和速度的修正選取壓力耦合方程的半隱式算法(semi-implicit-method for pressure linked equations,SIMPLE),網(wǎng)格離散格式選用二階迎風(fēng)格式。涉及非穩(wěn)態(tài)計算時,設(shè)置數(shù)值計算過程中的時間步長Δt=5×10-3s。對于湍流模型,選取Standardk-ω。試驗場地選擇室外波浪水池。全附體SUBOFF模型上浮試驗如圖10所示,數(shù)值計算結(jié)果與試驗結(jié)果的對比如圖11所示。
圖10 SUBOFF模型上浮試驗Fig.10 Floating test of SUBOFF model
(a) 縱向速度(a) Longitudinal speed
從圖11中可知,數(shù)值計算結(jié)果與上浮試驗結(jié)果的變化規(guī)律大致相同,數(shù)值計算結(jié)果與試驗結(jié)果具有較好的吻合度;數(shù)值計算的縱向速度和垂向速度與試驗結(jié)果相比很接近,縱傾角和橫傾角的數(shù)值計算結(jié)果與試驗結(jié)果存在有一定的誤差,但誤差較小。由此可知,該數(shù)值計算具有較好的可信度,從而驗證了該數(shù)值計算方法可以適用于潛艇上浮運動的研究。
本文介紹潛艇上浮運動的數(shù)值計算方法,并對數(shù)值計算中湍流模型、離散格式等進行了優(yōu)選,介紹了潛艇的模型上浮試驗,對數(shù)值計算的不確定度進行了分析,且與上浮模型試驗的結(jié)果進行了對比,主要得到以下結(jié)論:
1)通過分別對不同湍流模型以及不同離散格式的對比分析,優(yōu)選了Standardk-ω的湍流模型,選取了二階迎風(fēng)格式,由此確定了數(shù)值計算的基礎(chǔ);
2)基于網(wǎng)格收斂性和時間步長收斂性的分析,對數(shù)值計算的不確定度開展了研究,最終證明了數(shù)值模擬的可靠性;
3)針對潛艇六自由度上浮運動問題,采用重疊網(wǎng)格技術(shù),并結(jié)合6-DOF剛體運動模型和DFBI模型,對潛艇上浮運動數(shù)值計算流程進行了闡述,同時對潛艇的模型上浮試驗進行了介紹,該數(shù)值模擬方法和上浮試驗可支撐潛艇上浮運動的相關(guān)研究;
4)基于模型上浮試驗,數(shù)值計算結(jié)果與試驗結(jié)果具有較好的吻合度,且誤差較小,從而證明了該數(shù)值模擬方法的可行性和準(zhǔn)確性。