基于數(shù)值波浪的回轉(zhuǎn)體斜侵入水研究

郭張霞，羅 鵬，范光明，郝一凝

(1.中北大學(xué) 機(jī)電工程學(xué)院， 山西 太原 030051；2.中國(guó)第一重型機(jī)械股份公司，遼寧 大連 116113)

導(dǎo)彈、子彈等從高空進(jìn)入水中的現(xiàn)象稱為入水問題[1]。物體入水運(yùn)動(dòng)是一個(gè)涉及固、液、氣三相問題的三維非定常過程，并且在這一過程中涉及到自由面、波浪面的破碎，水空化和相變，受沖擊載荷引起的彈體形變和斷裂等復(fù)雜物理過程[2]。在當(dāng)前的國(guó)際環(huán)境下，研制高速突防的水下攻擊武器，對(duì)我國(guó)海洋防務(wù)具有重要意義，因而吸引著眾多的專家、學(xué)者投入大量精力對(duì)其進(jìn)行研究[3]。邱海強(qiáng)等對(duì)旋轉(zhuǎn)體垂直入水的沖擊載荷和空泡的形成進(jìn)行了研究，得到入水沖擊載荷峰值發(fā)生在入水的初期，且與速度的大小成正比，形成的空泡直徑、長(zhǎng)度則是按照平頭、尖頭和圓頭的順序依次減小，但是其研究并未涉及到旋轉(zhuǎn)體的傾斜入水[4]；魏照宇等對(duì)圓盤尖拱頭形回轉(zhuǎn)體垂直入水進(jìn)行了數(shù)值仿真，得到了其入水超空泡和垂直入水初期的彈道特性[5]；李佳川等則是對(duì)不同擾動(dòng)角下的射彈入水進(jìn)行了研究，得到了彈體入水軌跡、速度、俯仰角和俯仰角速度的變化規(guī)律，研究是在靜水條件下完成的，并沒有涉及到波浪存在的情況[6]；楊衡等利用試驗(yàn)的方法對(duì)不同頭型彈體低速入水的空泡現(xiàn)象進(jìn)行了研究，得到不同角度的錐頭彈在相同速度下入水產(chǎn)生不同狀態(tài)的空化現(xiàn)象[7]，其研究也是在靜水中完成的，沒考慮到波浪因素；然而絕對(duì)的靜水現(xiàn)象在自然界是不存在的，水面都會(huì)存在不同程度的波浪，筆者主要針對(duì)平頭型的回轉(zhuǎn)體傾斜侵入波浪過程中產(chǎn)生的空泡現(xiàn)象和波浪對(duì)平頭型回轉(zhuǎn)體傾斜入侵波浪過程中誘導(dǎo)氣相速度變化規(guī)律，以及回轉(zhuǎn)體頭部所受的壓強(qiáng)變化規(guī)律展開研究。研究結(jié)果可為彈丸入水的流動(dòng)特性提供一定的參考。

1 數(shù)學(xué)模型

1.1 控制方程

在文中將流體假設(shè)為不可壓縮流體，利用VOF(Volume of fluid)均質(zhì)多相流模型對(duì)旋轉(zhuǎn)體傾斜侵入波浪過程中氣相和液相的流動(dòng)進(jìn)行描述，VOF模型是將流場(chǎng)中的各部分流體看作是密度可以改變的單項(xiàng)流體介質(zhì)。通過對(duì)流體介質(zhì)混合連續(xù)性方程和動(dòng)量守恒方程的求解來獲得空泡流場(chǎng)的數(shù)值模擬。

連續(xù)性方程和動(dòng)量守恒方程[8]為

(1)

(2)

ρm=α1ρ1+αvρv+αgρg,

(3)

μm=α1μ1+αvμv+αgμg,

(4)

式中：xi和xj均代表笛卡爾坐標(biāo)，其i、j分別取1、2、3；t為時(shí)間；ρm為混合流體的密度；μm為混合流體的動(dòng)力粘度；μi、μj為旋轉(zhuǎn)體在i、j方向上的速度分量；ρ1、ρv、ρg分別為水相、蒸汽相、空氣相的密度；μ1、μv、μg分別為水相、蒸汽相、空氣相的粘度。

采用的湍流模型為RNG(重正化群)k-ε模型[9]，湍動(dòng)能k方程和湍能耗散率ε方程與標(biāo)準(zhǔn)k-ε方程相似，即：

(5)

(6)

式中：

(7)

(8)

(9)

(10)

式中：Gk為平均速度梯度所引起的湍動(dòng)能生成項(xiàng)；經(jīng)驗(yàn)常數(shù)C1ε、C2ε、η0、β取值分別為1.42、1.68、4.377、0.012.

1.2 計(jì)算模型、區(qū)域及計(jì)算條件

本文中所用的計(jì)算模型如圖1所示，頭型選擇平頭型，長(zhǎng)度a為40 mm，空化器直徑b為9 mm.

