基于概率密度演化方法的火炮輸彈過程不確定性分析

王明明，錢林方，2，陳光宋，劉太素

(1.南京理工大學 機械工程學院，江蘇 南京 210094；2.西北機電工程研究所，陜西 咸陽 712099；3.南京工程學院，江蘇 南京 211167)

0 引言

輸彈機是大口徑火炮的重要組成部分，其作用是將彈丸從輸彈出發(fā)位置沿炮膛軸線方向快速、可靠地推送到炮膛內。輸彈過程包括強制輸彈過程和慣性輸彈過程，強制輸彈過程是指彈丸在推彈板強制推動下沿炮膛軸線方向向前運動，直至彈丸底部脫離推彈板；慣性輸彈過程是指彈丸底部與推彈板脫離，在慣性作用下卡進坡膛。由于輸彈系統(tǒng)中存在輸彈機齒間間隙、推彈板速度、彈丸慣量等大量隨機因素，導致輸彈結束時彈丸速度和姿態(tài)存在不確定性，影響彈丸的卡膛一致性，進而通過火炮發(fā)射過程的傳播影響射擊精度。因此，對輸彈過程進行誤差傳遞分析，開展輸彈過程參數(shù)不確定性對輸彈一致性的影響研究，對進一步提高輸彈一致性、提升火炮射擊精度有著重要意義。

隨著火炮自動化的發(fā)展，國內相關學者對輸彈過程進行了研究。趙森等以試驗為基礎，分析研究了恒定輸彈力對初速的影響，提出了改進措施。薛百文等通過彈丸受力分析建立運動方程，研究了底盤運動狀態(tài)對彈丸在慣性輸彈段的影響。劉朋展等針對彈丸姿態(tài)影響輸彈可靠性與一致性問題，分析了影響彈丸姿態(tài)的因素。近年來，學者們開始更多地關注參數(shù)不確定性對輸彈結果的影響。李偉等考慮隨機參數(shù)漂移，探討了供輸彈系統(tǒng)的動作可靠性問題。高學星等提出了一種間接提高彈藥協(xié)調器定位精度可靠性的優(yōu)化方法，具有較好的優(yōu)化效果和較高的計算效率。劉太素等考慮輸彈協(xié)調機構中的眾多不確定性參數(shù)，提出了基于稀疏多項式混沌展開(SPCE)- 高維模型表示(HDMR)的穩(wěn)健優(yōu)化方法，實用有效。林通等基于ADAMS軟件建立考慮參數(shù)隨機性的輸彈過程動力學模型，進行了輸彈機穩(wěn)健優(yōu)化設計研究。張弘鈞等基于區(qū)間參數(shù)的不確定性分析方法對某大口徑火炮輸彈一致性進行研究，認為影響輸彈一致性的主要因素為慣性階段的參數(shù)不確定性。相關研究結果為分析參數(shù)不確定性對輸彈過程一致性的影響提供了重要參考。然而，如何考慮更多的參數(shù)以及提高分析的計算精度是需要進一步研究的熱點問題。針對不確定性參數(shù)數(shù)量眾多、計算效率低下的問題，發(fā)展一種考慮高維不確定性參數(shù)的輸彈不確定性分析方法，有助于進一步量化不確定性參數(shù)對輸彈過程的影響，也可為提高卡膛一致性研究提供良好的基礎。

輸彈過程的不確定性分析在理論上屬于不確定性傳播問題，常采用概率類不確定性分析方法進行研究，其數(shù)學理論較為完善且不確定性量化和傳遞的推導也更加準確便捷。概率方法通過隨機變量對不確定性進行描述分析，最經典的即蒙特卡洛(MC)方法，具有精度高、適用范圍廣的優(yōu)點，但其計算量龐大，一般僅用于驗證結果的準確性。實際上，大多概率類方法如代理模型、函數(shù)展開、降維積分等都是基于系統(tǒng)響應的統(tǒng)計矩計算其概率密度函數(shù)，隨著參數(shù)數(shù)量的增加，將面臨計算量龐大、高階矩計算精度不足、概率分布擬合不準確等問題。概率密度演化方法(PDEM)通過引入擴展狀態(tài)向量構造復合系統(tǒng)隨機動力學方程，建立系統(tǒng)響應的概率密度演化方程，直接積分求解系統(tǒng)響應的概率分布，精度高且對樣本量需求小，可以有效避免“維度災難”的發(fā)生。實際中，輸彈過程涉及到的不確定性參數(shù)較多，基于概率密度演化方法對其進行不確定性分析是較好的選擇。

