小口徑火炮彈炮耦合動態(tài)響應(yīng)有限元時程分析

蘇忠亭，徐 達，李曉偉，韓振飛

(1.裝甲兵工程學院 兵器工程系，北京 100072;2.北京京航計算通訊研究所，北京 100074)

火炮身管的柔性對射擊精度影響很大，彈丸在膛內(nèi)運動時，膛線賦予彈丸旋轉(zhuǎn)運動，以保證彈丸飛行穩(wěn)定。彈丸膛內(nèi)運動時，彈帶和前定心部與炮膛內(nèi)壁的接觸碰撞會導(dǎo)致身管激烈振顫，炮口振動引起彈丸跳角，直接影響射擊命中率。文獻[1]建立了彈丸膛內(nèi)運動引起炮管振動的方程，沒有考慮彈丸與身管的碰撞－接觸;文獻[2]基于多體接觸理論，考慮了彈丸彈帶、前定心部與膛壁的碰摩，但其有限元法本質(zhì)上仍將彈性身管簡單轉(zhuǎn)化為剛體接觸;文獻[3]采用等參六面體縮減積分單元，考慮了彈丸彈帶和前定心部與膛壁的碰撞及膛線對彈丸的扭轉(zhuǎn)作用，對彈丸在膛內(nèi)的運動過程進行了仿真計算，但膛線身管結(jié)構(gòu)復(fù)雜，很難建立準確的三維模型;文獻[4]研究了彈管間隙和加速移動質(zhì)量對身管振動的影響，采用懸臂梁的振型作為身管的近似模態(tài)，無法得到身管的真實模態(tài)。由于小口徑火炮身管短、膛壁薄，適合采用有限單元法編程求解。本文利用Dubowsky模型，基于有限單元法，將彈帶定為彈塑性體，前定心部視為剛體，身管為柔體，分別研究了彈帶凹槽部、凸起部和前定心部對身管的瞬態(tài)沖擊作用，根據(jù)彈帶行進過程確定身管各單元在各時間步長內(nèi)所受的載荷譜，最終建立了彈丸膛內(nèi)運動引起的身管動態(tài)響應(yīng)模型。

1 身管有限元模型

1.1 坐標系的建立

根據(jù)結(jié)構(gòu)動力學可知，身管屬于內(nèi)中空懸臂梁的空間軸對稱問題，采用圓柱坐標系比較方便。建立坐標系ozρθ，o為以炮膛軸線與膛線起始部端面的交點，oz軸為未變形的身管軸線，oρ軸為身管徑向，以o垂直向上為θ角零點，向右為起始方向，如圖1所示。

圖1 圓柱坐標系下的身管結(jié)構(gòu)Fig.1 Tube structure in cylindrical coordinate system

身管有n條膛線，膛線在oρθ平面內(nèi)展開曲線為f(z);身管內(nèi)壁半徑為 r1，外壁半徑為r2，膛線纏度為η，深度為d，陽線寬為a，彈丸彈帶寬為b，前定心部距彈帶前端距離為c，射角為α;在t時刻，彈丸質(zhì)量M在膛內(nèi)軸向運動速度為v(t)，沿膛線旋轉(zhuǎn)速度為ω(t)，彈丸膛內(nèi)運動距離為s(t)=，均由內(nèi)彈道計算得出。

圖2 有限單元劃分Fig.2 Division of finite element

1.2 有限單元劃分

身管單元包括兩部分，陽線部分和其余部分。如圖2(a)所示。首先取陽線進行網(wǎng)格劃分，軸向以彈帶寬b為單元長度的六面體網(wǎng)格，如圖2(b)所示。在ozρθ坐標系中，將陽線去掉，對身管進行離散，取軸向長度b，夾角θ，徑向dr的六面體網(wǎng)格單元，如圖2(c)所示。兩種單元通過公用面的節(jié)點進行載荷傳遞。它分為兩種情況:陰線內(nèi)壁的六面體網(wǎng)格單元直接承受彈帶凸起部載荷作用;與陽線單元直接連接及不與彈帶接觸的身管內(nèi)部結(jié)構(gòu)單元，其節(jié)點載荷為由陽線節(jié)點傳遞的載荷或陰線內(nèi)壁六面體網(wǎng)格單元節(jié)點傳遞的載荷。該六面體網(wǎng)格單元在圓柱坐標系中與三維六面體等參單元的局部坐標系中規(guī)則的母單元等效，采用八節(jié)點L型插值函數(shù)N=[N1N2N3N4N5N6N7N8]作為單元位移函數(shù)，單元節(jié)點位移向量為={u1，…，u8}，其中ui為三個元素的列向量。則單元撓度為:

