基于火藥燃氣壓力空間變化的火炮發射動力學研究

于情波， 楊國來， 葛建立， 孫全兆， 蕭 輝

(南京理工大學 機械工程學院， 南京 210094)

火炮射擊過程中彈丸與身管內壁的接觸碰撞作用以及高壓火藥燃氣對身管內壁的沖刷作用，使得身管產生一定的彈性振動，彈丸出炮口時身管炮口處的振動對彈丸初始運動姿態具有直接的影響[1-2]，從而影響火炮射擊精度。對于身管振動的研究，姜沐等[3]將身管簡化為一懸臂梁，彈丸的作用簡化為一加速移動的載荷，建立了身管的運動方程，求得了級數形式的解析解。劉寧等[4]把身管簡化成等截面懸臂梁，根據Bernoulli-Euler初等梁理論建立了身管振動方程，用模態分析法求解了身管振動方程。蘇忠亭等[5]將身管等效為變截面懸臂梁，并計算分析了內彈道時期彈丸移動對身管振動特性的影響。文獻[6]同樣研究了移動質量作用下懸臂梁的振動特性。該研究方法對身管實際結構過于簡化，并且忽略了火藥燃氣的沖擊作用,計算得到的解析解與實際情況誤差較大。

有限元分析方法作為火炮發射動力學研究的常用方法，能夠反映火炮的模態特征，應力，應變的分布情況，可以模擬彈丸在內彈道時期的運動形態以及火炮零部件的動力學響應。葛建立等[7]考慮了彈丸彈帶、前定心部與膛壁的碰撞及膛線對彈丸的扭轉作用，分別仿真計算了火炮在有無膛線的情況下彈丸在膛內的運動過程。許耀峰等[8]建立某大口徑火炮線膛身管與制導炮彈耦合動力學有限元模型，數值計算了膛線形式，膛線深度對彈炮動力學響應的影響。吳會民等[9]采用數值方法對彈丸身管耦合作用下系統動力響應進行了計算，得到了身管扭轉動力響應和橫向動力響應的基本規律。在現有的大部分針對彈丸在膛內的運動響應及身管振動特性的研究中，僅限于彈丸膛內響應規律研究，忽略了直接作用于身管內壁的火藥燃氣壓力對身管動態響應的影響。火藥燃氣壓力屬于強沖擊載荷，其壓力值和作用于身管內壁的分布情況均與時間有關，并且具有階躍型脈沖動態效果，是激發身管振動不可忽視的因素。可學為等[10]研究了火藥燃氣壓力對身管振動的激勵作用，利用 ABAQUS 軟件建立身管的三維有限元模型，將火藥燃氣壓力作用區域沿身管軸向離散成一系列較小的作用區域，并將對應的壓力曲線施加到身管各個離散區域。火藥燃氣壓力根據彈丸的運動位置實時作用于彈后身管內壁，利用此方法對身管施加燃氣壓力，不能準確模擬實際壓力場二維分布情況，并且沒考慮彈丸與身管的耦合作用，求得的結果與實際情況存在較大差別。

另外，關于身管強度失效的報道屢見不鮮，但現有關于身管強度分析的文獻少之又少。數值模擬實現火炮身管實際受載環境，得到身管強度動態變化情況，對身管強度失效分析具有重要的參考價值。部分文獻對身管內膛損傷機理進行了研究分析，曾志銀等[11]研究了彈丸裝填不到位造成的沖擊加載對身管材料動態性能的效應，分析了導致陽線起始段局部雙側棱邊斷裂損傷的可能原因。劉廣生等[12]對某型火炮多發射擊工況下內膛損傷破壞過程進行了數值模擬計算，分析了彈帶擠進內膛過程中身管內壁材料性能隨射彈發數變化的規律。文獻中數值計算了身管內膛在火炮射擊中的動力學響應，但均忽略了火藥燃氣壓力的影響，不能真實描述身管實際動力學響應。曾志銀等[13]運用 ANSYS 有限元分析軟件的APDL語言，計算分析身管在徑向受載下的動力學響應，通過載荷上升時間確定身管截面承受的階躍型脈沖，并沒有嚴格按照實時彈丸軸向運動位置定義燃氣壓力二維分布，不能準確模擬實際射擊環境下身管受載情況。

