末敏子彈組合體全姿態運動仿真

王海峰，賈 波，王志成

(中國人民解放軍63850部隊， 吉林 白城 137001)

末敏彈藥是新型智能彈藥，通過大口徑榴彈、火箭彈、機載布撒器等運載載體投放，能夠對集群裝甲武器和火炮目標高效毀傷。末敏彈藥集子母彈技術、敏感器技術和爆炸成形彈丸(EFP)技術于一身，能夠同時對多個目標自動探測、識別、攻擊，具備“打了不管”等優點。但這種打后不用管是指在末敏子彈掃描區域內的自動尋的、自動攻擊，如果投放不準，末敏子彈的掃描區域內沒有目標，其所有的能力、具備的優點都將無法實現。因此依賴精確的射表使末敏子彈準確投放成為末敏彈藥充分發揮性能的基本前提。末敏彈空中飛行機構動作頻繁，彈道特性復雜，特別是為減小末敏子彈開傘動載，大型火箭末敏彈在開艙拋撒前先將穩定尾翼分離，使多枚末敏子彈組成的戰斗部組合體進行大攻角減速運動。減速運動中，姿態角變化超過90°。用傳統的剛體運動方程[1]描述飛行姿態，微分方程出現奇點，計算發散溢出。質點彈道模型[2]忽略了姿態運動，有較好的計算穩健性，但由于建模假設與實際情況明顯不符，計算精度差。四元數法[3-7]從理論上講比較完美，但實際應用中存在較大的累積計算誤差。本文利用正反歐拉方程間精華區倒掛關系進行分區交替運算，描述末敏彈戰斗部組合體的全姿態運動，在此基礎上仿真分析質點運動假設的適用性，為此類彈箭彈道仿真與射表編擬模型的確定提供理論依據。

1 全姿態運動方程建立與求解

1.1 坐標系及相互關系

彈體運動可分為質心的平動運動和圍繞質心的轉動運動。質心運動可建立在速度坐標系中，繞心運動可建立在彈體坐標系中。定義地面坐標系以發射點O為原點，OX軸沿發射點水平線指向射擊方向，OY軸鉛直向上，OZ軸按右手法則確定。基準坐標系OC-XYZ可由地面坐標系平移至彈體質心形成，隨質心平動。速度坐標系OC-XCYCZC原點OC取在彈體質心，OCXC軸沿速度V方向，OCYC軸垂直OCXC，指向上方為正，速度坐標系可由基準坐標系經2次旋轉得到：先繞OCY軸旋轉ψ角，后繞OCZC軸旋轉θ角，逆時針旋轉為正。彈體姿態通常用固連在彈體上的彈體坐標系相對基準坐標系的歐拉角描述。彈體坐標系OC-ξηζ以彈體質心為原點，由基準坐標系經3次轉動而成，第1次繞第2軸(OCY)旋轉φ2，第2次繞第3軸(OCZ′)旋轉φ1，第3次繞第1軸(OCξ)旋轉γ，逆時針旋轉為正，如圖1所示。

圖1 基準坐標系與彈體坐標系關系示意圖

上述實例只是一種旋轉方案，由于3個坐標軸旋轉的先后次序有6種不同組合，即123、132、213、231、312、321方式，每個坐標軸旋轉方向有逆時針和順時針2種可能，使得同一彈體姿態可以有48種歐拉角描述方式[8]。

(1)

1.2 質心運動方程的建立

在速度坐標系中，質心動力學方程為：

(2)

質心運動學方程為：

式(2)中：m為彈體質量；FRXc、FRYc、FRZc分別為空氣阻力在速度坐標系3軸上的分量；FLXc、FLYc、FLZc為空氣升力在速度坐標系3軸上的分量；GXc、GYc為重力在速度坐標系前兩軸上的分量。阻力和升力隨飛行馬赫數和攻角δ變化。

1.3 繞心運動雙歐拉方程的建立

繞心動力學方程為：

(3)

繞心運動學正歐拉方程為:

sinδ1=sinφ1cosθ-cosφ1sinθcos(φ2-ψ)