本文通過分析某火炮輸彈機的結構與原理，利用ADAMS軟件建立輸彈機動力學模型，推導了輸彈響應的概率密度演化方程，進而通過Lax-Wendroff格式的有限差分法進行求解。將仿真結果與文獻[20]試驗結果對比，驗證了模型的準確性；將PDEM與MC方法進行對比，驗證了PDEM的適用性；對不同工況下的輸彈過程不確定性傳播進行了分析。

1 輸彈過程動力學建模

本文以某中大口徑火炮齒輪- 齒條式輸彈機為研究對象，對其進行動力學建模。

1.1 基本假設

為便于建立輸彈機動力學模型，分析輸彈過程，做如下假設：

1)在輸彈過程中彈丸不變形，可將其簡化成6自由度剛體；

2)將彈丸與托彈板、炮尾以及身管的碰撞考慮為彈性接觸碰撞；

3)協(xié)調到位后無其他外部激勵，輸彈機與身管相對于地面保持固定。

1.2 輸彈過程動力學建模

輸彈機包含協(xié)調油缸、輸彈油缸、齒輪、托彈板、推彈板、協(xié)調臂箱體等機構，具體結構如圖1所示。協(xié)調臂的一端安裝在耳軸上，協(xié)調油缸筒端與架體相連，協(xié)調油缸桿端與協(xié)調臂相連，三處連接均可旋轉；箱體固定在協(xié)調臂上，輸彈油缸整體安裝在箱體內，油缸筒端固定在箱體上，油缸桿與齒輪軸固連；齒輪通過軸承與齒輪軸連接，同時與上下兩根齒條嚙合，上齒條固連在托彈板上，下齒條與箱體固接；托彈板上有抱鎖裝置，用于防止彈丸掉落。

圖1 輸彈機結構圖Fig.1 Structural view of ammunition ramming mechanism

根據(jù)輸彈機結構以及各部件之間的約束關系，可分析得到輸彈機拓撲關系如圖2所示。

圖2 輸彈機拓撲結構Fig.2 Topological structure of ammunition ramming mechanism

輸彈機原理如下：通過電控輸彈油缸推動齒輪軸向身管方向平移，由于與齒輪嚙合的上下齒條分別固連在輸彈機箱體和托彈板上，而輸彈機箱體相對于地面固定不動，故此時承載著彈丸的托彈板向前平移運動。具體過程又分為強制輸彈和慣性輸彈兩部分，強制輸彈是指彈丸在托彈板的推動下向前加速平移運動達到最大速度的過程，慣性輸彈是指彈丸不再受到托彈板的推力作用，而僅靠慣性向前運動直至與托彈板分離，繼而飛入炮膛完成卡膛的過程。其中，強制過程中彈丸與托彈板相對靜止，是一個簡單加速運動過程；慣性過程中首先是彈丸在托彈板摩擦力與重力作用下作半約束慣性運動，其次是僅受重力作用的慣性飛行階段，最后是重力、炮膛摩擦力以及坡膛碰撞力作用下的膛內運動階段。

當彈丸與托彈板、身管坡膛發(fā)生摩擦時，摩擦力以滑動摩擦為主?？紤]到輸彈過程中的彈丸速度一般不超過5 m/s，屬低速摩擦，符合如下摩擦力經典公式：

=

(1)

