李嘉欣，沈大偉，馬鐵華，張煜，裴東興

(中北大學(xué) 儀器科學(xué)與動態(tài)測試教育部重點實驗室，山西 太原 030051)

彈丸在飛行的過程中，由于發(fā)射條件、物理環(huán)境以及彈丸自身的設(shè)計等原因，使彈丸外彈道產(chǎn)生誤差，影響飛行的穩(wěn)定性。其中，彈丸轉(zhuǎn)速及章動角成為影響彈丸飛行性能的最主要參數(shù)[1]。

隨著武器技術(shù)的發(fā)展，彈丸的外形變化很大，作用在彈丸本體的各種空氣動力以及力矩已經(jīng)不可忽略。因此，必須考慮影響彈丸在飛行過程中所有的空氣動力才能更加準(zhǔn)確地模擬彈丸的運動姿態(tài)，通過數(shù)學(xué)模型的建立，分析彈道和射表的設(shè)計編制以及實驗射程[2]。應(yīng)用理想軸對稱的彈丸，研究其轉(zhuǎn)速和章動角的變化，從而更加準(zhǔn)確地研究彈丸的飛行穩(wěn)定性[3]。

針對以上問題，筆者提出應(yīng)用MATLAB軟件求解無控旋轉(zhuǎn)彈丸外彈道的仿真模型，其結(jié)果與彈載測試儀現(xiàn)場試驗結(jié)果比對，找到彈丸的運動規(guī)律，為后續(xù)的武器系統(tǒng)的優(yōu)化設(shè)計以及新型武器的研發(fā)奠定理論基礎(chǔ)。

1 彈丸外彈道運動姿態(tài)模型

彈丸運動時受到擾動會形成不斷變化的攻角δ，當(dāng)攻角增大時，導(dǎo)致旋轉(zhuǎn)彈丸飛行不穩(wěn)定，因此研究運動中彈丸的姿態(tài)以及運動規(guī)律，將彈丸視為自由運動剛體建立一個數(shù)學(xué)仿真模型[4]，分析計算其外彈道參數(shù)以及穩(wěn)定性。

根據(jù)旋轉(zhuǎn)彈丸飛行的姿態(tài)建立坐標(biāo)系,如圖1所示，以彈丸質(zhì)心O為原點，射擊方向為Ox正向，建立地面坐標(biāo)系Oxyz，以速度v方向為Ox1正向，建立彈道速度坐標(biāo)系為Ox1y1z1，以穿過質(zhì)心彈丸軸線為Ox2正向，建立彈軸坐標(biāo)系為Ox2y2z2。根據(jù)彈丸受力可將其飛行運動的物理過程分解為質(zhì)心運動和繞質(zhì)心轉(zhuǎn)動兩種形式[5]，根據(jù)建立的坐標(biāo)系，列出主要運動方程。

1.1 彈丸質(zhì)心運動方程

當(dāng)彈丸在高速旋轉(zhuǎn)時，其圓錐運動服從質(zhì)心運動定理，為

(1)

式中：v為彈丸的速度；θa為速度俯仰角；θb為速度偏航角；Fx1、Fy1、Fz1分別為F在Ox1y1z1坐標(biāo)軸下的分力。

1.2 彈丸繞質(zhì)心運動方程

彈丸繞質(zhì)心的旋轉(zhuǎn)過程中，根據(jù)動量矩定理可得：

(2)

(3)

1.3 彈丸運動姿態(tài)模型

將彈丸的射程對時間t微分得到彈丸速度在地面坐標(biāo)系下的分量，即

(4)

(5)

式中：δ2為偏航攻角；δ1為俯仰攻角；β為坐標(biāo)系轉(zhuǎn)換角[6]。

將以上彈丸質(zhì)心運動方程(1),繞質(zhì)心運動方程(2)、(3)，彈丸速度在地面坐標(biāo)系下的分量(4),以及根據(jù)各方位角之間的關(guān)系得出的補(bǔ)充方程(5)聯(lián)立起來就可以得到彈丸運動姿態(tài)的數(shù)學(xué)模型。

2 彈丸外彈道仿真

為了方便求解并提高計算精度，采用龍格-庫塔法對仿真模型進(jìn)行解算。

2.1 龍格-庫塔法

龍格-庫塔法具有精度高且穩(wěn)定的特點，在工程技術(shù)上廣泛應(yīng)用的一種算法[7]，用于數(shù)值求解微分方程。針對彈丸運行姿態(tài)的數(shù)學(xué)模型，有

(6)