考慮到高壓火藥燃氣壓力二維分布特性對身管動力學響應的影響，本文在現有的火炮發射動力學研究的基礎上，建立了彈丸與身管耦合非線性動力學模型。借助有限元軟件提供的用戶自定義子程序，實現了彈丸膛內運動與經典內彈道火藥燃燒相互關聯的動力學過程，基于身管振動理論，數值計算過程計及了火藥燃氣壓力對身管振動的激勵作用，實現了膛內火藥燃氣壓力隨彈丸運動而動態變化的二維分布情況。結合實彈射擊實驗，有效驗證了建立的數值模型的正確性，對火炮發射過程的數值模擬提供了更接近實際的力學邊界條件，更真實地再現了身管動力學響應。

1 力學模型

借鑒歐拉-伯努利梁模型[14]，身管可以簡化為一端固定另一端自由的變截面懸臂梁。身管在整個內彈道過程中受到火藥燃氣膨脹波以及與彈丸接觸碰撞的作用，簡易力學模型如圖1所示。圖中M為彈丸的質量，Fcont為不計慣性效應時彈丸膛內運動作用于身管上的等效合外力，vt為不同時刻對應的彈丸運動速度。彈后空間為火藥燃氣作用區域，火藥燃氣壓力作用于身管內壁的荷載條件隨彈丸的運動而動態變化，屬于時變動力學體系。

根據達朗貝爾原理[15]，計及阻尼和彈丸慣性效應，移動彈丸以及火藥燃氣壓力作用下身管徑向運動方程為

(1)

式中：EI(x)為身管的抗彎剛度；H為身管速度衰減系數；y(x,t)為身管x位置處t時刻的撓度；F為身管受載情況，可以近似為

F=δ(x-s(t))Fcont+∑ξ(x-s(t))P(x,t)ds

(2)

(3)

(4)

借助有限元分析軟件建立物理離散模型，并利用中心差分法對相關力學模型進行計算分析。

2 彈炮耦合動力學有限元模型

本文以某122 mm口徑榴彈炮為研究對象，建立基于彈炮耦合的火炮發射動力學模型，并計及火藥燃氣壓力二維分布對身管動態響應的影響。

2.1 有限元模型

對火炮主要零部件進行離散建模，主要包括上架，搖架，后坐部分與彈丸等，對于不重要的附屬機構以集中質量單元模擬，通過剛性單元與相鄰節點連接模擬質量分布對結構強度的影響。上架與搖架結構主要采用等參四邊形單元和等參六面體單元，平衡機簡化為彈簧單元，高低機齒輪齒弧表面的接觸采用剛性單元連接進行處理。身管膛線為混合膛線，用分段掃描拉伸方式沿膛線空間曲線掃描獲得身管網格模型，身管坡膛處局部網格圖如圖2所示。彈丸彈體的變形不作為主體研究對象，故將彈體做剛性化處理，彈丸彈帶處理為彈塑性，在彈帶表面創建與身管膛線形狀匹配的溝槽，以彈帶完全擠進膛線時刻為數值計算起始時刻，預設溝槽彈帶局部網格圖如圖3所示。

圖2 線膛身管局部網格圖Fig.2 Local view of rifled barrel mesh

圖3 預設溝槽彈帶局部網格圖Fig.3 Local view of band mesh with predefined groove

有限元模型中彈帶與身管內壁之間定義為點-面接觸，彈體前后定心部與身管內壁以及后坐部分與搖架導向部定義為面-面接觸，均采用罰函數法進行處理。炮尾與反后坐裝置之間則采用自由度耦合的方式模擬后坐部分在炮膛軸線方向的滑移運動，僅釋放炮膛軸線方向的平移自由度，其他自由度均約束不動。搖架耳軸與上架耳軸室之間同樣采用自由度耦合的方式定義，僅釋放耳軸室中心軸線方向的轉動自由度。

2.2 材料模型