式中：為摩擦力；為接觸正壓力；為摩擦系數(shù)，可以直接按照黃銅- 鋼的摩擦系數(shù)設定彈帶與托彈板、身管坡膛的摩擦系數(shù)，按照鋼- 鋼的摩擦系數(shù)設定彈丸定心部與托彈板、身管坡膛的摩擦系數(shù)。

當彈丸與托彈板以及坡膛發(fā)生接觸碰撞時，按照彈性碰撞假設可直接引用Hertz接觸定律計算其接觸剛度

(2)

(3)

、分別為兩個物體在接觸點的曲率半徑；、為兩種材料的泊松比；、為兩種材料的彈性模量。

輸彈開始時，彈丸在重力作用下靜止在托彈板上，彈底與推彈板接觸，彈帶、定心部與托彈板接觸；然后通過輸彈油缸施加載荷壓力，推動齒輪- 齒條運動，進而帶動托彈板、彈丸向前運動，直至彈丸抵達卡膛起始位置。參考基本假設，本文根據(jù)某型中大口徑火炮輸彈機構尺寸利用ADAMS軟件建立輸彈過程動力學模型如圖3所示，該系統(tǒng)可輸出彈丸的位移、速度和加速度曲線等。

圖3 輸彈機ADAMS模型Fig.3 ADAMS model of ammunition ramming mechanism

2 輸彈過程概率密度演化方程及求解

2.1 輸彈過程概率密度演化方程

實際工程中影響輸彈精度的不確定性因素眾多，為方便建模分析，可將輸彈過程涉及到的不確定性因素簡單分為2類：系統(tǒng)隨機參數(shù)和輸入隨機參數(shù)。由于本文研究對象是一種彈射式輸彈機，故系統(tǒng)隨機參數(shù)主要指彈丸和輸彈機的物理、幾何參數(shù)，具體包括彈丸的質量、質心位置、3個方向的轉動慣量、前定心部半徑、前定心部與托彈板間的摩擦系數(shù)、彈帶半徑、彈帶與托彈板間的摩擦系數(shù)、托彈板半徑、托彈板軸線與身管軸線的傾角以及豎直距離，共12個參數(shù)；輸入隨機參數(shù)即輸彈油缸的壓力參數(shù)，主要包含最大加載壓力、最大壓力時刻、最小加載壓力、最小壓力時刻，共4個參數(shù)。

考慮上述不確定性參數(shù)，可將輸彈隨機動力學方程表示成如下形式

(4)

令=[,]，記其聯(lián)合概率密度函數(shù)為()，可見同時包含了輸彈過程中的16個隨機參數(shù)，則彈丸響應可表示如下：

=(,)=[(,),(,),…,(,)]

(5)

(6)

由(4)式可知，上述輸彈隨機動力學方程包含了輸彈過程中的所有隨機因素，可以看作一個保守的隨機系統(tǒng)，符合概率密度守恒原理的描述。因此，輸彈過程中的彈丸隨機響應與隨機變量的聯(lián)合概率密度函數(shù)(,,)應滿足廣義概率密度演化方程，可寫成如下形式：

(7)

即輸彈過程廣義概率密度演化方程。

2.2 方程求解

常用的邊界條件為

(,,)|→±∞=0,=1,2,…,6

(8)

初始條件為

(,,)|=0=(-)()

(9)

式中：(·)為狄拉克函數(shù)；為初始時刻的彈丸狀態(tài)變量值，即彈丸在托彈板上達到靜平衡時的位姿參數(shù)。由于彈丸靜置在托彈板上時僅彈底與托彈板接觸、前定心部和彈帶與托彈板接觸，故可以根據(jù)彈丸、托彈板尺寸直接計算。若以托彈板軸線與推彈板交點為坐標原點，彈丸前進方向為軸方向，軸方向垂直軸向上，軸方向由右手法則確定，則=[0344 3,-0588 3,0,0,-1004,0]。

求解(7)式可得到的概率密度函數(shù)：

(10)

式中：為隨機變量的參數(shù)空間。

