孔祥韶，石 干，王旭陽(yáng),3，周紅昌，吳衛(wèi)國(guó)

（1. 武漢理工大學(xué)綠色智能江海直達(dá)船舶與郵輪游艇研究中心，湖北 武漢 430063；2. 武漢理工大學(xué)交通學(xué)院，湖北 武漢 430063；3. 西安航天動(dòng)力研究所液體火箭發(fā)動(dòng)機(jī)技術(shù)重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710100）

大型水面艦艇通常采用液艙來(lái)防護(hù)反艦武器爆炸破片對(duì)內(nèi)部結(jié)構(gòu)和人員的威脅[1]。彈體在侵徹液艙過(guò)程中與液體的作用機(jī)理十分復(fù)雜，彈體在液體介質(zhì)中的貫徹距離、速度衰減規(guī)律以及高速?gòu)楏w引起的水錘效應(yīng)，是廣大學(xué)者長(zhǎng)期關(guān)注的問(wèn)題，對(duì)此已開(kāi)展了大量的研究。Lee 等[2]在經(jīng)典沖擊壓力模型的基礎(chǔ)上，將彈體簡(jiǎn)化為點(diǎn)源，對(duì)高速?gòu)楏w在液體中引起的空腔效應(yīng)展開(kāi)了研究，推導(dǎo)建立了空腔動(dòng)力學(xué)方程。李營(yíng)等[3]在考慮長(zhǎng)徑比的影響下，將阻力因數(shù)作為與初速度有關(guān)的常量來(lái)處理，擬合得到了阻力因數(shù)計(jì)算經(jīng)驗(yàn)公式。沈曉樂(lè)等[4]考慮彈體的墩粗變形進(jìn)行了實(shí)驗(yàn)研究，對(duì)方形彈體入水的阻力因數(shù)進(jìn)行了修正。郭子濤[5]分析了彈頭形狀對(duì)彈體水彈道穩(wěn)定性以及速度衰減規(guī)律的影響，提出了一種簡(jiǎn)化的彈體頭型阻力因數(shù)估算公式。Zhao 等[6]在考慮雷諾數(shù)的基礎(chǔ)上對(duì)球形彈入水阻力因數(shù)進(jìn)行擬合，提出了球形彈入水阻力因數(shù)經(jīng)驗(yàn)公式。Zhang 等[7]在Lee 等[2]研究工作的基礎(chǔ)上針對(duì)長(zhǎng)徑比在0.1～0.5 之間的餅形彈進(jìn)行研究，并對(duì)餅形彈入水阻力因數(shù)計(jì)算公式進(jìn)行了修正。這些研究普遍簡(jiǎn)化了阻力因數(shù)的影響因素，與實(shí)際物理過(guò)程有較大的區(qū)別。

本文中，擬采用理論分析和數(shù)值模擬的方法，針對(duì)以上存在的問(wèn)題，開(kāi)展截錐形彈體垂直穿透液體介質(zhì)時(shí)的速度衰減特性研究，重點(diǎn)分析不同彈體頭形因數(shù)對(duì)其速度衰減規(guī)律的影響，以期建立阻力因數(shù)與彈體的瞬時(shí)速度之間的關(guān)系，并引入彈體形狀參數(shù)的影響，建立截錐形彈體在液體中的速度衰減分析模型。

1 彈體侵徹液體介質(zhì)速度衰減規(guī)律

當(dāng)彈體在液體介質(zhì)中高速運(yùn)動(dòng)時(shí)，彈體損失的動(dòng)能主要轉(zhuǎn)化為液體的動(dòng)能及空泡的壓力勢(shì)能，而空泡的存在極大地減小了彈體與液體的接觸面積，從而彈體運(yùn)動(dòng)所受的摩擦阻力減小，因此壓差阻力是影響彈體阻力的主要因素。則根據(jù)牛頓第二定律建立彈體運(yùn)動(dòng)方程如下[8]：

2 數(shù)值計(jì)算

2.1 數(shù)值計(jì)算工況設(shè)計(jì)及模型

