呂紅慶,許 磊,王振清
(1.煙臺哈爾濱工程大學研究院, 山東 煙臺 264000; 2.哈爾濱工程大學 航天與建筑工程學院, 哈爾濱 150001)
入水過程是跨介質武器從空中飛行轉入水下航行的一個重要環(huán)節(jié),旋成體以一定初速度撞擊自由液面,將自身一部分動能傳遞給液體,液體發(fā)生復雜流動并與航行體表面分離。隨著入水深度的增加,旋成體兩側形成了穩(wěn)定的空泡區(qū)域,空泡區(qū)域主要由液體相變產(chǎn)生的水蒸氣和高速誘導進入的空氣組成。此外,在入水過程中,旋成體頭部產(chǎn)生具有強瞬時性、非定常性及高載荷性的沖擊載荷,同時伴隨著相變、湍動等非定常流體動力學現(xiàn)象以及彈體彈道難可控性等問題[1]。因此研究旋成體不同頭型在入水過程的空泡形態(tài)演化與瞬間沖擊載荷特性對跨介質兵器的發(fā)展具有指導意義。
早在十九世紀末期,學者們就已針對入水問題開展了試驗和理論研究,Worthington[2]利用剛興起的閃光攝影技術對小球入水進行了試驗研究,通過改變液體介質,觀察到入水過程中液體飛濺與空泡演化現(xiàn)象,同時分析了液體屬性、表面張力等因素對液體飛濺現(xiàn)象的影響,為之后的入水試驗奠定了基礎。Watanabe[3]通過對尖錐體進行入水試驗研究,測量出尖錐體在不同入水時刻所受到?jīng)_擊載荷大小。在早期獲得的試驗數(shù)據(jù)基礎上,最早進行入水問題理論研究的是T.von Karman[4],他假定流體無黏、無旋和不可壓縮,物體在穿越氣-水界面時液面不發(fā)生變化,同時引入速度勢理論,設定自由液面速度勢為零,提出附加質量法,基于動量守恒定律推導出入水沖擊載荷計算公式,該計算公式對于二維小底升角(α)楔形體入水沖擊載荷預測較為準確。Wagner[5]提出漸進匹配方法,針對物體入水過程中液面抬升現(xiàn)象,將流體分為內外流域,內流域中進行局部射流分析,外流域提出近似平板理論,使物體入水過程的數(shù)學描述更加符合實際問題。Wagner的理論研究思想為后續(xù)入水沖擊理論奠定了堅實的基礎。
隨著入水問題數(shù)學模型的完善,Logvinovich[6]在他的專著中提出了空泡獨立膨脹原理,根據(jù)相似理論給出了計算空泡尺寸的半經(jīng)驗公式,成為空泡理論的奠基人。Lee等[7]在空泡理論基礎之上,推導了普適性的數(shù)學模型,能夠解釋部分簡單的入水現(xiàn)象。與此同時,隨著高性能計算機和先進計算方法的出現(xiàn),對入水問題的研究迎來了黃金階段。Park等[8]提出了計算高速入水沖擊載荷以及跳彈現(xiàn)象的數(shù)值方法,假定入水沖擊在極短的時間內發(fā)生,入水區(qū)域的流動為無粘勢流,采用源面板法進行求解,計算獲得的數(shù)值數(shù)據(jù)與試驗數(shù)據(jù)具有良好的一致性。Neaves等[9]采用有限體積法,引入自然空化模型和Tait狀態(tài)方程模擬了射彈高速垂直入水空泡演化規(guī)律,同時考慮了流體介質的壓縮性和空化潛熱效應,獲得了高速射彈流場結構,數(shù)值結果與實驗結果一致。Abraham等[10]對小球在不同雷諾數(shù)、入水速度和液體表面張力情況下展開了數(shù)值模擬,認為動量傳遞的原因主要是拖曳力而非摩擦阻力。
蔣運華等[11]研究了彈體在不同通氣量和不同傾斜角度入水情況下近水表面空泡特性,包括表面封閉、頸縮、深閉合以及壁面波動等現(xiàn)象,同時針對超音速彈體串并聯(lián)入水過程中空泡幾何形狀和阻力特性展開了研究[12]。王聰?shù)萚13]對不同角度錐頭彈體高速垂直入水進行了數(shù)值模擬,得到不同頭型條件下高速入水運動參數(shù)及空泡形態(tài)發(fā)展規(guī)律、流場的壓力分布及速度分布規(guī)律,分析了頭型對入水空泡流場的影響。段文洋等[14]采用試驗與數(shù)值模擬相結合的辦法,研究了彈體入水初期的空泡流動與演化過程,分析了頭部形狀對該過程的影響。潘光等[15]提出一種非對稱機頭的AUV,采用數(shù)值模擬的方法分析了在不同的入水速度、機頭的形狀對AUV入水過程的彈道特性的影響。
目前大多數(shù)的入水問題研究主要集中于旋成體低速入水過程,對于旋成體由低速到高速入水的研究還較少。因此本文對跨介質旋成體低速到高速入水過程開展數(shù)值模擬研究,基于雷諾平均N-S方程、VOF多相流模型、realizablek-ε湍流模型,引入Schnerr and Sauer空化模型,建立能較好描述入水過程的數(shù)值計算模型,采用重疊網(wǎng)格技術進行仿真計算。對不同頭型旋成體在不同工況下的入水過程展開研究,分析旋成體在入水過程中流場參數(shù)、空泡演化、彈體動力學特性與瞬間沖擊載荷特性。所得研究成果可為跨介質兵器設計提供一定參考。
物理模型采用旋成體構型,以空投魚雷、跨介質射彈等幾何特征的跨介質武器作為目標模型。彈體長度設定為航行體直徑的3倍且不考慮尾翼的影響,計算所采用的具體模型如圖1所示。
圖1 旋成體物理模型示意圖Fig.1 Physical model of the vehicle
旋成體的直徑為324 mm,彈身長度為972 mm,球頭半徑為162 mm,航行體材料為普通碳素鋼,密度為7.85 g/cm3。具體幾何尺寸如表1所示。由于針對旋成體入水過程流場參數(shù)變化與空泡演化特性開展研究,所以將模型設定為剛體結構不考慮旋成體變形與破壞。
表1 旋成體參數(shù)
本文采用均質平衡多相流理論,選取Volume of Fluid(VOF)模型[16],模擬兩相或多相具有明顯分界面的流體運動問題。將氣相與液相作為混合的單一介質,通過各組分的所占混合物的體積分數(shù)比值來描述各組分的流動狀態(tài),不同相共用同一套基本控制方程組,根據(jù)質量守恒方程,連續(xù)性方程可表示為:
(1)
(2)
(3)
(4)
式中:P為壓力;S為源項;τij為剪切應力,表達式如下:
(5)
湍流模型選取為realizablek-ε模型[17],該模型引入了旋轉和曲率相關的數(shù)值,采用渦旋粘性各向同性假設,在模擬流動分離、旋轉流動、圓孔射流和射流撞擊等問題中具有良好適應性[18]。特別是對于射流曲率較大的情況,realizablek-ε模型具有很好的表現(xiàn)。該模型湍動能k和湍動能的耗散率ε的輸運方程表達為:
(6)
式(6)中,我們通常取C1ε=1.44、C2=1.9、σk=1、σε=1.2。
realizablek-ε模型是對standardk-ε模型的修正,也被稱之為可實現(xiàn)的k-ε模型。在計算過程中,為了保證計算結果的可實現(xiàn)性(Realizability),Cμ不再視為常數(shù),來保證正應力數(shù)學約束關系的實現(xiàn)。
預測常數(shù)Cμ的公式如下:
(7)
本文為確保仿真更加符合物理實際,所以采用Schnerr and Sauer模型[19]對彈體入水過程中的空化現(xiàn)象進行求解。水蒸氣相的質量輸運方程如下式所示:
(8)
式中:RB為氣穴半徑,RB=1×10-6m。不可凝結氣體體積分數(shù)αnuc=5×10-4、預測常數(shù)Fvap、Fcond的取值分別為50與0.001。
本文在計算過程中采用有限體積法對構建的控制方程進行離散化處理,利用PISO算法,對壓力和速度耦合問題求解,該算法對SIMPLE算法進行了改進,是一種基于壓力的隱式算子分裂方法,計算收斂性健壯且計算效率高。控制方程中對流項采用二階離散格式,擴散項采用二階中心差分格式,混合相體積率的離散采用CICSAM格式,時域采用隱式離散方法進行求解。
網(wǎng)格劃分采用重疊網(wǎng)格技術[20],該技術能夠確保在彈體入水過程中網(wǎng)格不發(fā)生變形,重疊區(qū)域網(wǎng)格尺寸與背景網(wǎng)格處尺寸基本保持一致,這樣物體在運動過程中運動網(wǎng)格與背景網(wǎng)格插值區(qū)域的網(wǎng)格質量保持最優(yōu)。本研究的計算域由彈體運動域和背景區(qū)域兩部分構成,如圖2所示。
圖2 計算域示意圖Fig.2 Schematic diagram of calculation domain
圖2(a)所示,運動域為圍繞彈體的鈍頭體圓柱,圓柱直徑為彈體直徑的2倍,計算域長度為1 200 mm,約為4倍彈體直徑。圖2(b)所示,為計算背景區(qū)域,該區(qū)域采用立方體結構,長為6 000 mm,寬為3 000 mm,高為6 000 mm,上半部分為空氣域,高度為1 500 mm,下半部分為水域,水深設定為4 500 mm,彈體的初始位置設定在空氣域中,距離水面高度為300 mm。由于本次計算的模型的為軸對稱模型,所以采用一半的模型進行計算以節(jié)約計算的資源和時間。
采用結構網(wǎng)格劃分運動區(qū)域與背景域,如圖3所示。計算域y-z平面的網(wǎng)格的剖分視圖,z軸的負方向為重力加速度的方向,坐標的原點為旋成體模型質心位置。在旋成體周圍進行局部的加密處理,以更好捕捉入水過程中空泡的演化與瞬間的沖擊載荷。臨近壁面處,法向速度會存在較大的梯度,對于該區(qū)域的計算采用壁面函數(shù)法。本次計算過程中雷諾數(shù)較高且存在分離與射流流動的現(xiàn)象,我們使用非平衡壁面函數(shù)法,采用結構化網(wǎng)格進行劃分,運動域網(wǎng)格數(shù)量為60萬,背景域網(wǎng)格數(shù)量為240萬。背景計算域的四周邊界條件與底面設定為壁面(wall),它的上表面設定為壓力出口(pressure-outlet),在旋成體周圍的區(qū)域邊界設定為重疊網(wǎng)格可以識別的over-set邊界條件,下落過程中,采用六自由度運動方式,背景域網(wǎng)格與運動域網(wǎng)格進行插值計算傳遞數(shù)據(jù)。
圖3 y-z平面剖分網(wǎng)格示意圖Fig.3 y-z plane mesh generation
有必要驗證本文中所采用的數(shù)值計算方法的可靠性,以確保研究結果的可信性。南京理工大學瞬態(tài)物理國家重點實驗室[21]關于不同頭型回轉體入水的試驗數(shù)據(jù)為參考,將數(shù)值模擬得到的空泡演化和回轉體入水速度變化與試驗結果進行對比。
驗證計算的試驗裝置模型如圖4所示,試驗中回轉體為長度40 mm、直徑8 mm的實心圓柱體,材料為普通碳素鋼,密度為7.85 g/cm3。液體介質為室溫下的水,密度為1 g/cm3,水溫為25 ℃,通過高速攝像機記錄回轉體在入水過程中對自由液面的擾動。
圖4 試驗裝置模型示意圖Fig.4 Schematic diagram of the test equipment
試驗在高度為500 mm、底部尺寸為250 mm×250 mm且設有防護層的水槽中進行。在水槽的側面貼有坐標紙,記錄回轉體運動的位置變化?;剞D體的初始位置在水槽的上方,被電磁鐵吸附固定,釋放后自由落體垂直入水,接觸水面的入水速度為1.67 m/s。通過圖5分析,可以清楚地看到,入水仿真所獲得的回轉體被空泡包裹、空泡的形成與演化、空泡的頸縮與閉合以及產(chǎn)生液體的射流等現(xiàn)象與試驗結果基本一致,驗證了本文數(shù)值計算方法的適用性。
圖5 數(shù)值仿真與試驗現(xiàn)象過程示意圖Fig.5 Comparison between numerical simulation and experimental phenomena
在數(shù)值求解的過程中,離散點分布在網(wǎng)格的節(jié)點上,對于物理量過渡平緩的位置,網(wǎng)格節(jié)點布置可以較為稀疏,對于梯度變化較大區(qū)域,需要增加網(wǎng)格節(jié)點數(shù)量。選擇適當?shù)木W(wǎng)格疏密程度不僅可以保證仿真結果的精確度,還能夠節(jié)約數(shù)值計算資源與時間。使用2.1~2.3節(jié)中給定的計算模型與邊界條件,以半球頭旋成體100 m/s的初速進行垂直入水計算,選取出三套不同密度的網(wǎng)格模型,網(wǎng)格數(shù)量見表2。通過驗證計算,最終確定出適合本次計算的網(wǎng)格尺度。
表2 網(wǎng)格劃分數(shù)量
計算過程中設置監(jiān)視點,得到旋成體入水過程中頭部最大壓力和速度隨時間變化的曲線,以一個標準大氣壓P0=101 325 Pa作為沖擊壓力峰值的基準,對頭部壓力峰值進行無量綱化的處理。如圖6與圖7所示??梢钥吹剑谟嬎氵^程中三套網(wǎng)格壓力與速度隨時間的變化趨勢基本一致,在入水初期,三套網(wǎng)格對于壓力峰值與后期壓力變化的計算曲線基本一致。觀察旋成體入水速度變化,可以看到稀疏網(wǎng)格的速度曲線雖然趨勢與另外2種網(wǎng)格密度保持一致,但在不同時刻時與其他2種情況的航行體速度存在較大差異。通過對計算的數(shù)據(jù)差異以及所需的計算資源進行比較,最終確定在后續(xù)模擬中選用運動域60萬網(wǎng)格、背景域240萬網(wǎng)格作為本次計算的網(wǎng)格數(shù)量。
圖6 旋成體無量綱壓力峰值隨時間變化曲線Fig.6 Dimensionless pressure peaks of axisymmetric body varying with time
圖7 旋成體速度隨時間變化曲線Fig.7 Variation of speed of axisymmetric body with time
影響旋成體的入水流場特性的因素有很多,本小節(jié)主要研究相同長徑比、相同質量密度情況下,在100 m/s垂直入水速度下,頭型對入水過程的影響。圖8給出了不同頭型在相同入水深度時空泡形態(tài),橫坐標為空泡半徑(R)與旋成體半徑(R0)的比值,縱坐標為空泡長度(H)與旋成體直徑(D)的比值,我們觀察可知3種頭型旋成體空泡均在旋成體的肩部產(chǎn)生,空泡的整體發(fā)展輪廓基本一致。
空泡沿著徑向逐漸向外擴張,空泡的半徑隨之增大,自由液面處的空泡半徑最大。不同頭型旋成體入水初期,在相同深度的情況下,平頭體產(chǎn)生的空泡半徑最大,半球頭體次之,尖錐頭體的空泡尺寸最小。對于不同頭型入水空泡大小差異這一現(xiàn)象,我們可以根據(jù)文獻[7],假設流體介質為無黏、無旋、不可壓縮的理想流體,根據(jù)能量守恒定律,彈體在入水過程中,自身動能的損失都轉化為液體壓力勢能和動能,由此可以得到入水空泡半徑預測公式,經(jīng)過一系列推導,并略去高階小量,可知空泡半徑的大小主要與彈體頭型在空化數(shù)為零時的阻力系數(shù)及入水時間有關,阻力系數(shù)越大那么在自由液面處形成的空泡尺寸就越大,在未發(fā)生空泡閉合時,旋成體入水時間越長,自由液面處空泡半徑越大。
圖8 不同入水深度下空泡形態(tài)曲線Fig.8 Comparison of cavitation shape curves at different water entry depths
入水初期,旋成體接觸水面后頭部產(chǎn)生巨大壓力差,尖錐頭斜面可以更好地將壓差轉化為水的動能產(chǎn)生液體的飛濺,所以尖錐頭體對自由液面擾動最大,濺起的水花區(qū)域最大,平頭體對水面擾動最小,所濺起的水花區(qū)域也最小,半球頭體介于兩者之間。
旋成體頭部外形影響著入水初期瞬間的壓力載荷峰值的大小,如圖9(a)所示為3種頭型在入水初速100 m/s時的壓力峰值曲線。由圖可知壓力峰值出現(xiàn)在入水的初期,即當旋成體頭部接觸到水面時,沖擊載荷達到峰值,但壓力峰值持續(xù)時間較短,隨著入水深度增加,旋成體被空泡包裹,頭部壓力逐漸趨于穩(wěn)定。在相同的入水速度下,可以看出不同的旋成體頭型的局部載荷具有不同的特點,平頭體的瞬間沖擊壓力峰值要明顯大于另外2種頭型的旋成體,半球頭體次之,尖錐頭體的沖擊壓力峰值最小。平頭體的峰值寬度最窄,半球頭次之,尖錐頭峰值寬度最大,由此可以得出旋成體頭部壓力峰值與峰值寬度成反比例的關系。入水一段時間后,3種頭型旋成體頭部的峰值逐漸趨于穩(wěn)定,入水沖擊載荷的波動情況表現(xiàn)基本一致。
如圖9(b)所示為不同頭型旋成體撞擊水面后速度變化曲線,當未接觸到水面時,3種頭型受重力加速度作用,三條速度曲線的變化趨勢基本一致。當旋成體接觸到水面時,可以很明顯的看到平頭體在撞水的瞬間速度曲線出現(xiàn)明顯的拐點,受到的阻力和加速度最大,在接觸水面的初期,尖錐頭體的速度要比半球頭體的速度衰減的慢,但隨著入水深度的增加,尖錐頭是速度衰減率要大于半球頭體。在整個計算時間域內半球頭旋成體速度曲線過渡較為平滑,速度逐漸下降,近似于線性速度下降趨勢。隨著時間推移,旋成體逐漸被空泡所包裹,3種頭型旋成體的速度下降速率趨于相似。
圖9 不同頭型入水參數(shù)曲線Fig.9 Comparison of water entry parameters of different head types
本小節(jié)探究不同頭型旋成體在不同入水速度下空泡演化及沖擊載荷特性,入水初速度由低速向高速過渡,具體計算工況參數(shù)見表3所示。
表3 入水速度工況參數(shù)
如圖10所示為不同入水初速度下同一時刻的無量綱入水空泡演化過程。由圖可知,相同入水速度下,在入水時間相同時,不同頭型的入水深度存在明顯差別。隨著入水初速度的提高,自由液面濺起的水花速度方向逐漸偏離彈體軸線向外擴展。入水初速度越大,液體射流速度越快,自由液面受到的擾動越大,濺起的水花面積越大。因為航行體入水初速度越大,自身的動能越大,撞擊水面后液體獲得的動能與壓力勢能也越大,所以液體飛濺范圍擴大。同一入水時刻下半球頭旋成體的入水深度要大于尖錐頭與平頭旋成體的入水深度。
圖10 不同初速度下相同時刻空泡演化過程曲線Fig.10 Comparison of cavitation evolution at different initial velocities at the same time
對比旋成體在不同初速度下入水過程的衰減速率,由于初始速度相差較大,不易對比觀察,故采用歸一化方法。通過對比各個時刻旋成體速度與入水初速度的比值,繪制出不同初速度下旋成體相對速度隨時間的衰減變化率,如圖11所示半球頭旋成體在不同入水初速度下相對速度變化率曲線??梢钥闯?,在旋成體接觸水面后,不同初始速度下旋成體相對速度變化率出現(xiàn)明顯的拐點,入水初速度越大,相對速度的衰減率越大,呈現(xiàn)正比例的關系。航行體接觸水面的瞬間,所受到的反向加速度急劇增加,這個過程的持續(xù)時間很短,且入水初速度越快,這個過程持續(xù)的時間就越短,旋成體更快被空泡所包裹。隨著入水深度的增加,在浮力、慣性力與重力的共同作用下,旋成體相對衰減速率逐漸趨于穩(wěn)定,速度按一定比例的下降。半球頭航行體在初速度50 m/s和100 m/s時的相對速度衰減曲線間隔要大于 100 m/s和150 m/s時曲線間隔,并且變化明顯。
圖11 半球頭旋成體相對速度變化率曲線Fig.11 Relative velocity change rate of a blunt head body
對于同一頭型的旋成體而言,入水初速度決定了瞬間沖擊壓力峰值的大小。從圖12可以看出,半球頭旋成體瞬間壓力峰值隨著入水初速度的增加而變大;對比初速度為50 m/s、100 m/s、150 m/s時不同頭型旋成體壓力峰值,可知相同頭型在150 m/s速度下的壓力峰值約為50 m/s速度下的8~10倍,100 m/s速度下的壓力峰值約為50 m/s速度下的3~5倍,可知旋成體在入水所受到的壓力峰值隨著初始入水速度的增加呈現(xiàn)近似平方關系的正相關增長。將壓力大小為峰值80%以內載荷分布的寬度定義為沖擊壓力峰值的幅值寬度,可以看到,隨著入水速度的增大,壓力峰值的幅值寬度也會變大。
圖12 半球頭旋成體無量綱壓力變化曲線Fig.12 Curve of dimensionless pressure of a blunt head body
圖13為不同入水初速度下半球頭旋成體的阻力系數(shù)變化曲線,由圖可得當旋成體接觸自由液面時阻力系數(shù)達到峰值且3種速度下峰值基本相近。隨著旋成體入水深度的增加,氣相與液相的流場基本趨于穩(wěn)定,同時不同初速度時的阻力系數(shù)趨于穩(wěn)定??梢钥闯鋈胨跛俣葘π审w的阻力系數(shù)影響不大。
為了更好的分析不同頭型旋成體在不同初速度下垂直入水時的沖擊壓力大小,提取不同頭型旋成體不同初速度下的沖擊壓力峰值及其在旋成體進入空泡后的穩(wěn)定值,提取所得具體數(shù)值如表4所示。
圖13 不同入水速度阻力系數(shù)變化曲線Fig.13 Variation curve of resistance coefficient at different water entry velocities
表4 不同入水速度下沖擊壓力峰值與穩(wěn)定值Table 4 Peak value and stable value of pressure at different water entry velocities
對于不同頭型旋成體而言,初速度為50 m/s時,旋成體所受到的沖擊壓力峰值及其入水之后的穩(wěn)定值最小。不同頭型在不同初速度下頭部沖擊壓力峰值的無量綱變化情況如圖14所示。由圖可知,隨著入水初速度的增大,平頭體的沖擊壓力峰值與穩(wěn)定值始終大于半球頭體和尖錐頭體。隨著旋成體入水速度增大,尖錐頭體與半球頭體頭部沖擊壓力峰值的變化幅度相較平頭體變化不明顯。
圖14 不同頭型旋成體無量綱壓力峰值隨速度的變化曲線Fig.14 Variation of dimensionless pressure peak value with velocities for different head-shaped bodies
采用數(shù)值模擬方法,研究了3種不同頭型旋成體初速度從低速到高速入水過程中空泡形態(tài)與沖擊載荷特性,得到以下結論:
1) 平頭旋成體入水過程形成的空泡半徑最大,半球頭次之,尖錐頭空泡半徑最??;旋成體入水接觸自由液面瞬間,頭部產(chǎn)生巨大壓力差,尖錐頭斜面可以更好將壓差轉化為水的動能將水排開,所以尖錐頭旋成體在入水過程中濺起的水花面積最大,半球頭旋成體次之,平頭旋成體濺起的水花面積最小。
2) 同一入水時刻下,半球頭旋成體的入水深度最深,尖錐頭次之,平頭體入水深度最淺;不同頭型旋成體從低速到高速的入水過程中,將旋成體不同初速度下入水速度變化進行歸一化處理,3種頭型旋成體的速度衰減特性表現(xiàn)基本一致,旋成體初速度越大,入水過程中速度衰減率就越大。
3) 不同頭型旋成體所受到的沖擊壓力峰值隨著入水初速度增加呈現(xiàn)近似平方關系的正相關增長,平頭旋成體沖擊壓力峰值變化幅度受入水初速的影響較另外2種頭型旋成體變化更加劇烈。隨著入水初速度的升高,沖擊壓力峰值的幅值寬度也會變寬,且入水初速度由低速到高速的變化對旋成體的阻力系數(shù)影響不大。