特別地，當只關心第個響應時，對(,,)求除以外的多重積分，得到如下邊緣概率密度

(11)

同理，對(7)式進行多重積分，得

(12)

由此可見，當考察的物理量為某單一響應時，其與隨機參數(shù)的聯(lián)合概率密度函數(shù)(,,)關于時間的變化率與關于空間的變化率呈正比，比例系數(shù)由該項位移響應的瞬時速度決定，體現(xiàn)了輸彈過程中概率演化與物理演化之間的內在規(guī)律。

(12)式是1階偏微分方程，由于無法給出顯式的表達式，只能通過數(shù)值算法進行求解。對隨機參數(shù)進行采樣后獲得一系列確定值，此時(12)式可以寫成

(13)

(13)式符合1階偏微分方程的一般形式

(14)

式中：(,)看作(13)式中的聯(lián)合概率密度函數(shù)；()為速度項；為位移項。此時可直接采用有限差分法對其進行求解。

首先對狀態(tài)空間域和時間域進行離散化，

=·Δ;=0,±1,±2,…

(15)

=·Δ;=0,1,2,…

(16)

式中：Δ、Δ分別表示空間離散步長和時間離散步長；為第步的空間坐標；為第步的時間坐標，不同的離散化格式對應不同的差分格式。

本文采用如下基于中心差分格式的Lax-Wendroff格式：

(17)

同時，為了計算的收斂性和穩(wěn)定性，空間以及時間離散步長的選取還需始終滿足

0<|()|≤1

(18)

這樣才能保證數(shù)值計算的傳播始終收斂，即著名的Courant-Friedrichs-Lewy(CFL)條件。

綜上，本文采用基于Lax-Wendroff格式的有限差分算法計算輸彈過程彈丸響應分布情況，具體步驟如下：

在輸彈過程隨機變量的參數(shù)空間內選取樣本點。根據(jù)各隨機變量的分布類型及統(tǒng)計信息，利用拉丁超立方采樣選取個離散代表點集={,,…，，…,}，并計算每個代表點的概率權重，又稱賦得概率。

建立時空差分網絡，利用基于Lax-Wendroff格式的有限差分法求解(12)式，得到第個彈丸響應的第(=1,2,…,)條軌跡的邊緣概率密度函數(shù),(,,)。

累計求和，將步驟3中的邊緣概率密度函數(shù)累加，得到第個彈丸響應隨時間演化的概率密度函數(shù)

(19)

3 輸彈過程不確定性分析

3.1 輸彈模型驗證

通過試驗測量某大口徑車載火炮在30°射角時的輸彈油缸壓力隨時間變化的曲線，結果如圖4(a)所示，油缸壓力變化大體分成兩個階段，即首次加載- 卸載階段以及二次加載- 穩(wěn)定階段。實際輸彈過程中，彈丸往往在油缸首次加載- 卸載階段即獲得最大速度，之后伴隨著托彈板減速，彈丸底部與推彈板分離，進入慣性輸彈階段。

圖4 30°角輸彈壓力曲線Fig.4 Pressure curve of ammunition ramming at 30°

為方便研究與建模加載，本文將輸彈壓力簡化為線性加載過程，如圖4(b)所示，、分別為首次加載、卸載的結束時間，、分別為首次加載、卸載結束時的壓力值，這4個參數(shù)就是4個輸入不確定參數(shù)。

為驗證輸彈機模型的準確性，將圖4(b)中的簡化油缸壓力曲線加載到ADAMS模型中，輸出彈丸軸向位移隨時間變化的曲線，并與試驗測試結果進行對比，結果如圖5所示。測試結果僅包含彈丸在1 m以內的變化曲線，主要是由于試驗中所用到的激光位移傳感器量程有限；盡管如此，圖5中的彈丸仿真位移曲線與試驗結果在有限范圍內還是非常吻合的，這也恰好驗證了輸彈機ADAMS模型的準確性。

圖5 30°偏角輸彈時的彈丸軸向位移仿真與試驗結果對比Fig.5 Comparison between simulated and test results of projectile axial displacement at 30°

