梁景奇,徐保成,王 瑞,祁曉斌,李瑞杰

(西北機(jī)電工程研究所, 陜西咸陽(yáng) 712099)

0 引言

超空泡航行器大部分表面被空泡包覆,航行阻力顯著減小,航行速度大幅提高。超空泡射彈就是一種利用了超空泡減阻技術(shù)的動(dòng)能武器,具備水下高亞聲速甚至超聲速航行的能力,可以對(duì)魚雷、水雷等水下目標(biāo)實(shí)現(xiàn)快速毀傷[1]。

超空泡射彈的重要特征之一是射彈的流體動(dòng)力受到超空泡與彈體之間相對(duì)位置耦合關(guān)系的影響。當(dāng)航速在300～1 200 m/s時(shí),受到初始擾動(dòng)后,包裹于超空泡內(nèi)的射彈開始周期性地與空泡壁面發(fā)生碰撞,即尾拍運(yùn)動(dòng)。Kulkarni[2]等人基于Milwitzky[3]的理論采用加入附加質(zhì)量的方法研究超空泡運(yùn)動(dòng)體的“尾拍”力;Euteneuer[4]建立了超空泡運(yùn)動(dòng)體的動(dòng)力學(xué)模型,將后體滑行力作為外界高頻擾動(dòng)看待,忽略了空泡的時(shí)間延遲效應(yīng);王茂勵(lì)[5]基于彈道仿真技術(shù)研究了流體動(dòng)力參數(shù)的攝動(dòng)對(duì)彈道穩(wěn)定性的影響規(guī)律;魏英杰[6]等提出了超空泡射彈流體動(dòng)力特性的動(dòng)態(tài)描述方法,建立了射彈尾拍運(yùn)動(dòng)的無(wú)約束動(dòng)力學(xué)方程;張廣[7]等運(yùn)用CFX建立了超空泡航行體三維6DOF數(shù)值模型,發(fā)現(xiàn)無(wú)控高速超空泡航行體具有運(yùn)動(dòng)穩(wěn)定性;時(shí)素果[8]等開展了自由飛試驗(yàn),獲得了不同預(yù)置舵角下超空泡航行體水平運(yùn)動(dòng)的彈道特性。水下高速超空泡射彈的彈道穩(wěn)定性與其尾拍運(yùn)動(dòng)具有重要關(guān)系,彈道失穩(wěn)會(huì)造成射擊精度和射程不足,影響射彈的作戰(zhàn)效能。因此,對(duì)射彈空泡流場(chǎng)與彈體耦合的尾拍特性的進(jìn)一步研究極具必要性。

依托CFD軟件FLUENT18.0及其二次開發(fā),結(jié)合動(dòng)網(wǎng)格移動(dòng)計(jì)算域技術(shù),忽略重力因素,在考慮液體壓縮性的基礎(chǔ)上建立了高速射彈剛體運(yùn)動(dòng)與空泡流場(chǎng)耦合數(shù)值模型,從理論及實(shí)驗(yàn)方面驗(yàn)證了數(shù)值模型的準(zhǔn)確性;求解超空泡射彈水下縱平面運(yùn)動(dòng),探究水下運(yùn)動(dòng)穩(wěn)定規(guī)律,分析了初始速度對(duì)射彈尾拍特性的影響。

1 數(shù)值模型

1.1 控制方程

采用VOF多相流模型模擬相界面運(yùn)動(dòng)。VOF是一種在固定Euler網(wǎng)格下的界面捕捉法,常用于由兩種及以上不相混液體組成的流體中,適用于多相間有清晰界面的流動(dòng)。

1.1.1 連續(xù)性方程

(1)

(2)

(3)

式中:ρ是流體混合密度;u是混合速度;n1是相數(shù),由于考慮了不凝氣體,在這里,n1=3;αk、ρk和μk分別為第k相的體積分?jǐn)?shù)、密度和速度。

1.1.2 動(dòng)量方程

勻質(zhì)平衡流模型認(rèn)為各相間不存在速度差,在流體微元中基于平均密度和平均動(dòng)力粘度的動(dòng)量守恒方程可描述為:

(4)

式中：ρm、μm分別為流體微元的各相平均密度和動(dòng)力粘度,且有:

(5)

