周 杰，徐勝利

(清華大學(xué)航天航空學(xué)院，北京 100084)

彈丸入水特性的SPH計(jì)算模擬*

周 杰，徐勝利

(清華大學(xué)航天航空學(xué)院，北京 100084)

應(yīng)用SPH方法研究彈丸入水過(guò)程中的動(dòng)力學(xué)特征。利用拉格朗日形式的N-S方程自編SPH程序，建立彈丸入水的計(jì)算模型，賦予相應(yīng)的材料參數(shù)及狀態(tài)方程，研究彈丸外形、入水速度和角度等因素對(duì)入水過(guò)程的影響。模擬結(jié)果表明：空化泡的形態(tài)及發(fā)展規(guī)律主要由彈丸的運(yùn)動(dòng)姿態(tài)決定；彈道越穩(wěn)定，阻力因數(shù)就越小，彈丸的存速就越大。SPH方法具有較強(qiáng)的自適應(yīng)性，適用于研究彈丸入水的流固耦合問(wèn)題。

流體力學(xué)；水下彈道；SPH；阻力因數(shù)

高速?gòu)椡枞胨诜礉搾呃椎确矫嬗兄匾獞?yīng)用前景。該問(wèn)題屬于典型的流固耦合問(wèn)題，涉及自由面破碎及其追蹤、水空化和相變、受沖擊載荷作用引起的彈丸形變和斷裂、彈水動(dòng)態(tài)流固耦合(水動(dòng)彈性)等復(fù)雜物理過(guò)程[1]。采用簡(jiǎn)化的理論模型只能分析空化泡形態(tài)等非耦合過(guò)程[2]。Logvinovich[3]提出空穴截面擴(kuò)張的獨(dú)立原理，這個(gè)原理是近似的，但用它可以很簡(jiǎn)單地確定各種情況下空化泡的外形。曹偉等[4]通過(guò)實(shí)驗(yàn)研究了高速射彈的自然超空泡形態(tài)和發(fā)展規(guī)律。易文俊等[5]基于Rayleigh-Plesset單一介質(zhì)可變密度混合多相流模型，計(jì)算分析了不同射彈的超空泡減阻特性。安偉光等[6]依據(jù)氣體泄漏規(guī)則建立空泡內(nèi)氣體平衡方程，聯(lián)合空泡截面擴(kuò)展及運(yùn)動(dòng)體姿態(tài)方程，數(shù)值研究了運(yùn)動(dòng)體帶空泡高速入水非定常過(guò)程。數(shù)值模擬非常適合求解上述多個(gè)耦合過(guò)程，其中無(wú)網(wǎng)格的光滑粒子流體動(dòng)力學(xué)(smoothed particle hydrodynamics, SPH)方法適合求解彈丸入水過(guò)程中介質(zhì)大變形和流固耦合問(wèn)題[7-8]。

本文中基于SPH方法，分析彈丸外形、入水速度和角度等因素對(duì)彈丸入水過(guò)程的影響，對(duì)彈丸入水過(guò)程中產(chǎn)生的空泡形態(tài)、彈道軌跡和彈丸阻力因數(shù)等進(jìn)行模擬計(jì)算，研究彈丸入水的機(jī)理，為分析彈丸穩(wěn)定性和阻力等提供理論依據(jù)。

1 計(jì)算模型和數(shù)值方法

1.1 計(jì)算模型

設(shè)計(jì)3種外形的彈丸，以不同角度、速度入水，研究其入水過(guò)程的動(dòng)力學(xué)特性及水中彈道規(guī)律。圖1給出了3種彈丸及入水模型的示意圖。彈丸以初速度v0、角度θ進(jìn)入H×L的水域的計(jì)算模型，水域的邊界條件設(shè)置為固壁邊界。彈丸有3種類型：(A) 平頭彈丸的長(zhǎng)度L1=40 mm，直徑D=12 mm；(B) 尖頭彈丸的長(zhǎng)度L2=20 mm、L3=20 mm，直徑D=12 mm；(C) 截?cái)囝^彈丸的長(zhǎng)度L2=20 mm、L4=10.44 mm，彈丸直徑D=12 mm，截?cái)嗝嬷睆紻1=6 mm。彈丸入水初速度v0分別為1 200、200 m/s，入水角度θ取90°、60°、30°。