內彈道過程中復雜的動態物理化學場使得彈帶材料常常處在大應變、高溫和大應變速率等綜合因素下發生彈塑性應變。采用Johnson-Cook模型描述彈帶材料塑性階段力學行為[16]

(5)

(6)

式中：Tr為參考溫度，一般取室溫；Tm為常態下材料的熔化溫度。

身管采用炮鋼材料PCrNi3MoVA，彈帶材料為H96黃銅，對應的材料參數設置如表1所示。

表1 數值計算涉及的材料相關參數Tab.1 Material parameters value for numerical computation

2.3 子程序實現內彈道物理場變化過程

火炮內彈道時期伴隨復雜的物理場變化過程，主要體現在火藥燃氣壓力變化以及對身管作用區域的變化，即二維分布特性。火藥燃氣一方面沿身管軸向流動，推動彈丸以及后坐部分完成內彈道運動過程。另一方面，火藥燃氣壓力作用于彈后身管內壁，壓力值與作用區域隨著彈丸的運動而動態變化，可以根據實時測得的彈丸位置實現燃氣壓力二維分布情況。利用VUAMP以及VDLOAD子程序接口可以實現內彈道物理場變化過程。

2.3.1 推力子程序實現火藥燃燒過程

多數學者采用彈底加載實測內膛壓力曲線的方法實現彈丸主動力的加載，采用這種方法相當于計及了兩次次要功，其彈丸出膛參數的準確程度不高[17]，并且該壓力曲線并不能表征數值模擬對應的實際內彈道過程。應用推力子程序(VUAMP)可以實現內彈道火藥燃燒規律與彈丸膛內運動的關聯效應。針對整個內彈道耦合過程做如下假設：

(1) 忽略內彈道起始時刻對應的彈丸卡膛以及彈帶擠進過程，假定擠進完成時膛內平均壓力為30 MPa，設為膛壓初始值。

(2) 以炮膛底面中心點為坐標原點，身管軸線指向炮口方向為x軸正向，基于拉格朗日模型假設得到的氣體壓力分布作為彈后空間身管內壁的載荷分布，彈后空間的氣體壓力呈拋物線分布，如下所示

(7)

式中：st為彈丸運動位移；φ為二次功系數；Pt為火藥氣體平均壓力。

由上式可知，彈底壓力寫為平均壓力的函數為

(8)

膛底壓力寫為平均壓力的函數為

(9)

利用Fortran語言對火藥燃燒過程編譯，在每個計算增量步，數值計算模型實時傳輸彈丸運動參量，子程序計算得到相應的實時內膛平均壓力Pt，并傳入有限元模型。彈底壓力Pd由式(9)計算得到，推動彈丸完成膛內運動過程。膛底壓力由式(10)計算得到，并作為后坐運動主動載荷作用于膛底表面。

本計算模型采用的是模塊化裝藥，相應的經典內彈道方程組可表示如下

式中：z為藥粒相對燃燒厚度；ψ為相對已燃體積；χ，λ，μ為與藥粒形狀相關的內彈道參數；w為裝藥量；f為火藥力；l0，l(t)為藥室縮頸長、彈丸行程；u1，e1，n為火藥燃燒系數，火藥弧厚，燃燒指數；θ為熱力指數；S為炮膛斷面面積。

2.3.2 荷載條件子程序實現身管內壁受載過程

膛內燃氣壓力作用于彈后身管內壁，受載區域隨彈丸運動而動態變化，利用荷載條件子程序(VDLOAD)實現身管實際受載情況。以時間增量步為單位，主程序對子程序進行調用運算，子程序主要實現步驟如下：

步驟1主程序向子程序傳輸實時彈丸運動位置l(t)，實時火藥燃氣壓力Pt。

步驟2以身管內表面所有積分點為加載單元，子程序依據主程序獲取所有積分點的軸向坐標。

步驟3子程序進行相關判斷分析以及賦值燃氣壓力。依次對每個積分點進行賦值，在滿足對應積分點的軸向坐標小于彈丸運動位移的情況下，將對應位置的實時火藥燃氣壓力賦予該積分點。

步驟4所有積分點賦值結束，返回主程序。

在每次子程序調用過程中，整個子程序實現流程示意圖如圖4所示。