(6)

1.2 空化模型

采用Schnerr and Sauer空化模型模擬超空泡射彈的水下穩(wěn)定空化繞流,Schnerr and Sauer空化模型將汽相體積分?jǐn)?shù)和單位體積流體含有的空泡數(shù)量聯(lián)系起來(lái),表達(dá)式為:

(7)

(8)

式中:ρl為液體密度;ρv為水蒸氣密度;ρ為混合相密度；pv為飽和為蒸氣壓;p為當(dāng)前壓力;αnuc為氣核體積分?jǐn)?shù);RB為氣核空泡直徑;n2為單位體積內(nèi)的空泡數(shù)量。

1.3 湍流模型

Realizabledk-ε湍流模型主要針對(duì)充分發(fā)展的湍流,穩(wěn)定性好,近壁面區(qū)域的流動(dòng)狀態(tài)使用壁面函數(shù)預(yù)測(cè),對(duì)邊界層網(wǎng)格要求較為寬松;結(jié)合尺度化壁面函數(shù),其在不增加計(jì)算量和保持模型穩(wěn)定性的前提下增加了模型的適用范圍,模擬效果好。

湍流強(qiáng)度k方程:

(9)

湍流耗散率ε的方程:

(10)

式中:μt為湍動(dòng)粘度;μ是流體的時(shí)均速度;σk、σε分別為k、ε方程的湍流能量普朗特?cái)?shù);C1、C2為經(jīng)驗(yàn)常數(shù);E為時(shí)均應(yīng)變率;ν為運(yùn)動(dòng)粘度;xi、xj為各方向距離。

1.4 可壓縮液體模型

Tait方程是通過(guò)采用非線性回歸的方法,對(duì)能夠反應(yīng)p-v-T三者關(guān)系的試驗(yàn)數(shù)據(jù)進(jìn)行擬合,而得到的液體狀態(tài)方程,廣泛應(yīng)用于描述可壓縮液體的物性。沒(méi)有溫度修正的簡(jiǎn)化Tait液體狀態(tài)方程可描述為:

(11)

K=K0+n3ΔP

(12)

Δp=p-p0

(13)

(14)

式中:p0為參考?jí)毫?ρ0為參考?jí)毫ο碌囊后w密度;K0為參考?jí)毫ο碌囊后w體積彈性模型;n3為密度指數(shù);p為當(dāng)前壓力;ρ為當(dāng)前壓力下的液體密度;K為當(dāng)前壓力下的液體體積彈性模量;c為水中聲速。

1.5 動(dòng)網(wǎng)格耦合運(yùn)動(dòng)

耦合運(yùn)動(dòng)即射彈運(yùn)動(dòng)與空泡流場(chǎng)的計(jì)算相互耦合。已知前一時(shí)刻剛體重心位置和偏轉(zhuǎn)角,通過(guò)對(duì)物體表面壓力和剪切應(yīng)力積分得到流體動(dòng)力和力矩,再根據(jù)剛體運(yùn)動(dòng)方程計(jì)算物體運(yùn)動(dòng)的平移速度和角速度,然后重新計(jì)算重心位置和偏轉(zhuǎn)角。

(15)

(16)

式中:G為變換矩陣,剛體的位置與方向根據(jù)每時(shí)間步線速度和角速度的變化而變化,c.g.表示質(zhì)心。

采用移動(dòng)計(jì)算域技術(shù)模擬彈體運(yùn)動(dòng),計(jì)算過(guò)程中僅涉及計(jì)算域的移動(dòng),不存在網(wǎng)格的變形與重構(gòu),計(jì)算效率高,結(jié)果一致性好。

2 計(jì)算模型及方法驗(yàn)證

2.1 模型、網(wǎng)格劃分及邊界條件

文中研究的射彈模型如圖1所示,射彈采用圓盤空化器,前部為兩段錐段,中間為圓柱段,柱段尾部安裝有6片尾翼。前端空化器直徑3.2 mm,圓柱段最大直徑15 mm,質(zhì)量為0.23 kg。

圖1 射彈模型圖