本節(jié)設(shè)計(jì)4 種不同頭形因數(shù)的截錐形彈體，采用顯式動(dòng)力學(xué)非線性有限元程序Autodyn 建立彈體的仿真模型并開(kāi)展系列工況的數(shù)值計(jì)算，分析頭形因數(shù)對(duì)彈體垂直入水后速度衰減的影響規(guī)律。

圖1 為截錐形彈體示意圖，定義彈體頭部圓臺(tái)直徑D2與主體最大直徑D1的比值D2/D1為截錐形彈體頭形因數(shù) ψ 。各工況中彈體主體長(zhǎng)度L1=30 mm ，錐體長(zhǎng)度L2=10 mm ，主體直徑D1=12 mm 保持不變，4 種彈體頭部圓臺(tái)直徑分別為0、4、8、12 mm，對(duì)應(yīng)頭形因數(shù) ψ1～ψ4分別為0、1/3、2/3、1。為了控制變量的數(shù)量，不同頭形因數(shù)截錐形彈體的質(zhì)量均調(diào)整為32.5 g，有限元模型如圖2 所示。本次數(shù)值模擬采用二維軸對(duì)稱建模方法，水域網(wǎng)格尺寸為400 mm×400 mm，采用邊長(zhǎng)為1 mm 的四邊形歐拉單元進(jìn)行離散，四周設(shè)置透射邊界條件，模擬無(wú)限水域，彈體用拉格朗日網(wǎng)格離散。

數(shù)值模擬中，采用S h o c k 狀態(tài)方程和Johnson-Cook 本構(gòu)關(guān)系描述彈體材料的動(dòng)態(tài)力學(xué)行為[12]：

圖 1 截錐形彈體示意圖Fig. 1 Schematic diagram of a truncated cone-shaped projectile

式中： σy為材料的動(dòng)態(tài)屈服應(yīng)力；A為材料的靜態(tài)屈服應(yīng)力；B為材料的硬化參數(shù)； εp為等效塑性應(yīng)變；n為硬化指數(shù)；C為應(yīng)變率參數(shù)；ε ˙ 為等效塑性應(yīng)變率；ε ˙0為參考應(yīng)變率，取 ε ˙0=1 s-1；T*=(T-Tr)/(Tm-Tr) ，T為溫度，Tr為室溫，Tm為材料熔化溫度，m為溫度指數(shù)。本文中彈體材料選用45 鋼[13]，A為506 MPa，B為320 MPa，n為0.28，C為0.064，m為1.06，Tm為1 750 K。

圖 2 不同頭形因數(shù)截錐形彈體有限元模型Fig. 2 Finite element models for truncated cone-shaped projectiles with different head type coefficients

表 1 水的狀態(tài)方程各項(xiàng)參數(shù)[14]Table 1 Parameters of equation of state for water[14]

本文中對(duì)每種不同頭形因數(shù)的截錐形彈體各設(shè)置了5 種不同的入水初速度開(kāi)展數(shù)值模擬，共包含20 種工況條件，如表2 所示。

表 2 數(shù)值計(jì)算工況Table 2 Numerical calculation conditions

2.2 有限元方法驗(yàn)證

采用本文的數(shù)值模擬方法對(duì)試驗(yàn)過(guò)程[5]開(kāi)展計(jì)算，通過(guò)對(duì)比分析數(shù)值計(jì)算結(jié)果和試驗(yàn)數(shù)據(jù)來(lái)驗(yàn)證數(shù)值方法對(duì)高速?gòu)楏w入水過(guò)程分析的可靠性和準(zhǔn)確性。文獻(xiàn)[5]中使用輕氣炮發(fā)射圓柱形平頭彈體，該彈體直徑D=12.65 mm，長(zhǎng)度L=25.4 mm，入水初速度為603 m/s，試驗(yàn)中使用高速攝影儀記錄了彈體運(yùn)動(dòng)軌跡及空泡變化，測(cè)量了不同時(shí)刻的空泡尺寸以及彈體速度等數(shù)據(jù)。將數(shù)值模擬得到的高速?gòu)楏w入水過(guò)程中與試驗(yàn)[5]中相應(yīng)的數(shù)據(jù)進(jìn)行對(duì)比，彈體入水后空泡尺寸對(duì)比如圖3 所示，可以發(fā)現(xiàn)在t=0.222 ms和t=0.440 ms 時(shí)刻相同位置處的空泡直徑數(shù)值計(jì)算結(jié)果與試驗(yàn)測(cè)量數(shù)據(jù)吻合較好