(7)

式中，點m處的值為(tm,y1n,y2n,…,ymn),參量定義詳見文獻(xiàn)[7]。

2.2 MATLAB程序的設(shè)計以及編寫

應(yīng)用龍格-庫塔法在MATLAB軟件中對微分方程進(jìn)行求解，程序框圖如下：

2.3 仿真結(jié)果及分析

對彈丸及環(huán)境的理想假設(shè)：

1)理想彈丸的慣性主軸為彈丸軸線且外形、質(zhì)量相對于彈軸分布是對稱的；

2)彈丸是靜不穩(wěn)定的狀態(tài)；

3)地表重力加速度鉛直向下，大小不變；

4)統(tǒng)一使用我國炮兵標(biāo)準(zhǔn)氣象條件。

設(shè)置初始速度俯仰角θa=21°，速度偏航角θb=0°，初始章動角δ0=3.92°，假設(shè)彈軸俯仰角φ1=24.92°，彈軸偏航角φ2=0°，在MATLAB中仿真計算，得到轉(zhuǎn)速和章動角隨時間變化曲線。

彈丸的初始轉(zhuǎn)速為390 r/min，在極阻尼力矩的作用下，自轉(zhuǎn)速度逐漸減小，彈丸轉(zhuǎn)速變化曲線如圖3所示。

彈丸的俯仰攻角以及偏航攻角的變化規(guī)律如圖4所示，其峰值近似呈指數(shù)規(guī)律減小，后期減緩趨于穩(wěn)定。

由定義可知，彈丸攻角就是彈丸章動角[8]，其彈丸的章動角變化如圖5所示。

將圖5的前3 s波形依次展開，如圖6～8所示。1次章動運動即1個周期的變化，由圖6～8可知彈丸章動周期在1、2、3 s內(nèi)分別是32.5、34.0、36.0 ms，彈丸章動周期逐漸增大，飛行趨于穩(wěn)定。

3 實測信號分析

3.1 彈丸的轉(zhuǎn)速

彈載測試儀[9]所采集到的彈丸轉(zhuǎn)速信息經(jīng)濾波去噪處理后如圖9所示。

根據(jù)圖9，按轉(zhuǎn)速測試原理得到的實測彈丸轉(zhuǎn)速值在0、1、2、3 s時刻分別為390.0、378.8、370.4、344.8 r/min。

3.2 彈丸的章動角

提取圖9中轉(zhuǎn)速曲線中包絡(luò)線，如圖10所示，包絡(luò)線中波動就代表彈丸的章動運動。根據(jù)T=t/n對濾波后的曲線進(jìn)行解算即可得到彈丸在飛行過程中的章動周期[10]。將圖中信號分別在1、2、3 s附近展開，如圖11～13所示。數(shù)據(jù)表明彈丸章動在第1、2、3 s內(nèi)的周期分別是33.34、33.40、33.8 ms。

濾去無效信號，通過計算獲得彈丸在不同時刻的章動角。根據(jù)測試儀不同時刻的輸出(mV)及對應(yīng)轉(zhuǎn)速，可得圖14，彈丸在出炮口時刻的地磁感應(yīng)電動勢為E=56 mV，彈丸章動變化的電動勢為ΔE=6.5 mV，彈丸軸線與當(dāng)?shù)氐卮艌鲋g的夾角為30°，彈丸在剛出炮口時兩電動勢之間有如下關(guān)系[8]：

(8)

根據(jù)式(8)可以計算得出彈丸剛出炮口時刻的章動角δ=3.92°。依次求解出彈丸在1、2、3 s時的章動角分別為2.464°、2.329°、2.205°。

將彈丸轉(zhuǎn)速、章動角、章動周期的仿真結(jié)果與彈載測試儀實測進(jìn)行對比，如圖15所示。

通過對圖15的數(shù)據(jù)計算，可得到仿真值與實測值誤差在10%以內(nèi)，相關(guān)系數(shù)均大于0.94，因此所建立的彈丸外彈道運動姿態(tài)仿真模型是合理可行的。

4 結(jié)束語

筆者以無控旋轉(zhuǎn)彈丸為主體，建立了外彈道數(shù)學(xué)模型，并基于MATLAB軟件對彈丸的運動飛行姿態(tài)進(jìn)行了仿真分析，得到了轉(zhuǎn)速和章動曲線，并將速度、章動角以及章動周期與實測數(shù)據(jù)進(jìn)行對比，所得結(jié)果基本一致，驗證了筆者所建立姿態(tài)模型的合理性。