圖4 子程序實現流程示意圖Fig.4 The flow diagram of subroutine operation

3 數值仿真結果分析以及實驗對比

3.1 再現身管應力響應動態過程

圖5顯示了不同時刻對應的身管截剖面動力學響應情況，圖中隱藏彈帶部件是為了清晰觀察身管的應力分布。從圖中可清晰看到不同時刻彈丸在膛內的運動位置，身管主要應力分布在彈后區域，身管的受載面積隨彈丸的運動而動態變化，與實際發射環境相一致，實現了火藥燃氣壓力作用于身管內壁的實際荷載條件。

圖5 身管軸向截剖面應力分布云圖動態變化過程Fig.5 Dynamic change process of stress distribution in axial cross section of barrel

3.2 實彈射擊實驗簡述

為了驗證所建立的非線性動力學模型的正確性，進行了該型號火炮實彈射擊試驗，將火炮上架及以上部分固結于發射臺。采用七維航測科技公司的SDI-ARG-720型角陀螺傳感器，用于測量炮口高低方向的運動角參量，數據采集系統為DEWETRON 1201系統，試驗中角陀螺儀在身管上的布置位置如圖6所示。

采用非接觸式測試系統獲取炮口振動位移參量，在炮口外表面布置和標定專用的測試標識，運用高速攝影系統記錄到測試標識的運動軌跡，并通過專用的分析系統處理進而得到相應的振動參數。高速攝影設備采用IDT公司生產的Y3-S2高速攝影機，光學鏡頭為Nikkor AF-S 400mm F2.8D ED鏡頭，采用Xcitex公司開發的ProAnalyst軟件對圖像進行跟蹤處理。高速攝影系統及其支架固定在炮口側面的地面上，其布置示意圖如圖7所示。射擊工況與數值模擬工況相同，為0°射角。

圖6 角陀螺儀裝配實物圖Fig.6 The assemble paragraph of angle gyroscope

圖7 高速攝影系統布置示意圖Fig.7 The arrangement diagram of high-speed camera system

3.3 炮口振動參量對比分析

數值模擬火炮射擊過程，在與實驗對應的炮口振動實測位置設置細小的桿單元，并通過耦合建模與身管進行連接，得到炮口振動量的時程變化曲線，并與試驗測得的數據進行對比。圖8和圖9為通過實彈射擊試驗獲取的炮口測試點振動垂向位移和角位移與仿真數據的對比。從圖中可以看出，仿真得到的炮口測試點的位移和角位移時程變化情況與試驗測得的數據吻合較好。本文以彈丸彈帶脫離時刻為彈丸出炮口的度量標準，彈丸出炮口時刻對應的炮口振動量結果對比如表2所示，誤差均在可接受范圍之內，有效驗證了本文計及身管實際受載火炮發射動力學模型的正確性。

圖8 炮口垂向線位移Fig.8 The muzzle displacement in vertical direction

圖9 炮口垂向角位移Fig.9 The muzzle angular displacement in vertical direction

表2 炮口初始擾動測試與仿真結果對比Tab.2 Comparison of muzzle initial vibration by test and simulation

3.4 身管動力學響應

現有的關于火炮身管振動以及身管強度數值分析的文獻中，由于燃氣壓力徑向效應加載不準甚至無法加載的問題，大部分僅考慮了燃氣壓力的軸向效應，而忽略了燃氣壓力二維分布特性對身管動態響應的影響。為有效表述燃氣壓力徑向效應對身管動力學響應的影響，預先設置了兩種數值模擬環境：計算工況一，僅考慮高壓火藥氣體軸向效應，不考慮身管徑向受載；計算工況二，綜合彈炮耦合與身管徑向受載，計及火藥燃氣壓力二維分布特性，由子程序根據彈丸實時運動位置定義載荷大小與分布。

3.4.1 身管振動響應分析