圖 3 針對(duì)彈體直徑為12.65 mm，長(zhǎng)度為25.4 mm，入水初速度為603 m/s 時(shí)，試驗(yàn)與數(shù)值計(jì)算的空泡尺寸Fig. 3 Comparison of cavitation sizes obtained experimentally and numerically for the projectile with the diameter of 12.65 cm and the length of 25.4 cm, water entering at 603 m/s

此外，進(jìn)一步開(kāi)展尺寸D=12.65 mm，L=25.4 mm 的平頭圓柱形彈體入水速度603 m/s 及397m/s 工況；D=12.65 mm，L=38.1 mm 的平頭圓柱形彈體入水速度為498 m/s 及414 m/s 工況的數(shù)值計(jì)算，將彈體位移隨時(shí)間的變化以及彈體速度衰減的數(shù)值計(jì)算結(jié)果與相應(yīng)的試驗(yàn)測(cè)量數(shù)據(jù)進(jìn)行對(duì)比，如圖4～5 所示?？梢钥闯觯疚臄?shù)值模擬方法在計(jì)算彈體的運(yùn)動(dòng)特性方面具有較好的精度，可為本文的進(jìn)一步的研究工作開(kāi)展提供了準(zhǔn)確可靠的分析手段。

圖 4 長(zhǎng)度為25.4 mm 的彈體在兩種入水速度工況下位移和速度的變化Fig. 4 Changes of displacement and velocity of the projectile with the length of 25.4 mm at two initial velocities of water entry

圖 5 長(zhǎng)度為38.1 mm 的彈體在兩種入水速度工況下位移和速度的變化Fig. 5 Changes of displacement and velocity of the projectile with the length of 38.1 mm at two initial velocities of water entry

2.3 數(shù)值模擬結(jié)果及分析

采用驗(yàn)證后的數(shù)值模擬方法開(kāi)展表2 中20 種工況的數(shù)值計(jì)算，將頭形因數(shù)不同的截錐形彈體在不同入水速度下的速度衰減與入水距離的關(guān)系繪制于圖6，通過(guò)對(duì)比分析圖中數(shù)據(jù)可以發(fā)現(xiàn)：

（1）當(dāng)截錐形彈體的頭形因數(shù)不變時(shí)，彈體入水后速度變化曲線的斜率隨著入水初始速度的升高而逐漸增大，說(shuō)明入水速度越高，速度衰減越快；

（2）彈體入水初速度不變、頭形因數(shù)增大時(shí)，彈體剩余速度并非單純隨頭形因數(shù)增大而降低，這一特性在v0=1 400 m/s 時(shí)表現(xiàn)最明顯。4 種頭形因數(shù)的彈體均以1 400 m/s 的初速度垂直入水時(shí)，運(yùn)動(dòng)400 mm后彈體剩余速度與初速度的比值vb/v0分別為0.699 0、0.771 0、0.709 9、0.512 7。說(shuō)明入水初速度不變時(shí)，隨著彈體頭形因數(shù)的增大，vb/v0呈現(xiàn)先增大后減小的趨勢(shì)，在 ψ =1/3 時(shí)達(dá)到最大，為0.771。

圖 6 不同頭形因數(shù)的截錐形彈體在不同入水速度工況下速度隨入水距離的衰減Fig. 6 Velocity attenuation of truncated cone-shaped projectiles with different head coefficients with water-entry distance at different initial velocities of water entry

3 速度衰減分析模型

3.1 模型建立

根據(jù)式（4）可以計(jì)算得出截錐形彈體垂直入水后的瞬時(shí)阻力因數(shù)Cd，并將其作為縱坐標(biāo)，將彈體的瞬時(shí)速度vb與入水初速度v0的無(wú)量綱比值作為橫坐標(biāo)，繪制Cd- (vb/v0) 關(guān)系圖。由于篇幅有限，本文中僅繪出ψ=0 的彈體Cd- (vb/v0) 關(guān)系圖，如圖7 所示。從圖7 可以發(fā)現(xiàn)，阻力因數(shù)Cd隨著彈體速度的衰減呈振蕩趨勢(shì)，且當(dāng)截錐形彈體速度較低時(shí)，隨著侵徹過(guò)程的進(jìn)行，速度衰減率出現(xiàn)較小幅度的上升。究其原因是因?yàn)殡S著速度的降低，彈體誘導(dǎo)空泡的尺寸也在不斷減弱，彈液間的接觸面積增大，摩擦阻力增大，導(dǎo)致彈體的速度衰減率和阻力因數(shù)Cd小幅增大。