采用圓柱形計(jì)算域,直徑取50倍尾截面空泡直徑,計(jì)算域軸向長(zhǎng)度為11倍彈長(zhǎng),前端邊界距離空化器4倍彈長(zhǎng),后端邊界距離彈尾6倍彈長(zhǎng),該計(jì)算域徑向尺度可以忽略空泡阻塞效應(yīng)[9]。

針對(duì)所建立的三維計(jì)算域,采用ICEM軟件的O-Block技術(shù)劃分全結(jié)構(gòu)化網(wǎng)格,如圖2所示在彈體周圍3 mm范圍內(nèi)的流域劃分外O-block用于設(shè)置邊界層網(wǎng)格,近壁面添加邊界層網(wǎng)格,并根據(jù)y+值對(duì)網(wǎng)格進(jìn)行優(yōu)化。劃分網(wǎng)格時(shí)特別注意在空泡兩相交界面位置進(jìn)行網(wǎng)格加密,最終劃分的網(wǎng)格總數(shù)約80萬(wàn),網(wǎng)格質(zhì)量均在0.6以上。

圖2 剖面網(wǎng)格劃分圖

計(jì)算域四周邊界均設(shè)置為壓力入口條件(pressure-inlet),設(shè)置靜壓值且靜壓值隨深度變化,射彈表面的邊界條件設(shè)置為壁面(wall),并且壁面與臨界網(wǎng)格相對(duì)靜止。計(jì)算域和邊界條件設(shè)置如圖3所示。

圖3 計(jì)算域及邊界條件圖

2.2 數(shù)值模型驗(yàn)證

在超空泡射彈運(yùn)動(dòng)中,運(yùn)動(dòng)特性取決于流體動(dòng)力特性,流體動(dòng)力取決于彈體與空泡的相對(duì)位置關(guān)系,因此空泡形態(tài)計(jì)算準(zhǔn)確性是運(yùn)動(dòng)模擬準(zhǔn)確性的前提。將文獻(xiàn)[10]經(jīng)驗(yàn)公式和文獻(xiàn)[11]試驗(yàn)結(jié)果與數(shù)值模型的空泡外形作對(duì)比,驗(yàn)證數(shù)值模型準(zhǔn)確性。

2.2.1 經(jīng)驗(yàn)公式驗(yàn)證

Savechenko[10]等在烏克蘭國(guó)家科學(xué)院水動(dòng)力實(shí)驗(yàn)室進(jìn)行了空泡截面獨(dú)立擴(kuò)張?jiān)碛行则?yàn)證試驗(yàn),給出了空泡外形描述經(jīng)驗(yàn)公式:

(17)

(18)

(19)

式中:Rn為空化器半徑;Rc,max、Lc,max分別為超空泡最大截面半徑和空泡全長(zhǎng);x為空泡截面到空化器的距離;σ為空化數(shù),σ=0.012～1.057。

對(duì)比同工況下經(jīng)驗(yàn)公式和數(shù)值模型計(jì)算的空泡輪廓,驗(yàn)證數(shù)值模型的準(zhǔn)確性。取水深1 m,航速100 m/s,空化器半徑1.6 mm,前段空泡輪廓對(duì)比如圖4所示。

圖4 空泡輪廓對(duì)比圖

圖4對(duì)比結(jié)果顯示,同工況下經(jīng)驗(yàn)公式得到的空泡輪廓在徑向尺度略大于CFD計(jì)算的空泡輪廓,計(jì)算誤差小于12%,可以證明數(shù)值模型的準(zhǔn)確性。

2.2.2 試驗(yàn)結(jié)果驗(yàn)證

Hrubes[11]采用試驗(yàn)方法,在水深4 m,速度970 m/s條件下水平發(fā)射平頭錐形射彈,利用高速攝像機(jī)捕捉射彈入水后的空泡形態(tài)。試驗(yàn)結(jié)果如圖5所示,射彈完全被空泡包裹。運(yùn)用建立的數(shù)值模型,采用與文獻(xiàn)[11]相同的彈體模型及初始條件,模擬射彈水下直航運(yùn)動(dòng),空泡輪廓對(duì)比如圖6所示。

圖5 970 m/s初速試驗(yàn)與數(shù)值模擬空泡形態(tài)圖

圖6 數(shù)值模擬與試驗(yàn)空泡輪廓對(duì)比圖