sinδ2=cosφ1sin(φ2-ψ)/cosδ1

cosδ=cosδ1cosδ2

繞心運動學反歐拉方程為:

式(3)中：ωξ、ωη、ωζ為彈體角速度在彈體坐標系三軸上的分量；Mxzξ為極阻尼力矩；Mzη、Mzζ為翻轉力矩在第2、3軸上的分量；Mzzη、Mzzζ為赤道阻尼力矩在第2、3軸上的分量；Myη、Myζ為馬格努斯力矩在第2、3軸上的分量；C、A為彈體極、赤道轉動慣量。

1.4 兩種歐拉角互換

從繞心運動方程可以看出，當彈體坐標系第2次旋轉角度為±90°時，其余弦函數值為0，方程存在奇異點，計算發散，在±90°附近由于余弦函數值很小，作為除數放大了計算誤差，因此無法準確確定第2次及第3次旋轉的角度。

若彈體坐標系按231旋轉方式，基準坐標系到彈體坐標系轉換矩陣見式(1)，按321方式，基準坐標系到彈體坐標系轉換矩陣為：

(4)

彈體的每次轉動帶來彈體姿態的變化，但對應某一固定彈體姿態，可以找到多個旋轉角度值，這些值以2π為周期，為確保求解的唯一性，將旋轉角度定義在[-π，π]，將正反歐拉方程轉換的時機定為第二次旋轉的姿態角位于[-3π/4，-π/4)∪[π/4，3π/4)區間。則由231旋轉方式轉換至321旋轉方式的反歐拉角為：

由321旋轉方式轉換至231旋轉方式的正歐拉角為：

K取0和1，對應于[-π，π]內的反三角函數的兩組解。

彈體姿態變化具有連續性，所以當殘差和Δφ10+Δφ20+Δγ0>Δφ11+Δφ21+Δγ1時，取K為1時的那組解，反之取K為0時的那組解。

1.5 求解流程

某火箭末敏彈點火后，彈載時間裝置開始計時，母彈在空中自由飛行至預置分離時間t0，穩定尾翼與戰斗部組合體分離，戰斗部組合體大攻角飛行固定時間ta后，開艙拋撒末敏子彈，進入末敏子彈飛行段。采用4階龍格-庫塔法數值求解彈道微分方程，積分步長為h，正歐拉方程解算區間為D=[-π，-3π/4)∪[-π/4，π/4)∪[3π/4，π]。求解流程如圖2所示。

圖2 全姿態運動方程求解流程框圖

2 全姿態運動的仿真與阻力特性分析

2.1 姿態運動規律仿真分析

利用上述的全姿態運動方程與某型箭體氣動數據，對戰斗部組合體從母彈分離(拋掉尾翼穩定裝置)后的大攻角減速運動進行仿真。30°仰角射擊時正歐拉方程中的彈軸擺角、攻角與飛行速度計算結果見圖3、圖4。從圖3可以看出，在5 s運動中，組合體出現了翻滾運動。剛分離后，彈體能夠保持原來的姿態運動，當組合體出現擾動產生攻角后，由于壓心處于質心與彈頂之間，在翻轉力矩作用下，彈體迅速翻轉，0.76 s左右沿第3軸正向擺動達到180°。當攻角達到180°時翻轉力矩為0，但擺動速度達到極值，此后，翻轉力矩改變方向，阻止組合體擺動，擺動速度逐漸減小，當擺動速度為0時，在翻轉力矩作用下彈體沿第3軸負向擺動。在大攻角運動下，飛行速度從1.5個馬赫數降到0.5個馬赫數，衰減迅速。而在低速、大攻角運動中，彈體的壓心在質心前后來回變動，導致第3軸呈周期性振蕩，攻角在100°附近作周期性變化。

圖3 231旋轉方式下擺角、攻角變化規律曲線

圖4 全姿態運動飛行速度變化規律曲線

采用不同仰角進行彈道仿真，考察不同仰角情況下戰斗部組合體攻角變化情況，結果見圖5。