文中采用二維模型，將三維問題轉(zhuǎn)化為二維問題，以達(dá)到提高計(jì)算效率的目的。邊界條件的設(shè)置與計(jì)算域的大小如圖2所示，不同頭型旋轉(zhuǎn)體以60°射角侵入波浪；左側(cè)邊界為壓力入口，設(shè)置為壓力輸入，并利用profile文件賦予彈丸相應(yīng)的速度輸入；右側(cè)邊界為壓力出口，設(shè)置為自由流出；上下邊界設(shè)置為壁面。計(jì)算域的水平長(zhǎng)度為750 mm，水深為900 mm，上部高度為350 mm，旋轉(zhuǎn)體的中心坐標(biāo)為(70,1 128.8)(坐標(biāo)原點(diǎn)為左下角頂點(diǎn))。

1.3 網(wǎng)格劃分

由于在計(jì)算過程中采用動(dòng)網(wǎng)格中的局部網(wǎng)格重組法，故而采用非結(jié)構(gòu)網(wǎng)格對(duì)計(jì)算模型進(jìn)行網(wǎng)格劃分。為了保證不同相交面的空泡流場(chǎng)的計(jì)算精度，并且由于存在網(wǎng)格無關(guān)性，故對(duì)旋轉(zhuǎn)體附近進(jìn)行加密處理，其網(wǎng)格模型如圖3所示，網(wǎng)格數(shù)據(jù)為：網(wǎng)格數(shù)143 980個(gè)，節(jié)點(diǎn)數(shù)為48 289個(gè)，最大網(wǎng)格面積為7.73 mm2，最小網(wǎng)格面積為0.787 mm2.

1.4 波浪設(shè)置

海面或者江面等水面由于受到空氣流動(dòng)的影響，空氣與水面相對(duì)運(yùn)動(dòng)而產(chǎn)生摩擦力，摩擦力導(dǎo)致平靜的海面發(fā)生形變，加之重力會(huì)使水面產(chǎn)生回復(fù)原有狀態(tài)的趨勢(shì)，因此水面會(huì)不斷上下起伏，進(jìn)而形成風(fēng)浪，風(fēng)浪的等級(jí)受風(fēng)對(duì)海洋結(jié)構(gòu)物的影響而決定的，根據(jù)蒲福風(fēng)級(jí)[10]，選取輕風(fēng)作為本文所用波浪，其具體信息如表1所示。

表1 蒲福風(fēng)級(jí)表(2級(jí))

在ANSYS-Fluent 15.0中采用明渠流邊界條件模擬明渠效應(yīng)進(jìn)行造波，采用Open Channel Wave BC來設(shè)置入射波，波浪參數(shù)按照表1波浪參數(shù)來設(shè)置，其波形如圖3所示，波形與理論波形相吻合，說明這種造波方式是可行的。

1.5 計(jì)算方法

在ANSYS-Fluent 15.0中利用profile文件和動(dòng)網(wǎng)格技術(shù)賦予彈丸以60°射角和100 m/s的勻速運(yùn)動(dòng)的運(yùn)動(dòng)狀態(tài)，湍流模型選用上述的RNG(重正化群)k-ε模型。綜合考慮計(jì)算的效率和精度，壓力和速度耦合方法采用SIMPLEC算法，Green-Gauss Node Based壓力方程利用PRESOT!方法離散，動(dòng)量方程、湍動(dòng)能耗散率等均采用一階迎風(fēng)方程(First Order Upwind)來計(jì)算。

2 計(jì)算結(jié)果與分析

2.1 平頭型回轉(zhuǎn)體入水過程

圖4是平頭型旋轉(zhuǎn)體入水過程，在圖4中紅色(下部)為水相，藍(lán)色(上部)為空氣相。在回轉(zhuǎn)體運(yùn)動(dòng)前，即t≤300 ms時(shí)，為回轉(zhuǎn)體傾斜侵入波浪的造波期，保證回轉(zhuǎn)體能夠傾斜侵入到波浪的前波面。時(shí)間超過300 ms后，回轉(zhuǎn)體開始運(yùn)動(dòng)，302.78 ms時(shí)，旋轉(zhuǎn)體開始侵入波浪的前波面，在旋轉(zhuǎn)體侵入波面后，在波浪力、重力以及旋轉(zhuǎn)體產(chǎn)生的沖擊波的作用下，產(chǎn)生以旋轉(zhuǎn)體中心不對(duì)稱的水霧化區(qū)，并且在波浪力的影響下，左側(cè)水面產(chǎn)生向上的運(yùn)動(dòng)，右側(cè)的水面則延續(xù)波浪力的影響，繼續(xù)向右運(yùn)動(dòng)并形成霧化。在旋轉(zhuǎn)體入水2 ms旋轉(zhuǎn)體開始脫離，而浸入形成的水霧化區(qū)在波浪力的作用下繼續(xù)增大，并在重力的影響下逐漸下降。

2.2 回轉(zhuǎn)體入水時(shí)氣相變化規(guī)律