圖1 彈丸外形及入水示意圖Fig.1 Schematic diagram of the projectile shape and projectile entry into the water

1.2 狀態(tài)方程

沖擊載荷作用下，通常采用 Mie-Grüneisen狀態(tài)方程來(lái)描述水的動(dòng)力學(xué)性質(zhì)[9]。水的狀態(tài)方程取決與水的狀態(tài)，在壓縮和膨脹狀態(tài)下水的壓力分別為：

(1)

式中：p為水的壓力；ρ0是初始密度；η是擾動(dòng)前后的密度比，μ=η-1，當(dāng)μ>0時(shí)，水處于壓縮狀態(tài)，當(dāng)μ<0時(shí)，水處于膨脹狀態(tài)；γ0為Grüneisen常數(shù)，a為體積修正系數(shù)；c為水中的聲速，S1、S2、S3均為實(shí)驗(yàn)擬合系數(shù)。水的相關(guān)參數(shù)如表1所示。

表1 Mie-Grüneisen狀態(tài)方程的材料參數(shù)Table 1 Material parameters of Mie-Grüneisen equation of state

彈丸假定為理想彈塑性體，靜水壓力與體積變化率之間呈線性變化，相關(guān)方程如下[10]：

dps=KdεV

(2)

式中：ps為靜水壓力；εV為材料應(yīng)變；K為體積模量，K=E/[3(1-2ν)]，E為彈性模量，ν為泊松比。彈丸的密度為7 800 kg/m3，彈性模量為210 GPa，泊松比為0.3。

1.3 SPH方法

光滑粒子流體動(dòng)力學(xué)(SPH)方法是由L.B.Lucy等[11]提出的一種解決天體物理學(xué)問(wèn)題的純拉格朗日方法，后來(lái)該方法被用于研究介質(zhì)大變形及流固耦合問(wèn)題[12]，在求解過(guò)程中展示出無(wú)網(wǎng)格方法所特有的強(qiáng)大的自適應(yīng)性。SPH方法的核心思想：用一系列任意分布的粒子來(lái)表示問(wèn)題域，用積分表示法來(lái)近似場(chǎng)函數(shù)，應(yīng)用粒子來(lái)對(duì)核近似方程進(jìn)一步近似，在每個(gè)時(shí)間步內(nèi)都要進(jìn)行粒子近似的過(guò)程，將粒子近似法應(yīng)用于描述場(chǎng)函數(shù)的偏微分方程，得到只與時(shí)間相關(guān)的離散化形式的常微分方程，應(yīng)用積分法來(lái)求解常微分方程，得到粒子的場(chǎng)變量。以上特點(diǎn)結(jié)合起來(lái)，使得SPH方法成為具有無(wú)網(wǎng)格、自適應(yīng)、穩(wěn)定等特點(diǎn)以及拉格朗日性質(zhì)的動(dòng)力學(xué)問(wèn)題求解方法，非常適合求解彈丸入水過(guò)程中介質(zhì)大變形和流固耦合問(wèn)題。

與有限元方法不同，SPH方法的節(jié)點(diǎn)是離散的。該方法采用光滑長(zhǎng)度內(nèi)的粒子代替有限元方法中的影響節(jié)點(diǎn)集。因?yàn)闆](méi)有固定的單元體連接，SPH方法中的每個(gè)粒子周圍光滑長(zhǎng)度內(nèi)的粒子個(gè)數(shù)和分布是不固定的，采用J.J.Monaghan[13]提出的B-樣條函數(shù)的光滑函數(shù)，用積分表示法來(lái)近似某一點(diǎn)的場(chǎng)函數(shù)值。為了使SPH算法適合模擬沖擊問(wèn)題，引入Monaghan型人工黏度[14]。彈丸入水問(wèn)題中水域有2種不同邊界：自由邊界和固壁邊界。由于SPH方法的自適應(yīng)性，自由液面處的粒子可以自動(dòng)跟隨液面運(yùn)動(dòng)，不需要做特殊處理。固壁邊界采用排斥邊界條件，在壁面外設(shè)置一組虛擬粒子，虛擬粒子參與流體粒子計(jì)算，粒子自身密度也不斷更新，但位置和速度保持不變[14]。采用蛙跳法求解具有拉格朗日性質(zhì)的SPH方法的N-S方程[12]，時(shí)間步長(zhǎng)滿足CFL條件[15]。