圖10與圖11分別顯示了兩種數值模擬環境下炮口在垂向以及水平方向的振動情況。綜合垂向以及水平方向炮口振動曲線可知，計及燃氣壓力作用下，身管振動較劇烈。高壓氣體沖擊線膛身管內表面，膛線的存在使得身管整體受載不對稱，在水平面內的振動情況呈現一定的波動，且振動幅度明顯強于僅考慮燃氣壓力軸向效應的情況。由得到的炮口振動曲線可知，膛內燃氣壓力對身管振動具有直接影響，因此在數值模擬研究身管振動問題時，準確描述火藥燃氣壓力二維分布特性是必要的。

圖10 炮口垂向線位移Fig.10 The muzzle displacement in vertical direction

圖11 炮口水平方向線位移Fig.11 The muzzle displacement in horizontal direction

3.4.2 身管動態強度分析

為捕捉身管應力響應隨時間的變化情況，在身管內壁沿軸向等間距選取四個積分單元，A積分點靠近藥室部，并依次向炮口方向選取B，C，D三個積分點，得到對應的應力時程曲線。圖12對應僅考慮彈炮耦合作用的情況，身管主要受到彈丸與身管相互接觸碰撞產生的載荷，身管局部指定區域應力響應短暫。運動的初始時刻，彈丸慣性效應不明顯，與身管接觸碰撞產生的載荷較小，與彈丸接觸區域應力幅值變化不明顯。隨著運動的進行，身管應力幅值曲線的峰值與彈丸運動速度具有正相關。圖13對應計及火藥燃氣作用于身管內壁的情況，數值計算得到的身管指定位置應力幅值曲線反映了火藥氣體實際作用情況，與不計及火藥燃氣作用于身管內壁的情況相比，更能反映出身管實際動力學響應。

圖12 僅考慮燃氣壓力軸向效應下身管應力時程變化曲線Fig.12 Stress-time curves of barrel under axial effect of gases pressure

傳統身管設計依據靜強度設計理論，通過最大徑向壓力確定身管各部位的安全系數，而身管徑向受載屬于與二維分布有關的動態過程,文獻[13]指出身管強度設計按動態問題處理才是符合實際的設計。圖14顯示了兩種數值模擬環境下得到的身管應力峰值時程分布曲線。僅考慮彈炮耦合作用下，身管應力響應隨彈丸速度的增加呈現遞增的趨勢，取決于彈丸的慣性力。計及火藥燃氣壓力與彈炮耦合作用下，身管應力響應較劇烈，在身管強度分析方面更具說服力。整個動力學過程身管等效應力峰值曲線最大值為508 MPa，位于彈丸導向部，對應火藥燃氣壓力最大值時刻。由Von Mises屈服準則可知，身管材料的屈服應力為1 040 MPa，安全系數大于2，符合身管強度要求，高壓火藥燃氣以及彈炮耦合作用下身管處于較安全狀態。

圖13 計及徑向受載身管應力時程變化曲線Fig.13 Stress-time curves of barrel considering gas pressure

圖14 身管應力峰值時程分布曲線Fig.14 Stress peak distribution curves of barrel in time domain

4 結 論

本文在現有的基于彈炮耦合火炮發射動力學模型的基礎上，考慮到高壓火藥燃氣的沖擊作用對身管動力學響應的影響，借助有限元軟件提供的用戶自定義子程序，實現了彈丸膛內運動與經典內彈道火藥燃燒相互關聯過程，同時實現了膛內火藥燃氣壓力二維分布特性，解決了以往燃氣壓力徑向效應加載不準甚至無法加載的問題。結合實彈射擊實驗，驗證說明了所建立的數值模型是可取的，為火炮發射過程的數值模擬提供了更接近實際的物理場邊界條件。通過對比分析身管在有無燃氣壓力徑向效應數值模擬環境下的動力學響應可得出如下結論：

(1) 由于燃氣壓力對線膛身管的沖擊效應，身管在計及燃氣壓力二維分布特性條件下的振動響應明顯強于僅考慮彈炮耦合的情況，說明在數值模擬研究身管振動問題時，準確描述火藥燃氣二維分布特性是必要的。

(2) 計及燃氣壓力二維分布的火炮發射動力學數值模擬可以得到接近實際物理場作用下的結構強度響應，為今后火炮設計以及故障分析提供了更具說服力的參考依據。