圖 7 不同入水速度下， ψ =0 的截錐形彈體阻力因數(shù)與 v b/v0 的關(guān)系Fig. 7 Resistance factor varying with v b/v0 for the truncated cone-shaped projectile with ψ=0 at different water-entry velocities

從完整的侵徹過(guò)程來(lái)看，阻力因數(shù)Cd變化過(guò)程較復(fù)雜，阻力因數(shù)隨著彈體速度的降低，呈現(xiàn)如上述圖像所示的震蕩變化趨勢(shì)。在以往的研究方法中將其作為常量來(lái)處理，不具備廣泛的適用性。從簡(jiǎn)化計(jì)算的角度出發(fā)，本文中采用線性關(guān)系進(jìn)行擬合：

式中：a0和a1為未知參數(shù)。

針對(duì)圖7 中的數(shù)據(jù)，采用最小二乘法進(jìn)行擬合，擬合的結(jié)果如圖7 中紅線所示。同時(shí)，對(duì)頭形因數(shù)ψ=1/3 ， ψ =2/3 ， ψ =1 的 截錐形彈體數(shù)據(jù)按照同樣方式進(jìn)行擬合，確定的參數(shù)a0和a1列于表3。

表 3 不同工況下的參數(shù) a0 和 a1 的數(shù)值Table 3 Values of parameters a0 and a1 under different working conditions

水中聲速取1 400 m/s，用最小二乘法對(duì)參數(shù)a0和a1進(jìn)行擬合：

式中：a01、a02、a03、a04、a05、a11、a12、a13、a14、a15為待定參數(shù)，結(jié)果如圖8 所示。

由擬合結(jié)果來(lái)看，四階方程可以較好的擬合a0、a1關(guān)于無(wú)量綱數(shù)M的關(guān)系，可以發(fā)現(xiàn)同一頭形因數(shù)下參數(shù)a0及a1的高度相關(guān)，但變化趨勢(shì)相反，近似關(guān)于某一水平線y=n1對(duì)稱。隨著頭形因數(shù)的增大，參數(shù)a0及a1的變化逐漸趨于復(fù)雜，這是因?yàn)楫?dāng)頭形因數(shù)增大時(shí)，彈體的垂直迎流面積增大，侵徹時(shí)鐓粗變形效應(yīng)變得明顯，導(dǎo)致阻力因數(shù)的變化趨于復(fù)雜。參數(shù)a0、a1的擬合結(jié)果列于表4。

圖 8 對(duì)于頭形因數(shù)不同的截錐形彈體，參數(shù)a0 和a1 與馬赫數(shù) Ma 的關(guān)系Fig. 8 Changes of parameters a0 and a1 with Mach number Ma for the truncated cone-shaped projectiles with different head shape factors

表 4 參數(shù) a0 、 a1 的擬合結(jié)果Table 4 Fitting results of parameters a0 and a1

圖 9 a 系列參數(shù)與頭形因數(shù) ψ 的關(guān)系Fig. 9 Series parameters of a varyied with head shape factor ψ

表 5 ai j 系列參數(shù)的擬合結(jié)果Table 5 Fitting results of ai j series parameters

綜上所述，當(dāng)頭形因數(shù) ψ 不同的截錐形彈體垂直侵徹入水時(shí)可以通過(guò)如下公式計(jì)算實(shí)時(shí)阻力因數(shù)Cd：

3.2 算例分析及驗(yàn)證

通過(guò)2.2 節(jié)中已驗(yàn)證的計(jì)算方法開(kāi)展系列模擬計(jì)算，通過(guò)對(duì)公式（18）的編程求解，得到彈體的速度衰減。將數(shù)值計(jì)算結(jié)果與上述速度衰減模型計(jì)算結(jié)果進(jìn)行對(duì)比，驗(yàn)證其可靠性。