3.2 方法適用性驗證

在輸彈過程的系統(tǒng)隨機參數(shù)中選取12個代表性參數(shù)作為本文的研究重點，其統(tǒng)計信息相對更加容易獲取，對輸彈結果的影響相對也更為顯著。具體的參數(shù)統(tǒng)計信息如表1所示。

表1 輸彈過程中的系統(tǒng)參數(shù)不確定性Tab.1 Uncertainty of system parameters during ammunition ramming

此外，由于輸彈與卡膛過程是火炮發(fā)射的前兩個連續(xù)階段，沒有明確的界面，本文設置彈丸軸向位移達到坡膛前的某一點時輸彈仿真結束，認為此時的彈丸狀態(tài)變量值即輸彈動力學模型的輸出，同時也是后續(xù)卡膛過程的輸入。

考慮到輸彈過程，尤其是慣性輸彈階段彈丸進入內膛后的運動狀態(tài)參數(shù)難以測量，這里選擇與蒙特卡洛方法的計算結果對比來驗證上述方法的有效性。以0°輸彈角為例，仿照圖4中的簡化方法，根據(jù)0°角輸彈壓力測試得到系統(tǒng)輸入的參數(shù)不確定性如表2所示。

表2 0°角輸彈時的輸入?yún)?shù)不確定性Tab.2 Uncertainty of input parameters in ammunition ramming process at 0°

同時考慮這16個系統(tǒng)隨機變量和輸入隨機變量，分別通過MC方法和PDEM計算輸彈結束時的彈丸概率密度函數(shù)。其中，MC方法通過拉丁超立方抽樣方法直接采樣10 000組進行仿真計算，對結果統(tǒng)計得到彈丸響應的概率密度函數(shù)；PDEM則按照3.2節(jié)的步驟，利用概率密度演化的方法采樣500組樣本點，同樣可以計算獲得彈丸響應的概率密度函數(shù)。以彈丸軸、軸方向的位移和自轉角為對象，分別計算其概率密度函數(shù)，結果如圖6～圖8所示。

圖6 擬合概率密度- 彈丸y軸方向位移曲線Fig.6 Projectile displacement in y direction

由圖6、圖7可以看出：0°角輸彈時的彈丸軸方向位移明顯偏大，這是重力作用的結果；但軸、軸方向位移的波動范圍是接近的，其概率密度函數(shù)均為單峰分布，而通過PDEM擬合的概率密度函數(shù)曲線與MC方法統(tǒng)計結果幾乎一致，僅在邊緣以及峰部存在微小誤差。圖8所示為擬合概率密度函數(shù)- 彈丸軸方向偏角曲線。由圖8可見，0°角輸彈時彈丸自轉角的概率密度函數(shù)是雙峰分布，PDEM較為準確地擬合出雙峰趨勢，曲線也基本吻合。

圖7 擬合概率密度- 彈丸z軸方向位移曲線Fig.7 Projectile displacement in z direction

圖8 擬合概率密度- 彈丸x軸方向偏角曲線Fig.8 Spin angle of projectile in x direction

獲得概率密度函數(shù)，不僅可以直觀判斷彈丸響應的具體分布情況，還可以進一步分析其統(tǒng)計特性。根據(jù)概率密度函數(shù)計算彈丸響應的前2階原點矩如下：

(20)

(21)

式中：E(·)表示期望。進而有

=E()

(22)

(23)

式中：為的均值；為的均方差。彈丸響應統(tǒng)計特性結果如表3所示。

表3 彈丸響應統(tǒng)計特性Tab.3 Statistical characteristics of projectile response

概率密度演化方法僅通過500組樣本即可準確地擬合出輸彈響應概率密度函數(shù)的趨勢，無論是單峰還是多峰，這是通過低階矩擬合概率密度函數(shù)的方法無法比擬的，表明其對于輸彈過程高維不確定性分析問題的適用性。另一方面，表3對比了該方法的計算精度，僅需500組樣本點即可獲得誤差不超過5%的均值、均方差，在保證高效率的同時又滿足了工程實際需求。