圖6給出了Hrubes試驗(yàn)數(shù)據(jù)點(diǎn)和數(shù)值模擬的空泡外形曲線。從中可以看出,數(shù)值模擬結(jié)果與試驗(yàn)數(shù)據(jù)點(diǎn)基本吻合,均能有效包裹射彈,最大誤差不超過(guò)10%,表明文中數(shù)值模型合理可行。

3 計(jì)算結(jié)果及分析

采用變參數(shù)法探究不同初速下超空泡射彈尾拍的流體動(dòng)力特性及運(yùn)動(dòng)特性,水深恒定1 m,忽略重力、海流和波浪影響,考慮液體壓縮性,模擬定深繞質(zhì)心旋轉(zhuǎn)的縱平面尾拍運(yùn)動(dòng)[12]。先定常計(jì)算初始空泡流場(chǎng),再基于三維動(dòng)網(wǎng)格模型并結(jié)合二次開發(fā)技術(shù)模擬自由減速運(yùn)動(dòng),探究運(yùn)動(dòng)穩(wěn)定性及流體動(dòng)力、彈道參數(shù)的變化規(guī)律,初始條件如表1所示。

表1 初始條件表

3.1 流體動(dòng)力特性

探究1 m水深、30 rad/s角速度擾動(dòng)條件下不同初速射彈的流體動(dòng)力變化特性,以計(jì)算時(shí)長(zhǎng)0.015 s為計(jì)算終止條件,時(shí)間步長(zhǎng)根據(jù)計(jì)算殘差在1e-08到1e-05之間逐漸增加,計(jì)算結(jié)果如圖7、圖8所示。

圖7 尾拍升力對(duì)比圖

圖8 尾拍阻力對(duì)比圖

圖7、圖8分別為3種速度工況尾拍升、阻力對(duì)比曲線。初始擾動(dòng)使彈體軸線轉(zhuǎn)動(dòng),并與空泡壁面發(fā)生周期性碰撞,彈體的尾拍升、阻力均呈現(xiàn)周期性振蕩特性。根據(jù)曲線特征,可將受擾動(dòng)后超空泡射彈的振蕩過(guò)程分為發(fā)散階段和收斂階段。在第一階段,由于彈體穿刺空泡深度不斷增加,尾拍升、阻力曲線均呈發(fā)散特性,振蕩幅值迅速增大并達(dá)到最大,振蕩周期不斷減小;在第二階段,由于速度衰減,升、阻力曲線呈收斂特性,振蕩幅值逐漸減小,振蕩周期緩慢增大,兩個(gè)階段升、阻力曲線的幅值和周期均同步變化,同一時(shí)刻升力幅值約為阻力幅值的3.5倍。

尾拍力與速度的平方正相關(guān),初速越大,振蕩周期越小,尾拍升、阻力越快達(dá)到最大幅值,最大幅值越大,均值越大;隨著速度的衰減,3種工況下阻力和升力均值的差距逐漸減小;0.01 s后時(shí)間步長(zhǎng)增至1e-05,振蕩曲線呈階梯形,但依然能客觀反映尾拍力振蕩特性。

3.2 空泡形態(tài)

以1 200 m/s工況為例,分別取振蕩發(fā)散階段和振蕩收斂階段的一個(gè)尾拍周期T,分析尾拍過(guò)程中空泡形態(tài)變化及彈體與空泡的位置關(guān)系,結(jié)果如圖9、圖10所示。

圖9 振蕩發(fā)散階段尾拍動(dòng)態(tài)圖(第一周期)

圖9為振蕩發(fā)散階段(第一個(gè)尾拍周期)內(nèi)的空泡動(dòng)態(tài)變化過(guò)程。彈體向下偏轉(zhuǎn),尾翼局部穿刺空泡下壁面后產(chǎn)生二次空泡,升、阻力達(dá)到最大,彈體回彈并觸碰上壁面,最終回到起始位置。在該階段,尾拍角幅值隨穿刺深度的增加呈增大趨勢(shì),尾拍角為0°時(shí),彈體除空化器外不沾濕,空泡軸線與彈體軸線基本重合。