2 計(jì)算結(jié)果分析和討論

本文中研究3種外形的彈丸，以不同初速度v0、角度θ入水，計(jì)算空泡形態(tài)、彈道軌跡、及速度的變化規(guī)律，了解影響彈丸動(dòng)力學(xué)特征的相關(guān)因素。

2.1 彈道軌跡

圖2(a) 尖頭彈丸入水的空化泡形狀發(fā)展(θ=90°)Fig.2(a) Shape formation of cavitation bubble during the cuspidal projectile entry into the water (θ=90°)

圖2給出了尖頭彈丸以初速度v0=1 200 m/s，不同角度θ入水的空化泡形狀的發(fā)展規(guī)律。由圖可知：(1)彈丸垂直入水時(shí)，空化泡的長(zhǎng)度、寬度是隨著運(yùn)動(dòng)時(shí)間(t)增加而變大的，空化泡的長(zhǎng)度的變化率要大于空化泡寬度，各時(shí)刻的空化泡形態(tài)保持幾何相似，且關(guān)于y軸對(duì)稱，在計(jì)算時(shí)間內(nèi)空化泡的沒(méi)有出現(xiàn)閉合現(xiàn)象，如圖2(a)所示；(2)彈丸以θ為60°和30°入水時(shí)，空化泡的長(zhǎng)度、寬度隨著運(yùn)動(dòng)時(shí)間增加而變大，空化泡沒(méi)有出現(xiàn)閉合現(xiàn)象，由于受到偏轉(zhuǎn)力矩作用，彈丸的運(yùn)動(dòng)軌跡是弧形的，導(dǎo)致空化泡不再是軸對(duì)稱的形狀，各時(shí)刻空化泡頭部的形狀、及運(yùn)動(dòng)規(guī)律主要是由彈丸姿態(tài)決定的，如圖2(b)和2(c)所示。

圖2(b) 尖頭彈丸入水的空化泡形狀發(fā)展(θ=60°)Fig.2(b) Shape formation of cavitation bubble during the cuspidal projectile entry into the water (θ=60°)

圖2(c) 尖頭彈丸入水的空化泡形狀發(fā)展(θ=30°)Fig.2(c) Shape formation of cavitation bubble during the cuspidal projectile entry into the water (θ=30°)

圖3是尖頭彈丸以v0=1 200 m/s，θ=60°入水的計(jì)算結(jié)果。由圖可知：(1)t=0.1 ms時(shí)，彈丸高速進(jìn)入水域，彈丸的運(yùn)動(dòng)方向保持與入水前一致；(2)t=0.2和0.3 ms時(shí)，由于彈丸入水角度θ=60°，因此在水中運(yùn)動(dòng)時(shí)會(huì)受到一個(gè)偏轉(zhuǎn)力矩作用，引起彈丸的運(yùn)動(dòng)姿態(tài)發(fā)生變化，導(dǎo)致彈丸運(yùn)動(dòng)方向相對(duì)于入水前偏轉(zhuǎn)了較大的角度(彈道失穩(wěn))，且彈丸周圍產(chǎn)生空化泡；(3)t=0.5和0.7 ms時(shí)，隨著彈丸的穩(wěn)定性的下降，彈道偏轉(zhuǎn)現(xiàn)象越來(lái)越明顯，彈丸周圍形成了一個(gè)不閉合的弧形空泡，且速度急劇下降。

圖3 尖頭彈丸的入水軌跡(v0=1 200 m/s，θ=60°)Fig.3 Trajectory of cuspidal projectile entry into the water (v0=1 200 m/s，θ=60°)