1.3 單元編碼

對于變速彈丸依托膛線旋轉(zhuǎn)前進的身管結(jié)構(gòu)，單元的編碼非常重要。結(jié)合圓柱坐標系和時間軸，建立四維向量為單元編碼。設(shè)某單元編碼為 ei，j，k，t，其中，i，j，k分別代表整數(shù)表示的周向、徑向和軸向的遞進關(guān)系，t為實數(shù)，表征在t時刻與彈丸發(fā)生觸碰的單元。

對于i，由于彈丸彈帶的凹槽部和凸起部在內(nèi)膛圓周對各單元進行觸碰，首先對彈帶各部編號，i={0，1，2，…，2n－1}，以θ角零點為0，順時針遞進。彈丸沿膛線運動時，與彈帶各編號對應(yīng)的陽線和陰線在i上不變;

2 彈炮耦合過程中身管動態(tài)響應(yīng)模型

2.1 彈丸作用下身管的徑向動力學方程

將身管簡化為膛線起始部一端固定的變截面彈性懸臂梁[5]，根據(jù) Euler-Bernoulli梁理論，假設(shè)橫向剪切應(yīng)變和橫向線應(yīng)變?yōu)榱悖砉軓澢匠虨閇6]:

其中:f(z，t)為身管所受的接觸力，EI(z)為身管(梁)的抗彎剛度，h為身管速度衰減系數(shù)，由身管的結(jié)構(gòu)和材料性能決定。u(z，t)為身管z位置處t時刻的撓度，Me為身管在單元e的質(zhì)量分布。

初始條件為:

邊界條件為:

將彈丸看作移動質(zhì)量－彈簧－阻尼，彈簧和阻尼分別通過彈帶凸起部和前定心部與身管發(fā)生觸碰。剛度系數(shù)和阻尼系數(shù)根據(jù)接觸體的材料彈性模量和幾何形狀合適選取。如圖3所示，彈丸質(zhì)量M以速度v(t)在膛內(nèi)運動，通過彈簧K、阻尼C與身管(梁)相互作用，陰線面所受彈丸彈帶或前定心部載荷為以身管軸線即懸臂梁中性軸為起點的徑向力，大小為:

其中:δ為Dirac delta函數(shù);s(t)為t時刻彈丸膛內(nèi)運動距離，s(t)=v(t)dt，由內(nèi)彈道計算得出;ρ(t)為彈丸徑向位移為彈簧剛度系數(shù)為阻尼系數(shù)，塑性變形后的彈帶與內(nèi)膛的接觸剛度根據(jù)銅－鋼的材料彈性模量和幾何形狀合適選取，前定心部與內(nèi)膛的接觸剛度根據(jù)鋼－鋼的材料彈性模量和幾何形狀合適選取，阻尼系數(shù)采取工程上的經(jīng)驗值。

圖3 彈丸移動彈簧－阻尼－質(zhì)量系統(tǒng)作用下身管懸臂梁模型Fig.3 Tube cantilever beam model under the moving spring，damping＆ mass system of projectile

彈丸(彈帶或前定心部)徑向運動方程為:

整理得:

由于離散體系的平衡方程是在節(jié)點上進行分析的，因此對身管劃分的三維單元與平面微分方程之間主要依靠單元的等效節(jié)點載荷來銜接，根據(jù)靜力等效原則，將作用在身管(梁)上的非節(jié)點面力載荷轉(zhuǎn)化為等效節(jié)點載荷，三維單元的面力等效節(jié)點載荷為:

根據(jù)彈丸與身管相互作用的不同，將彈丸膛內(nèi)運動過程分為三個部分:階段一為彈丸彈帶由坡膛擠進膛線過程;階段二為直膛段彈丸彈帶和前定心部與身管共同作用過程;階段三為直膛段彈丸彈帶單獨作用過程。如圖4所示。

圖4 移動彈丸各階段膛內(nèi)運動過程Fig.4 Moving process in chamber of projectile

2.2 階段一身管內(nèi)壁單元受力分析

在階段一，彈丸彈帶前端面由0到2b。當彈帶完全嵌入膛線時，彈帶首先為彈性變形，逐步達到屈服極限后過渡到塑性變形，當?shù)竭_2b點時，彈帶塑性變形結(jié)束[7]。因本文主要對彈炮系統(tǒng)的徑向應(yīng)力應(yīng)變進行分析，彈帶擠進階段彈丸對膛線的沖擊主要體現(xiàn)在彈性變形階段[8]，在此階段，可采用彈性變形階段的彈簧－阻尼－質(zhì)量系統(tǒng)模擬彈丸與身管的相互作用。