不同仰角射擊，戰斗部組合體運動攻角隨時間的變化特性整體上有一定規律，剛分離時較小，而后彈體翻滾運動，攻角迅速上升，以100°為中心反復振蕩；但局部差異較大，振蕩起始時間不同，分離時仰角越大、速度越高振蕩越早，振蕩周期與攻角幅值也不相同。

2.2 阻力特性分析

通過流場數值模擬得到戰斗部組合體氣動參數。阻力系數隨馬赫數、攻角的變化規律見圖6、圖7。與通常情況下彈箭阻力系數一樣，戰斗部組合體阻力系數與飛行馬赫數和攻角有關，阻力系數最大值是最小值的40余倍，峰值位于1.2馬赫、90°攻角左右，但攻角變化范圍很大，阻力系數隨攻角的變化遠遠大于其隨馬赫數的變化。全姿態運動情況下，攻角對阻力的影響將超過飛行速度，成為影響運動的最重要因素。

圖6 不同攻角下阻力系數隨馬赫數變化規律曲線

圖7 不同馬赫數下阻力系數隨攻角變化規律曲線

不同仰角情況下戰斗部組合體阻力系數變化情況見圖8、圖9。

不同仰角射擊時，相同飛行速度下攻角未必相同，所以阻力系數不同，阻力系數隨馬赫數的變化有較大差異；阻力系數隨時間的變化特性整體上有一定規律，剛分離時阻力較小，而后由于攻角的變化，阻力系數也隨之振蕩，但振蕩周期與幅值都不相同。不同仰角情況下，阻力系數隨時間的變化規律、隨馬赫數的變化規律均不相同。

圖8 不同仰角下阻力系數隨飛行時間變化規律曲線

圖9 不同仰角下阻力系數隨馬赫數變化規律曲線

2.3 戰斗部組合體質點運動模型計算誤差驗證

以往彈箭穩定飛行情況下，認為攻角很小，誘導阻力影響不大，很多情況下僅考慮阻力系數隨馬赫數的變化，可用質點運動模型與平均阻力系數進行彈道仿真。按這一通常做法，在不同仰角下用質點模型計算出的開艙點彈道諸元與全姿態運動方程仿真結果作差，得到各主要彈道諸元偏差，結果見表1-表3。其中開艙點坐標是對分離點的相對坐標。驗證中，將1.2節中的質心運動方程忽略升力后作為質點模型。將60°仰角下的阻力系數擬合成馬赫數的線性方程，作為平均阻力系數，如圖9所示。

表1 全姿態運動模型下開艙點諸元

表2 質點模型下開艙點諸元

表3 2種運動模型下彈道諸元偏差

從驗證結果可以看出，除橫偏外，質點模型誤差大，而且不同仰角下偏差變化較大，甚至符號相反。雖然可以通過放大或縮小阻力系數消除某一特定射擊條件下的彈道計算偏差，但是這種簡單調整帶來的效果難以具有普遍適用性，當射擊條件偏離特定條件時計算誤差又會突顯出來。戰斗部組合體減速運動過程中，攻角呈周期性運動，且攻角均值與幅值較大，攻角不能忽略，而且攻角運動規律也難以用簡單的形式準確定量描述，因此質點彈道模型無法適用這類彈道準確仿真。

3 結論

采用雙歐法建立了剛體全姿態運動方程，給出了坐標系2種旋轉方式下姿態角的互換公式，并進行了實例仿真驗證與運動規律分析。仿真結果表明：末敏子彈組合體從母彈分離后，彈體先翻轉，而后攻角在100°附近作周期性振蕩，飛行阻力大，組合體飛行速度衰減迅速；不同仰角下，彈體運動姿態存在明顯差異，攻角幅值與周期均不同，飛行過程中阻力系數變化劇烈，其規律難以用數學模型描述，因此，不考慮姿態運動，采用質點模型和平均阻力計算彈道必將產生大誤差。本文采用雙歐法有效解決了大攻角運動彈道仿真問題，為此類彈箭彈道仿真與射表編擬提供了必要的技術支撐。