圖4是平頭彈丸以v0=1 200 m/s，θ=90°入水的計(jì)算結(jié)果。由圖可知：(1)t=0.08 ms時(shí)，彈丸高速進(jìn)入水域，會(huì)在彈頭邊界形成應(yīng)力集中，導(dǎo)致邊界產(chǎn)生裂紋；(2)t=0.2 ms時(shí)，彈頭邊界部分出現(xiàn)斷裂現(xiàn)象，破片與彈丸分離，彈頭變化為錐形凸臺(tái)，且彈丸周圍有空化泡產(chǎn)生；(3)t=0.4和0.6 ms時(shí)，隨著彈丸在水中運(yùn)動(dòng)，空泡半徑越來(lái)越大，產(chǎn)生的破片以較大的速度向兩側(cè)運(yùn)動(dòng)，同樣也會(huì)產(chǎn)生空化泡，對(duì)彈丸周圍的空化泡產(chǎn)生影響；(4)t=0.8 ms時(shí)，彈丸的運(yùn)動(dòng)速度越來(lái)越小，彈丸周圍形成閉合空泡。

圖4 平頭彈入水運(yùn)動(dòng)軌跡(v0=1 200 m/s，θ=90°)Fig.4 Trajectory of blunt projectile entry into the water (v0=1 200 m/s，θ=90°)

圖5(a) 彈丸入水過(guò)程中的彈道軌跡(θ=90°)Fig.5(a) Ballistic trajectory of projectile during the process of entry into the water (θ=90°)

圖5給出了3種彈丸在高速和低速下的彈道軌跡，其中：A、B、C分別代表平頭、尖頭和截?cái)囝^的彈丸。圖5(a)為彈丸垂直入水的彈道軌跡：平頭彈丸在高、低速條件下，彈道軌跡均保持良好的穩(wěn)定性；尖頭彈丸在高速情況下彈道穩(wěn)定性良好，而低速情況下彈道則失穩(wěn)；截?cái)囝^彈丸在高、低速條件下，彈道軌跡均失穩(wěn)。圖5(b)和5(c)分別是彈丸入水角度θ=60°和30°的彈道軌跡：低速情況下，3種彈丸的彈道軌跡均能保證良好的穩(wěn)定性；而在高速情況下彈丸的彈道軌跡的穩(wěn)定性均會(huì)變的不理想，尤其是尖頭彈丸。實(shí)際上，彈丸在水中的彈道軌跡主要是由彈丸的質(zhì)心、壓心的位置關(guān)系決定的，設(shè)計(jì)彈丸外形時(shí)應(yīng)給予考慮，如對(duì)于大長(zhǎng)徑比的彈丸，應(yīng)該考慮通過(guò)擺動(dòng)(尾拍)保持彈丸運(yùn)動(dòng)平衡。

圖5(b) 彈丸入水過(guò)程中的彈道軌跡(θ=60°)Fig.5(b) Ballistic trajectory of projectile during the process of entry into the water (θ=60°)

圖5(c) 彈丸入水過(guò)程中的彈道軌跡(θ=30°)Fig.5(c) Ballistic trajectory of projectile during the process of entry into the water (θ=30°)

圖6 彈丸入水的速度變化規(guī)律Fig.6 Profile of the velocity variation during the projectile entry into the water

2.2 運(yùn)動(dòng)規(guī)律

圖6給出了3種彈丸在不同入水角度下的速度。圖6(a)為彈丸入水速度v0=1 200 m/s的運(yùn)動(dòng)規(guī)律。垂直入水的情況下，尖頭彈丸產(chǎn)生的空化泡有效的降低了水的阻力，而截?cái)囝^彈的彈道軌跡發(fā)生了偏轉(zhuǎn)現(xiàn)象，所以受到的阻力較大，速度下降的最為明顯；θ=60°和30°時(shí)，3種類型的彈丸在水中運(yùn)動(dòng)的穩(wěn)定性均比較差，受到的水的阻力均較大，尤其是尖頭彈丸，速度下降的最為明顯。圖6(b)為彈丸入水速度v0=200 m/s的運(yùn)動(dòng)規(guī)律：垂直入水時(shí)，3種彈丸的入水速度大小關(guān)系為vB>vC>vA；θ=60°和30°時(shí)，則為vB>vC≈vA。