將曲面內(nèi)壁簡化為平面，當彈丸運動到單元e時，聯(lián)立式(2)、(3)、(6)，此單元的外載荷向量為:

式中:Me為單元質(zhì)量。

式中:Δl為彈帶半徑強制量。

彈丸沒有直接作用的身管單元外載荷為自重:

2.3 階段二身管內(nèi)壁單元受力分析

在階段二，彈丸彈帶前端由2b到2ηr1－c。彈帶凸起部與陰線膛壁為彈性接觸，可采用彈簧－阻尼－質(zhì)量系統(tǒng)進行模擬，將膛線定義為主面，彈帶定義為從面，運用單純主從搜索算法確定接觸點位置;彈帶凹槽部與陽線的徑向受力為正壓力。彈帶在軸向受到彈底壓力作用、彈前空氣阻力作用及彈丸和身管單元的摩擦阻力，而本文的彈炮系統(tǒng)模型描述的是圓柱坐標系下身管徑向的受力及運動特征，故不考慮彈底壓力和摩擦阻力作用。由于彈管間隙的存在，彈丸前定心部與陽線產(chǎn)生彈性碰撞，之后分離與另一部分陽線產(chǎn)生碰撞，如此反復(fù)。

圖5 陽線受力分析Fig.5 Force analysis of rifling

陽線受力情況如圖5所示，陽線導(dǎo)轉(zhuǎn)側(cè)壓力使身管產(chǎn)生左旋回轉(zhuǎn)力矩，對徑向變形無影響，故不考慮。因彈帶材料采用理想彈塑性線性強化模型，在塑性變形后應(yīng)力應(yīng)變?nèi)猿市甭屎苄〉木€性關(guān)系。為簡化計算，取屈服應(yīng)力σs為彈帶凹槽部等效應(yīng)力。

故徑向正壓力為:

彈帶凸起部與陰線單元e的觸碰力模型與式(5)彈帶凸起部相同。

彈帶凹槽部與單元e的觸碰力模型為:

彈丸前定心部與單元e的觸碰力模型與式(5)相同，接觸分為兩部分:初始接觸采用單純主從搜索算法，膛線定義為主面，前定心部定義為從面，搜索初始接觸單元;前定心部與膛壁某單元碰撞后彈回，經(jīng)過兩個單側(cè)間隙B、軸向位移和彈丸旋轉(zhuǎn)角度，與下一單元碰撞[9]。設(shè)t時刻與彈丸前定心部碰撞的陽線單元為ei'，j，k'，t，則下一次與彈丸前定心部碰撞的陽線單元為ei'+1，j，k'+1，t，其中，

式中:B為前定心部與陽線單側(cè)間隙。

2.4 階段三受力分析

在階段三，彈帶與身管接觸模型與階段二相同。彈丸前定心部脫離炮口，與陽線不再接觸碰撞。故階段三彈炮耦合數(shù)值模型為式(7)與式(9)。

2.5 動力學系統(tǒng)建模及求解

為便于數(shù)值求解，綜合以上三個階段，編寫MCK矩陣[10]，建立系統(tǒng)動力學方程:

因矩陣[M]、[C]和[K]中有部分是時變的，需要以分塊矩陣形式表示，故:

式(15)中，分塊矩陣中的上標e即表示了在時間t參與接觸的時變單元，附加的剛度、阻尼只與移動荷載直接作用的單元位移有關(guān)。

彈丸膛內(nèi)運動屬于移動荷載問題，MATLAB軟件內(nèi)部函數(shù)中有專門針對移動荷載問題的ODE求解器[11]。在MCK矩陣基礎(chǔ)上，基于四階五次Ruge-Kutta微分方程進行二次開發(fā)，對彈炮耦合系統(tǒng)進行數(shù)值求解。

MATLAB提供的ODE45的調(diào)用格式為:

其中:odefun為指示微分方程等號右側(cè)的函數(shù)句柄，在迭代中計算MCK矩陣中的時變分塊。變量tspan為計算范圍，變量y0為系統(tǒng)初始狀態(tài)變量的值。options可以定義函數(shù)運行時的參數(shù)。T為時間點的列向量，Y為與時間T對應(yīng)的求解數(shù)組。

3 數(shù)值計算算例