當(dāng)回轉(zhuǎn)體開始運(yùn)動(dòng)后，它會(huì)將本身的動(dòng)能傳遞給附近的氣體，并使氣體以一定的速度伴隨回轉(zhuǎn)體運(yùn)動(dòng)。氣體運(yùn)動(dòng)最大速度曲線如圖5所示，由曲線可知：隨著回轉(zhuǎn)體速度的增加，氣體的最大速度也隨之增加，在回轉(zhuǎn)體速度穩(wěn)定后，其誘導(dǎo)氣相速度也在110 m/s大小附近波動(dòng)，其結(jié)果符合文獻(xiàn)[3]中所得到的“誘導(dǎo)氣相速度是入水速度的1.00～1.36倍”這一結(jié)果。當(dāng)回轉(zhuǎn)體繼續(xù)運(yùn)動(dòng)，傾斜入侵波浪時(shí)，誘導(dǎo)氣相速度在波浪力的作用下迅速增加，最大可達(dá)到250 m/s，為回轉(zhuǎn)體入水速度的2.5倍，在回轉(zhuǎn)體完全進(jìn)入波浪后，誘導(dǎo)氣相速度也開始下降，逐漸下降到118 m/s左右。

回轉(zhuǎn)體入水過程中回轉(zhuǎn)體周圍的速度云圖，如圖6所示。

由圖6(a)可知，當(dāng)回轉(zhuǎn)體速度增加時(shí)，由于氣體壓力的相互作用誘導(dǎo)氣相速度最大值在回轉(zhuǎn)體的尾部。由圖6(b)～(f)可知，回轉(zhuǎn)體侵入波浪后，由于波浪的波浪力，回轉(zhuǎn)體前部的誘導(dǎo)氣體速度增加，并快速向后方流動(dòng)，這也是產(chǎn)生以回轉(zhuǎn)體中心不對(duì)稱的水霧化區(qū)的一個(gè)原因。由圖6(g)可知，當(dāng)回轉(zhuǎn)體完全侵入到波浪后，前部的氣體被擠壓到后部，故而誘導(dǎo)氣相速度最大區(qū)域轉(zhuǎn)移到回轉(zhuǎn)體后部。

2.3 波浪對(duì)回轉(zhuǎn)體的沖擊影響

在回轉(zhuǎn)體傾斜侵入波浪的過程中，回轉(zhuǎn)體頭部會(huì)受到來自空氣阻力、水的阻力以及波浪的波浪力的作用，回轉(zhuǎn)體頭部所受到的壓強(qiáng)曲線如圖7所示。由圖7可知，在回轉(zhuǎn)體入水之前，回轉(zhuǎn)體主要受到空氣的阻力作用，由于回轉(zhuǎn)體速度并不大，故而在前期回轉(zhuǎn)體頭部所受壓強(qiáng)較小，在回轉(zhuǎn)體開始侵入波浪時(shí)，回轉(zhuǎn)體頭部受到的壓強(qiáng)急劇增加，最大達(dá)到6.87 MPa，由于波浪的波浪力存在，回轉(zhuǎn)體頭部所受到的壓強(qiáng)呈現(xiàn)上下波動(dòng)的狀態(tài)。

3 結(jié)束語

利用平頭型回轉(zhuǎn)體為計(jì)算模型，在ANSYS-Fluent 15.0中以動(dòng)網(wǎng)格技術(shù)和profile嵌入的方法賦予平頭型回轉(zhuǎn)體以60°射角和100 m/s的速度傾斜侵入波浪進(jìn)行了數(shù)值模擬，得到如下結(jié)果：

1)平頭型回轉(zhuǎn)體在傾斜侵入波浪過程中產(chǎn)生了以回轉(zhuǎn)體中心不對(duì)稱的水霧化區(qū)，并且回轉(zhuǎn)體中心上部的空泡區(qū)域更大。

2)回轉(zhuǎn)體侵入波浪過程中產(chǎn)生的誘導(dǎo)氣相速度在回轉(zhuǎn)體侵入波浪之前呈線性增加，之后以回轉(zhuǎn)體侵入速度的1.1倍大小波動(dòng)；在回轉(zhuǎn)體侵入波浪的過程中，誘導(dǎo)氣相速度急劇增加，最大值達(dá)到回轉(zhuǎn)體速度的2.5倍；回轉(zhuǎn)體進(jìn)入波浪后誘導(dǎo)氣相速度穩(wěn)定在回轉(zhuǎn)體速度的1.2倍。

3)回轉(zhuǎn)體頭部所受壓強(qiáng)在回轉(zhuǎn)體侵入波浪前期很??；在回轉(zhuǎn)體入侵波浪的過程中回轉(zhuǎn)體頭部所受壓強(qiáng)迅速增加到最大；回轉(zhuǎn)體侵入波浪后，回轉(zhuǎn)體頭部所受壓強(qiáng)大小變化成波浪狀態(tài)。