3.3 輸彈過程的概率密度演化

為更加全面地分析參數(shù)不確定性對輸彈過程系統(tǒng)響應的影響，以0°、30°以及51°輸彈角為例，直接通過概率密度演化方法擬合彈丸響應的概率密度函數(shù)，進行不確定性分析。系統(tǒng)參數(shù)統(tǒng)計特性仍見表1，0°輸彈角時的輸入?yún)?shù)不確定性見表2。表4、表5分別給出了30°、51°輸彈角時的輸入?yún)?shù)不確定性。

表4 30°角輸彈時的輸入?yún)?shù)不確定性Tab.4 Uncertainty of input parameters in ramming process at 30°

表5 51°角輸彈時的輸入?yún)?shù)不確定性Tab.5 Uncertainty of input parameters in ramming process at 51°

影響彈丸卡膛的主要因素是彈丸的位置、姿態(tài)以及卡膛速度，故本節(jié)以彈丸在軸、軸方向的位移、偏角以及軸方向的速度為對象，擬合的概率密度函數(shù)如圖9～圖13所示。

圖9 擬合概率密度- 彈丸y軸方向位移曲線Fig.9 Projectile displacement in y direction

圖10 擬合概率密度- 彈丸z軸方向位移曲線Fig.10 Projectile displacement in z direction

圖11 擬合概率密度- 彈丸y軸方向偏角曲線Fig.11 Deflection angle of projectile in y direction

圖12 擬合概率密度- 彈丸z軸方向偏角曲線Fig.12 Deflection Angle of projectile in z direction

圖13 擬合概率密度- 彈丸x軸方向速度曲線Fig.13 Velocity of projectile in x direction

與16個不確定性參數(shù)相似的是，對應3種輸彈角的彈丸響應分布都是對稱的，都屬于單峰類正態(tài)分布。圖9～圖12表明，隨著輸彈角的增大，彈丸位置、姿態(tài)的波動愈大，彈丸軸線與身管軸線的偏差增加，卡膛一致性變差。圖13表明，彈丸在軸方向的速度最低接近2 m/s，與卡膛速度要求的3 m/s相差較大。

表6所示為由(23)式計算的不同輸彈角下的彈丸響應均方差。對比表6中的數(shù)據(jù)可以發(fā)現(xiàn)：隨著輸彈角的增大，響應的均方差幾乎成倍遞增；在給定的參數(shù)不確定性條件下，軸方向的位移均方差明顯大于軸方向，最大可達62 mm；偏角均方差相差無幾，均不超過0.08°，卡膛速度均方差最大可達0.4 m/s。

表6 不同工況下的響應散布Tab.6 Dispersion of response under different working conditions

考慮到系統(tǒng)參數(shù)的誤差幾乎不變，為減小彈丸響應的波動，提高卡膛一致性，避免卡膛失敗，需要進一步減小輸入載荷的誤差。

4 結論

本文通過分析某火炮輸彈機的結構與原理，建立了輸彈機ADAMS仿真模型，基于概率密度演化理論推導了輸彈過程概率密度演化方程，并通過有限差分法對其進行求解；分析了不同工況下參數(shù)不確定性對輸彈響應的影響。得到如下主要結論：

1)通過仿真結果與試驗結果的對比，驗證了該火炮輸彈機動力學仿真模型的準確性。

2)在相同參數(shù)條件下，通過PDEM與MC方法的對比，驗證了本文方法對涉及高維不確定性參數(shù)的輸彈過程不確定性分析的適用性。

3)基于不同輸彈角的載荷條件擬合彈丸響應的概率密度函數(shù)，輸彈角越大，彈丸響應的波動越明顯，卡膛一致性越差。

4)在本文給定的參數(shù)不確定性條件下，輸彈結束時彈丸的位移均方差最大可達60 mm，卡膛速度均方差最大可達0.4 m/s。