取 θ=11.25°，M=84.2 kg，m=0.36 kg，r1=0.015 m，r2=0.025 m，d=0.000 3 m，a=0.002 m，b=0.001 5 m，c=0.08 m，η =30，B=0.000 5 m，E=300 GPa[12]，h=0.2[12]，dr=0.002 m，對于彈帶，對于前定心部，根據(jù)內(nèi)彈道計算得到的彈丸膛內(nèi)運動速度和旋轉(zhuǎn)速度，采用2.5節(jié)中數(shù)值求解方法，應(yīng)用MATLAB生成單元信息并進行彈炮耦合動態(tài)響應(yīng)計算，生成的單元個數(shù)為105605，節(jié)點個數(shù)為115 360，其中陽線單元為9 605。取炮口外壁周向θ=0°單元為研究對象，求解其徑向(即垂直向)位移。

實彈射擊時，選用技術(shù)狀態(tài)良好的武器系統(tǒng)進行炮口振動測試，高低、方位射角均為0°，彈種為曳光穿甲彈，單發(fā)射擊，在常規(guī)靶場進行。在炮口端正上方安裝3 000 g加速度傳感器，應(yīng)用NI公司高分辨率數(shù)據(jù)采集卡4472 B，在PXI－8105平臺上采用LabVIEW語言編寫測試框圖程序，測試與研究對象位置基本對應(yīng)的單元的垂直向振動位移，總共獲得了12組數(shù)據(jù)。12發(fā)彈的地面散布較小，12組振動數(shù)據(jù)的差異也非常小，由于彈丸出炮口時間較短(約6 ms)，該自動炮身管仍處于自由后坐階段，未壓縮緩沖器彈簧，從而與搖架、炮塔、車體等的相互作用并不明顯，對測試曲線中前7 s的幅值較低部分取均值并與仿真位移曲線對比如圖6所示。

由圖6可以看出，兩者非常接近，對兩者進行相關(guān)性分析，corrcoef=0.895，彈炮耦合數(shù)學模型較為準確。

分別對炮口內(nèi)壁周向 θ=0°、θ=90°、θ=180°、θ=270°的單元進行分析，得到彈丸膛內(nèi)運動時，身管的撓度曲線如圖7所示。

由圖7可以看出，彈丸繞膛線在膛內(nèi)運動一周，對內(nèi)壁各方向觸碰規(guī)律相似，但由于前定心部的碰撞，各方向又有所差別[13]。改動數(shù)學模型，不考慮階段二中前定心部的碰撞規(guī)律，取炮口內(nèi)壁周向θ=0°為研究對象，得出振動曲線如圖8所示。

圖6 仿真模型的實測曲線對比驗證Fig.6 Comparison of simulated model and tested model

圖7 炮口內(nèi)壁周向各單元撓度曲線Fig.7 Deflection curve of elements inside the muzzle

由圖8可以看出，不考慮前定心部的碰撞時，身管振動幅度和頻率都較小，可知前定心部的碰撞對彈炮耦合影響很大。

取炮口端面θ=0°處內(nèi)外壁單元為對象，時間歷程內(nèi)炮口振動曲線如圖9所示。

圖9 炮口內(nèi)壁與外壁單元振動曲線對比Fig.9 Comparison of elements vibration inside the muzzle and outside the muzzle

由圖9可以看出，與炮口內(nèi)壁相比，外壁單元的振動頻率和幅值都較低，表明與彈丸直接作用的炮口內(nèi)壁對彈丸影響更大。

4 結(jié)論

考慮彈丸彈帶擠進過程和剛體前定心部對內(nèi)膛的瞬態(tài)碰撞對身管振動的影響，將彈帶凹槽部視為剛塑性體，凸起部應(yīng)用彈簧－阻尼－質(zhì)量系統(tǒng)，身管簡化為柔性等截面懸臂梁，對從彈帶擠進膛線開始到彈帶出炮口結(jié)束的彈丸膛內(nèi)運動對身管陽線和陰線單元的觸碰作用進行了分析。由于計算量大，比較適合小口徑火炮。以某型火炮為例進行了數(shù)值計算，在實驗驗證的基礎(chǔ)上對炮口振動進行了分析，得出以下結(jié)論:

(1)彈炮耦合過程對炮口橫向振動和縱向振動影響規(guī)律相當。在一發(fā)彈的大約6 ms的膛內(nèi)運動過程中，炮口各周向單元的徑向位移變化規(guī)律相同，相對于沖擊碰撞，彈丸重力因素影響較小;

(2)彈丸前定心部對炮口振動影響很大，使振動幅值和頻率都增大。由于彈丸膛內(nèi)高速運動的章動性，使彈丸前定心部趨于不碰撞，而本文的數(shù)學模型和圖8證明了由于入膛時隨機碰撞引起的一系列規(guī)律碰撞是存在的。