圖10為振蕩收斂階段(第6個(gè)尾拍周期)的空泡動(dòng)態(tài)變化過(guò)程,在該階段,彈體穿刺空泡深度更大,空泡形態(tài)發(fā)生顯著改變。到達(dá)最大尾拍角時(shí),除尾翼外,彈體第二錐段和靠近尾翼的圓柱段也發(fā)生沾濕,且兩次沾濕區(qū)域具有不對(duì)稱性,而尾拍角為0°時(shí)彈體依然存在沾濕情況,這是由于空泡變化與彈體振蕩周期存在相位差,空泡軸線變化滯后于射彈軸線變化,即空泡時(shí)間延遲效應(yīng)。

3.3 運(yùn)動(dòng)特性

以地面系為參考坐標(biāo)系,分析3種初速射彈受擾動(dòng)后的尾拍運(yùn)動(dòng)特性,俯仰角、俯仰角速度、軸向速度、質(zhì)心位移對(duì)比曲線分別如圖11～圖14所示。

圖11、圖12為俯仰角和俯仰角速度隨時(shí)間變化曲線。不同速度下的俯仰角、俯仰角速度曲線均呈周期性振蕩特性,同樣可分為振蕩發(fā)散和振蕩收斂階段,變化規(guī)律與流體動(dòng)力曲線相似。由俯仰角曲線可知,初速越大,振蕩周期越小,振蕩發(fā)散階段俯仰角幅值越小;在振蕩收斂階段,3種工況俯仰角幅值均接近4°,基本保持恒定,這表明,振蕩收斂階段彈體俯仰角幅值受初速影響極小。

圖11 俯仰角隨時(shí)間變化曲線

圖12 俯仰角速度隨時(shí)間變化曲線

由俯仰角速度曲線可知,初速越大,俯仰角速度曲線振蕩周期越小,越快達(dá)到最大幅值,最大幅值越大,1 200 m/s工況最大俯仰角約為270 rad/s,增幅高達(dá)9倍;達(dá)到最大幅值后,俯仰角速度振蕩幅值均隨時(shí)間逐步衰減。通過(guò)與圖7、圖8對(duì)比可知,尾拍力曲線俯仰角曲線存在相位差,亦是由空泡延遲效應(yīng)所致。

圖13為X向速度隨時(shí)間變化曲線,3條曲線均呈單調(diào)遞減趨勢(shì),初速度越大,阻力越大,因此速度衰減越快,0.014 s內(nèi)1 200 m/s、1 000 m/s、800 m/s工況的速度衰減幅度分別為260 m/s、190 m/s和130 m/s;隨著時(shí)間的增加3種工況航速差距逐漸縮小,由此可知通過(guò)增加初速增加射程效果有限。

圖13 X向速度隨時(shí)間變化曲線

圖14為0.015 s內(nèi)3種工況射彈質(zhì)心位移變化曲線,質(zhì)心隨著彈體的振蕩均產(chǎn)生了側(cè)向位移,且位移量不斷增加。以1 200 m/s為例,軸向位移15 m時(shí),側(cè)向位移為0.02 m,與軸向位移量相比,側(cè)向位移量可忽略不計(jì),彈體基本保持了直線航行,這表明超空泡射彈的水下彈道具有動(dòng)態(tài)穩(wěn)定性。

圖14 質(zhì)心位移變化曲線

4 結(jié)論

文中基于建立的流場(chǎng)運(yùn)動(dòng)耦合數(shù)值模型,數(shù)值仿真了超空泡射彈的尾拍運(yùn)動(dòng),分析了初速對(duì)尾拍力的特性及運(yùn)動(dòng)特性的影響,結(jié)論如下:

1)受擾動(dòng)后超空泡射彈的水下運(yùn)動(dòng)具有動(dòng)態(tài)穩(wěn)定特性,尾拍力、俯仰角、俯仰角速度等參數(shù)均呈現(xiàn)準(zhǔn)周期振蕩特性,可將該過(guò)程分為振蕩發(fā)散階段和振蕩收斂階段；

2)尾拍運(yùn)動(dòng)使空泡軸線產(chǎn)生偏移,空泡軸線變化滯后于射彈軸線變化,同一尾拍周期的兩次沾濕區(qū)域具有不對(duì)稱性；

3)初速度越大,尾拍周期越小,彈體尾拍升、阻力及俯仰角速度振蕩幅值越大,俯仰角振蕩幅值受初速影響極小。