孔祥韶,石 干,王旭陽,3,周紅昌,吳衛(wèi)國
(1. 武漢理工大學(xué)綠色智能江海直達(dá)船舶與郵輪游艇研究中心,湖北 武漢 430063;2. 武漢理工大學(xué)交通學(xué)院,湖北 武漢 430063;3. 西安航天動力研究所液體火箭發(fā)動機(jī)技術(shù)重點(diǎn)實驗室,陜西 西安 710100)
大型水面艦艇通常采用液艙來防護(hù)反艦武器爆炸破片對內(nèi)部結(jié)構(gòu)和人員的威脅[1]。彈體在侵徹液艙過程中與液體的作用機(jī)理十分復(fù)雜,彈體在液體介質(zhì)中的貫徹距離、速度衰減規(guī)律以及高速彈體引起的水錘效應(yīng),是廣大學(xué)者長期關(guān)注的問題,對此已開展了大量的研究。Lee 等[2]在經(jīng)典沖擊壓力模型的基礎(chǔ)上,將彈體簡化為點(diǎn)源,對高速彈體在液體中引起的空腔效應(yīng)展開了研究,推導(dǎo)建立了空腔動力學(xué)方程。李營等[3]在考慮長徑比的影響下,將阻力因數(shù)作為與初速度有關(guān)的常量來處理,擬合得到了阻力因數(shù)計算經(jīng)驗公式。沈曉樂等[4]考慮彈體的墩粗變形進(jìn)行了實驗研究,對方形彈體入水的阻力因數(shù)進(jìn)行了修正。郭子濤[5]分析了彈頭形狀對彈體水彈道穩(wěn)定性以及速度衰減規(guī)律的影響,提出了一種簡化的彈體頭型阻力因數(shù)估算公式。Zhao 等[6]在考慮雷諾數(shù)的基礎(chǔ)上對球形彈入水阻力因數(shù)進(jìn)行擬合,提出了球形彈入水阻力因數(shù)經(jīng)驗公式。Zhang 等[7]在Lee 等[2]研究工作的基礎(chǔ)上針對長徑比在0.1~0.5 之間的餅形彈進(jìn)行研究,并對餅形彈入水阻力因數(shù)計算公式進(jìn)行了修正。這些研究普遍簡化了阻力因數(shù)的影響因素,與實際物理過程有較大的區(qū)別。
本文中,擬采用理論分析和數(shù)值模擬的方法,針對以上存在的問題,開展截錐形彈體垂直穿透液體介質(zhì)時的速度衰減特性研究,重點(diǎn)分析不同彈體頭形因數(shù)對其速度衰減規(guī)律的影響,以期建立阻力因數(shù)與彈體的瞬時速度之間的關(guān)系,并引入彈體形狀參數(shù)的影響,建立截錐形彈體在液體中的速度衰減分析模型。
當(dāng)彈體在液體介質(zhì)中高速運(yùn)動時,彈體損失的動能主要轉(zhuǎn)化為液體的動能及空泡的壓力勢能,而空泡的存在極大地減小了彈體與液體的接觸面積,從而彈體運(yùn)動所受的摩擦阻力減小,因此壓差阻力是影響彈體阻力的主要因素。則根據(jù)牛頓第二定律建立彈體運(yùn)動方程如下[8]:
本節(jié)設(shè)計4 種不同頭形因數(shù)的截錐形彈體,采用顯式動力學(xué)非線性有限元程序Autodyn 建立彈體的仿真模型并開展系列工況的數(shù)值計算,分析頭形因數(shù)對彈體垂直入水后速度衰減的影響規(guī)律。
圖1 為截錐形彈體示意圖,定義彈體頭部圓臺直徑D2與主體最大直徑D1的比值D2/D1為截錐形彈體頭形因數(shù) ψ 。各工況中彈體主體長度L1=30 mm ,錐體長度L2=10 mm ,主體直徑D1=12 mm 保持不變,4 種彈體頭部圓臺直徑分別為0、4、8、12 mm,對應(yīng)頭形因數(shù) ψ1~ψ4分別為0、1/3、2/3、1。為了控制變量的數(shù)量,不同頭形因數(shù)截錐形彈體的質(zhì)量均調(diào)整為32.5 g,有限元模型如圖2 所示。本次數(shù)值模擬采用二維軸對稱建模方法,水域網(wǎng)格尺寸為400 mm×400 mm,采用邊長為1 mm 的四邊形歐拉單元進(jìn)行離散,四周設(shè)置透射邊界條件,模擬無限水域,彈體用拉格朗日網(wǎng)格離散。
數(shù)值模擬中,采用S h o c k 狀態(tài)方程和Johnson-Cook 本構(gòu)關(guān)系描述彈體材料的動態(tài)力學(xué)行為[12]:
圖 1 截錐形彈體示意圖Fig. 1 Schematic diagram of a truncated cone-shaped projectile
式中: σy為材料的動態(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)方程各項參數(shù)[14]Table 1 Parameters of equation of state for water[14]
本文中對每種不同頭形因數(shù)的截錐形彈體各設(shè)置了5 種不同的入水初速度開展數(shù)值模擬,共包含20 種工況條件,如表2 所示。
表 2 數(shù)值計算工況Table 2 Numerical calculation conditions
采用本文的數(shù)值模擬方法對試驗過程[5]開展計算,通過對比分析數(shù)值計算結(jié)果和試驗數(shù)據(jù)來驗證數(shù)值方法對高速彈體入水過程分析的可靠性和準(zhǔn)確性。文獻(xiàn)[5]中使用輕氣炮發(fā)射圓柱形平頭彈體,該彈體直徑D=12.65 mm,長度L=25.4 mm,入水初速度為603 m/s,試驗中使用高速攝影儀記錄了彈體運(yùn)動軌跡及空泡變化,測量了不同時刻的空泡尺寸以及彈體速度等數(shù)據(jù)。將數(shù)值模擬得到的高速彈體入水過程中與試驗[5]中相應(yīng)的數(shù)據(jù)進(jìn)行對比,彈體入水后空泡尺寸對比如圖3 所示,可以發(fā)現(xiàn)在t=0.222 ms和t=0.440 ms 時刻相同位置處的空泡直徑數(shù)值計算結(jié)果與試驗測量數(shù)據(jù)吻合較好
圖 3 針對彈體直徑為12.65 mm,長度為25.4 mm,入水初速度為603 m/s 時,試驗與數(shù)值計算的空泡尺寸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)一步開展尺寸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ù)值計算,將彈體位移隨時間的變化以及彈體速度衰減的數(shù)值計算結(jié)果與相應(yīng)的試驗測量數(shù)據(jù)進(jìn)行對比,如圖4~5 所示。可以看出,本文數(shù)值模擬方法在計算彈體的運(yùn)動特性方面具有較好的精度,可為本文的進(jìn)一步的研究工作開展提供了準(zhǔn)確可靠的分析手段。
圖 4 長度為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 長度為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
采用驗證后的數(shù)值模擬方法開展表2 中20 種工況的數(shù)值計算,將頭形因數(shù)不同的截錐形彈體在不同入水速度下的速度衰減與入水距離的關(guān)系繪制于圖6,通過對比分析圖中數(shù)據(jù)可以發(fā)現(xiàn):
(1)當(dāng)截錐形彈體的頭形因數(shù)不變時,彈體入水后速度變化曲線的斜率隨著入水初始速度的升高而逐漸增大,說明入水速度越高,速度衰減越快;
(2)彈體入水初速度不變、頭形因數(shù)增大時,彈體剩余速度并非單純隨頭形因數(shù)增大而降低,這一特性在v0=1 400 m/s 時表現(xiàn)最明顯。4 種頭形因數(shù)的彈體均以1 400 m/s 的初速度垂直入水時,運(yùn)動400 mm后彈體剩余速度與初速度的比值vb/v0分別為0.699 0、0.771 0、0.709 9、0.512 7。說明入水初速度不變時,隨著彈體頭形因數(shù)的增大,vb/v0呈現(xiàn)先增大后減小的趨勢,在 ψ =1/3 時達(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
根據(jù)式(4)可以計算得出截錐形彈體垂直入水后的瞬時阻力因數(shù)Cd,并將其作為縱坐標(biāo),將彈體的瞬時速度vb與入水初速度v0的無量綱比值作為橫坐標(biāo),繪制Cd- (vb/v0) 關(guān)系圖。由于篇幅有限,本文中僅繪出ψ=0 的彈體Cd- (vb/v0) 關(guān)系圖,如圖7 所示。從圖7 可以發(fā)現(xiàn),阻力因數(shù)Cd隨著彈體速度的衰減呈振蕩趨勢,且當(dāng)截錐形彈體速度較低時,隨著侵徹過程的進(jìn)行,速度衰減率出現(xiàn)較小幅度的上升。究其原因是因為隨著速度的降低,彈體誘導(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
從完整的侵徹過程來看,阻力因數(shù)Cd變化過程較復(fù)雜,阻力因數(shù)隨著彈體速度的降低,呈現(xiàn)如上述圖像所示的震蕩變化趨勢。在以往的研究方法中將其作為常量來處理,不具備廣泛的適用性。從簡化計算的角度出發(fā),本文中采用線性關(guān)系進(jìn)行擬合:
式中:a0和a1為未知參數(shù)。
針對圖7 中的數(shù)據(jù),采用最小二乘法進(jìn)行擬合,擬合的結(jié)果如圖7 中紅線所示。同時,對頭形因數(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,用最小二乘法對參數(shù)a0和a1進(jìn)行擬合:
式中:a01、a02、a03、a04、a05、a11、a12、a13、a14、a15為待定參數(shù),結(jié)果如圖8 所示。
由擬合結(jié)果來看,四階方程可以較好的擬合a0、a1關(guān)于無量綱數(shù)M的關(guān)系,可以發(fā)現(xiàn)同一頭形因數(shù)下參數(shù)a0及a1的高度相關(guān),但變化趨勢相反,近似關(guān)于某一水平線y=n1對稱。隨著頭形因數(shù)的增大,參數(shù)a0及a1的變化逐漸趨于復(fù)雜,這是因為當(dāng)頭形因數(shù)增大時,彈體的垂直迎流面積增大,侵徹時鐓粗變形效應(yīng)變得明顯,導(dǎo)致阻力因數(shù)的變化趨于復(fù)雜。參數(shù)a0、a1的擬合結(jié)果列于表4。
圖 8 對于頭形因數(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ù)Cd:
通過2.2 節(jié)中已驗證的計算方法開展系列模擬計算,通過對公式(18)的編程求解,得到彈體的速度衰減。將數(shù)值計算結(jié)果與上述速度衰減模型計算結(jié)果進(jìn)行對比,驗證其可靠性。
驗證計算中,設(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)行驗證計算,結(jié)果如圖10~11 所示。
圖 10 入水速度為400 m/s,頭形因數(shù)不同的彈丸,速度隨入水距離衰減的模擬結(jié)果與公式計算結(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
由計算可以發(fā)現(xiàn),對于彈徑相同、入水初速相同而頭形因數(shù)不同的彈體,其速度衰減規(guī)律有著顯著差異。本文的計算分析模型中考慮了因頭形因數(shù)不同而引起的阻力因數(shù)變化對彈速衰減的影響,與數(shù)值計算結(jié)果吻合較好。當(dāng)彈體入水初速為400 m/s,頭形因數(shù)為1/2 時,數(shù)值計算數(shù)據(jù)與本文公式計算結(jié)果之間的誤差約為6.8%,其他兩種頭形因數(shù)工況下誤差更小;在彈體頭形因數(shù)為1/6、不同速度工況中,公式計算與數(shù)值模擬結(jié)果之間的最大誤差約為8.71%。
圖 11 頭形因數(shù)為1/6,入水速度不同的彈丸,速度隨入水距離衰減的模擬結(jié)果與公式計算結(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
從上述對比結(jié)果可以發(fā)現(xiàn),本文的速度衰減模型可對不同頭形因數(shù)的截錐形彈體垂直侵徹液體介質(zhì)過程中的速度衰減進(jìn)行快速預(yù)報,彌補(bǔ)了以往研究中將阻力因數(shù)Cd當(dāng)作常量計算的不足,更符合彈體垂直侵徹液體介質(zhì)的實際物理過程。
本文中提出了與彈體形狀和無量綱速度相關(guān)的阻力因數(shù)的計算方法,通過建立運(yùn)動方程得到了截錐形彈體垂直侵徹液體介質(zhì)速度衰減的計算分析模型。并通過算例分析,驗證了本文計算方法的可靠性。得到以下結(jié)論:
(1)高速彈體在液體介質(zhì)中運(yùn)動方程的阻力因數(shù)與彈體形狀和無量綱速度有關(guān);該因數(shù)對彈體高速運(yùn)動過程中與液體的接觸面積進(jìn)行了修正,將彈體頭形因數(shù)和速度因素考慮在內(nèi),更接近實際的物理過程。
(2)將阻力因數(shù)引入運(yùn)動方程后得到的計算結(jié)果與數(shù)值模擬數(shù)據(jù)吻合較好, 說明本文提出的分析模型可對不同頭形因數(shù)的截錐形彈體垂直侵徹液體介質(zhì)過程中的速度衰減進(jìn)行快速預(yù)報。
(3)由于彈體侵徹液體時的作用機(jī)理非常復(fù)雜,且彈體高速侵徹時會發(fā)生彈體鐓粗和變形,同時流體具有可壓縮性,對彈速衰減影響較大。本文的分析中未考慮其影響,在后續(xù)工作中將進(jìn)一步的深入研究。