驗(yàn)證計(jì)算中，設(shè)置v0=400 m/s 速度下1/6、1/2、5/6 等3 種頭形因數(shù)工況以及頭形因數(shù)為 ψ =1/6 的彈體在100～600 m/s 區(qū)間內(nèi)6 種入水速度工況進(jìn)行驗(yàn)證計(jì)算，結(jié)果如圖10～11 所示。

圖 10 入水速度為400 m/s，頭形因數(shù)不同的彈丸，速度隨入水距離衰減的模擬結(jié)果與公式計(jì)算結(jié)果的比較Fig. 10 Comparison of velocity attenuation with distance between numerical simulation and formula calculation for the projectiles with different head shape factors at the initial water-entry velocity of 400 m/s

由計(jì)算可以發(fā)現(xiàn)，對(duì)于彈徑相同、入水初速相同而頭形因數(shù)不同的彈體，其速度衰減規(guī)律有著顯著差異。本文的計(jì)算分析模型中考慮了因頭形因數(shù)不同而引起的阻力因數(shù)變化對(duì)彈速衰減的影響，與數(shù)值計(jì)算結(jié)果吻合較好。當(dāng)彈體入水初速為400 m/s，頭形因數(shù)為1/2 時(shí)，數(shù)值計(jì)算數(shù)據(jù)與本文公式計(jì)算結(jié)果之間的誤差約為6.8%，其他兩種頭形因數(shù)工況下誤差更小；在彈體頭形因數(shù)為1/6、不同速度工況中，公式計(jì)算與數(shù)值模擬結(jié)果之間的最大誤差約為8.71%。

圖 11 頭形因數(shù)為1/6，入水速度不同的彈丸，速度隨入水距離衰減的模擬結(jié)果與公式計(jì)算結(jié)果的比較Fig. 11 Comparison of velocity attenuation with distance between numerical simulation and formula calculation for the projectile with the head shape factor of 1/6 at different initial water-entry velocities

從上述對(duì)比結(jié)果可以發(fā)現(xiàn)，本文的速度衰減模型可對(duì)不同頭形因數(shù)的截錐形彈體垂直侵徹液體介質(zhì)過(guò)程中的速度衰減進(jìn)行快速預(yù)報(bào)，彌補(bǔ)了以往研究中將阻力因數(shù)Cd當(dāng)作常量計(jì)算的不足，更符合彈體垂直侵徹液體介質(zhì)的實(shí)際物理過(guò)程。

4 結(jié) 論

本文中提出了與彈體形狀和無(wú)量綱速度相關(guān)的阻力因數(shù)的計(jì)算方法，通過(guò)建立運(yùn)動(dòng)方程得到了截錐形彈體垂直侵徹液體介質(zhì)速度衰減的計(jì)算分析模型。并通過(guò)算例分析，驗(yàn)證了本文計(jì)算方法的可靠性。得到以下結(jié)論：

（1）高速?gòu)楏w在液體介質(zhì)中運(yùn)動(dòng)方程的阻力因數(shù)與彈體形狀和無(wú)量綱速度有關(guān)；該因數(shù)對(duì)彈體高速運(yùn)動(dòng)過(guò)程中與液體的接觸面積進(jìn)行了修正，將彈體頭形因數(shù)和速度因素考慮在內(nèi)，更接近實(shí)際的物理過(guò)程。

（2）將阻力因數(shù)引入運(yùn)動(dòng)方程后得到的計(jì)算結(jié)果與數(shù)值模擬數(shù)據(jù)吻合較好, 說(shuō)明本文提出的分析模型可對(duì)不同頭形因數(shù)的截錐形彈體垂直侵徹液體介質(zhì)過(guò)程中的速度衰減進(jìn)行快速預(yù)報(bào)。

（3）由于彈體侵徹液體時(shí)的作用機(jī)理非常復(fù)雜，且彈體高速侵徹時(shí)會(huì)發(fā)生彈體鐓粗和變形，同時(shí)流體具有可壓縮性，對(duì)彈速衰減影響較大。本文的分析中未考慮其影響，在后續(xù)工作中將進(jìn)一步的深入研究。