根據(jù)牛頓第二定律，在水中運(yùn)動(dòng)的彈丸質(zhì)心運(yùn)動(dòng)方程為：

(3)

式中：m為彈丸質(zhì)量，v為彈丸運(yùn)動(dòng)速度，Rx為彈丸受到的阻力，Ry為升力，Rz為馬格努斯力，Rg為重力，RC為科氏力升力，RΩ為地球自轉(zhuǎn)產(chǎn)生的慣性力。計(jì)算中只考慮阻力Rx=-0.5ρwv2CxS，ρw是水的密度，Cx是彈丸阻力因數(shù)，S是垂直中心軸的最大截面積。

圖7給出了不同速度下，阻力因數(shù)隨時(shí)間的變化規(guī)律。彈丸入水速度v0=1 200 m/s時(shí)：垂直入水的情況下，平頭彈入水的初始阻力因數(shù)較大，隨著彈丸的破壞，彈頭形狀變?yōu)殄F形凸臺(tái)，有效的降低了阻力因數(shù)，尖頭彈的阻力因數(shù)較小，隨著時(shí)間增長(zhǎng)保持在0.5左右，而截?cái)囝^彈丸的彈道軌跡發(fā)生了偏轉(zhuǎn)現(xiàn)象，所以阻力因數(shù)均較大；入水角度為60°、30°的情況下，3種類型的彈丸的阻力因數(shù)較垂直入水時(shí)大，3種彈丸的阻力因數(shù)從大到小依次為平頭彈、尖頭彈、截?cái)囝^彈。彈丸入水速度v0=200 m/s時(shí)，平頭彈丸的阻力因數(shù)明顯要高于其他的2種彈形；尖頭彈、截?cái)囝^彈的阻力因數(shù)均在1.2附近震蕩。

圖7 彈丸阻力因數(shù)隨時(shí)間的變化規(guī)律Fig.7 Variation of the projectile's drag coefficient with time

3 結(jié) 論

本文中采用具有無(wú)網(wǎng)格、自適應(yīng)、穩(wěn)定以及拉格朗日性質(zhì)的SPH方法，研究了不同外形彈丸以不同初速度、角度入水的動(dòng)力學(xué)過(guò)程：

(1) 彈丸入水時(shí)，如果彈道軌跡穩(wěn)定，會(huì)產(chǎn)生對(duì)稱的空化泡，否則會(huì)產(chǎn)生不規(guī)則的空化泡，空化泡的形態(tài)及發(fā)展規(guī)律主要由彈丸的運(yùn)動(dòng)姿態(tài)決定；

(2) 垂直入水時(shí)，高速條件下平頭和尖頭彈丸的彈道軌跡穩(wěn)定性較好，低速條件下只有平頭彈丸的彈道穩(wěn)定；以一定角度入水時(shí)，3種外形彈丸在低速條件下的彈道穩(wěn)定性要優(yōu)于高速條件；彈道穩(wěn)定性是由彈丸的質(zhì)心、壓心的位置關(guān)系決定的；

(3) 垂直高速入水時(shí)，尖頭彈丸的阻力因數(shù)最小，而以一定角度高速入水時(shí)，截?cái)囝^彈阻力因數(shù)最??；低速入水時(shí)，平頭彈的阻力因數(shù)較大，而尖頭彈、截?cái)囝^彈的阻力因數(shù)大小相當(dāng)；

(4) 彈道穩(wěn)定性越好，阻力因數(shù)就越小，彈丸的存速就越大；阻力因數(shù)的大小還與彈丸頭部空化器的尺寸相關(guān)。

[1] Putilin S I. Some features of dynamics of supercavitating models[J]. Applied Hydromechanics, 2000,2(74):65-74.

[2] Knapp R T, Daily J W, Hammitt F G. Cavitation[M]. NewYork: McGraw-Hill, 1970.