本文建立了適合小口徑火炮的較為精確的彈炮耦合數(shù)學模型，為建立發(fā)射動力學模型和射擊精度分析提供了彈炮方面的理論基礎(chǔ)，并對其他位置的碰撞分析具有一定的借鑒意義。

[1]周 叮，謝玉樹.彈丸膛內(nèi)運動引起炮管振動的小參數(shù)解法[J].振動與沖擊，1999，18(1):76 －81.ZHOU Ding，XIE Yu-shu.Solution of barrel vibration caused by motion of bullets in bore by use of small parameter method[J].Journal of Vibration and Shock，1999，18(1):76 －81.

[2]劉 雷，陳運生，楊國來.基于接觸模型的彈炮耦合問題研究[J].兵工學報，2006，27(6):985 －987.LIU Lei，CHEN Yun-sheng，YANG Guo-lai.A study on the projectile-barrel coupling based on contact model[J].Acta Armamentarii，2006，27(6):985 －987.

[3]葛建立，楊國來，陳運生，等.基于彈塑性接觸/碰撞模型的彈炮耦合問題研究[J].彈道學報，2008，20(3):103－106.GE Jian-li，YANG Guo-lai，CHEN Yun-sheng，et al.A study on projectile-barrel coupling problem based on elastoplastic contact/impact Model[J]. Journal of Ballistics，2008，20(3):103－106.

[4]劉 寧，楊國來.彈管橫向碰撞對身管動力響應(yīng)的影響[J].彈道學報，2010，22(2):67 －70.LIU Ning，YANG Guo-lai.Effect of lateral impact between projectile and barrel on dynamic response of tube[J].Journal of Ballistics，2010，22(2):67 －70.

[5]陳 強，楊國來，王曉鋒，等.變速移動彈簧阻尼質(zhì)量系統(tǒng)作用下梁的動態(tài)響應(yīng)[J].動力學與控制學報，2010，8(3):254－257.CHEN Qiang， YANG Guo-lai， WANG Xiao-feng，et al.Dynamic responses of beam subjected to a moving oscillator with non-uniform velocity[J].Journal of Dynamics and Control，2010，8(3):254 －257.

[6]Nakalswarny K K.Finite element analysis of a projectile during gun launch[D]. BachelorofEngineering in Mechanical Engineering B.I.E.T，Kuvempu University，Karnataka，2002.

[7]孫河洋，馬吉勝，張高明，等.彈帶擠進過程中膛線的應(yīng)力分析[J].軍械工程學院學報，2008，20(4):43 －46.SUN He-yang， MA Ji-sheng， ZHANG Gao-ming， et al.Analysis of the land stress during the engraving process[J].Journal of Ordnance Engineering College，2008，20(4):43－46.

[8]彭 濤，王學軍，黃善文.基于ANSYS的彈帶擠進變形及應(yīng)力分析[J].艦船電子工程，2009，29(11):156 －159.PENG Tao，WANG Xue-jun，HUANG Shan-wen.Analysis on material deformation and stress in the process of driving band engraving based on ANSYS[J].Ship Electronic Engineering，2009，29(11):156 －159.

[9]吳 宏.彈丸膛內(nèi)外運動及其對射擊精度影響研究[D].南京:南京理工大學，2001.

[10]王少欽，夏 禾，郭薇薇，等.變速移動荷載作用下簡支梁橋的動力響應(yīng)及共振分析[J].振動與沖擊，2010，29(2):26－28.WANG Shao-qin，XIA He，GUO Wei-wei，et al.Dynamic response and resonance analysis for a simply-supported bridge under speed-varying loads[J].Journal of Vibration and Shock，2010，29(2):26 －28.

[11]Liu K，Reynders E，Roeck G D，et al.Experimental and numerical analysis of a composite bride for hig-speed trains[J].Journal of Sound and Vibration，2009，320:201 － 220.

[12]史躍東，王德石.考慮慣性效應(yīng)的移動彈丸作用下身管振動特性[J].兵工學報，2011，32(1):414 －419.SHIYue-dong， WANG De-shi. Study on vibration characteristicsofbarrelsubjected to moving projectile considering inertia effect[J].Acta Armamentarii，2011，32(1):414－419.

[13]黃建亮，陳樹輝.縱向與橫向振動耦合作用下軸向運動梁的非線性振動研究[J].振動與沖擊，2011，30(8):24－27.HUANG Jian-liang，CHEN Shu-hui. Study on nonlinear vibration of an axially moving beam with coupled transverse and longitudinal motions[J].Journal of Vibration and Shock，2011，30(8):24 －27.