[3] Franc J-P, Michel J-M. Fundamentals of cavitation[M]. The Netherlands: Kluwer Academic Publishers, 2004

[4] 曹偉，王聰，魏英杰，等.自然超空泡形態(tài)特性的射彈試驗(yàn)研究[J].工程力學(xué),2006,23(12):175-187. Cao Wei, Wang Cong, Wei Yingjie, et al. High-speed projectile experimental investigations on the characteristics of natural supercavitation[J]. Engineering Mechanics, 2006,23(12):175-187.

[5] 易文俊,王中原,熊天紅,等.水下高速射彈超空泡減阻特性研究[J].彈道學(xué)報(bào),2008,20(4):1-4. Yi Wenjun, Wang Zhongyuan, Xiong Tianhong, et al. Research on drag reduction characteristics of a underwater high-speed supercavitation projectile [J]. Journal of Ballistics, 2008,20(4): 1-4.

[6] 安偉光,蔣運(yùn)華,安海.運(yùn)動(dòng)體高速入水非定常過(guò)程研究[J].工程力學(xué),2011,28(3):251-256. An Weiguang, Jiang Yunhua, An Hai. The unsteady water entry process study of high-speed vehicle[J]. Engineering Mechanics, 2011,28(3):251-256.

[7] Chen J K, Beraun J E. A generalized smoothed particle hydrodynamic method for nonlinear dynamic problems[J]. Computer Methods in Applied Mechanics and Engineering, 2000,190(1):225-239.

[8] Cleary P W, Prakash M, Ha J. Novel applications of smoothed particle hydrodynamics (SPH) in metal forming[J]. Journal of Materials Processing Technology, 2006,177(1):41-48.

[9] Shin Y S, Lee M, Lam K Y, et al. Modeling mitigation effects of watershield on shock wave[J]. Shock and Vibration, 1998,5(4):225-234.

[10] Libersky L D, Petschek A G. High strain Lagrangian hydrodynamics: A three-dimensional SPH code for dynamic material response[J]. Journal of Computational Physics, 1993,109(1):67-75.

[11] Lucy L B. A numerical approach to the testing of the fission hypothesis[J]. The Astronomical Journal, 1977,82(12):1013-1024.

[12] Liu G R, Liu M B. Smoothed particle hydrodynamics: A meshfree particle method[M]. German: Springer Berlin /Heidelberg, 2004:1-491．

[13] Monaghan J J. Particle methods for hydrodynamic[J]. Computer Physics Report,1985,3(2):71-124.

[14] Monaghan J J. On the problem of penetration in particle menthods[J]. Journal of Computer Physics, 1989,82(1):1-15.

[15] Monaghan J J. Smoothed particle hydrodynamics[J]. Reports on Progress in Physics, 2005,68(8):1703-1759.

(責(zé)任編輯 王小飛)

SPH simulation on the behaviors of projectile water entry

Zhou Jie, Xu Shengli

(SchoolofAerospaceEngineering,TsinghuaUniversity,Beijing100084,China)

In this work we investigated the dynamic behaviors of the projectile water entry using the SPH method. We developed our own SPH program based on the N-S equation of the Lagrange form and established a calculation model for the projectile water entry and, with corresponding material parameters and equation of state given, studied the influence of such factors as projectile shape, velocity and angle into the water on the process of the projectile water entry. The simulation results show that the formation and the development of the cavitation bubble are mainly determined by the projectile's state of motion: the more stable the projectile's trajectory, the smaller its drag coefficient, and the greater its sustained velocity. It is found that the SPH method has a high self-adaptability, for which it is applicable for studying the problems related with fluid-structure interaction occurring during the process of the projectile water entry.

fluid mechanics; underwater trajectory; SPH; drag coefficient

10.11883/1001-1455(2016)03-0326-07

2014-09-22；

2014-12-05

中國(guó)博士后科學(xué)基金面上項(xiàng)目(2015M581081)

周 杰(1986- )，男，博士，Beijihu1986@163.com。

O352國(guó)標(biāo)學(xué)科代碼：13025

A