煙幕初始云團爆炸分散模型建立及計算方法

陳 浩，高欣寶，李天鵬，張 倩，陳玉丹，楊 洋

(1.陸軍工程大學石家莊校區(qū)， 石家莊 050003； 2.陸軍北京軍代局駐763廠軍代室， 太原 030000)

為有效對抗制導武器的精確打擊，提高部隊戰(zhàn)場生存能力，發(fā)展與之對抗的干擾理論與技術，受到世界各國軍隊的重視[1]。煙幕無源干擾憑借成本低、效費比高的優(yōu)點已經成為干擾領域發(fā)展的重要研究方向。其中，分析煙幕初始云團的運動規(guī)律，計算煙幕初始云團的參數(shù)，是煙幕作戰(zhàn)效能評估研究的熱點之一[2-3]。爆炸分散型煙幕經中心擴爆管爆轟驅動，向四周迅速膨脹，當殼體破裂瞬間，形成了由爆轟產物與干擾劑小顆粒構成的氣/固混合物云團，將其稱為煙幕初始云團。從形成到煙幕粒子速度衰減為零的過程，稱為煙幕初始云團的爆炸分散過程。將煙幕初始云團近似看成不斷擴張的球體，其最大半徑用于表征煙幕遮蔽面積，作為煙幕評估中的一項重要指標[4-6]。

為了探索煙幕初始云團運動規(guī)律，提高爆炸分散的理論模型精度和適用性，近年來，學者們做了大量工作。Chen等[7]研究并建立了真空環(huán)境中煙幕的膨脹模型；Zhu等[8]分析了赤磷煙幕初始云團膨脹過程，指出空氣阻力對煙幕膨脹影響較大，利用常數(shù)變異法建立了基于膨脹力和空氣阻力作用下的煙幕初始云團膨脹模型；Xu等[9]考慮忽略等熵過程中氣動阻力作用的情況，建立煙幕初始云團膨脹模型，采用四階龍格—庫塔法對模型進行了計算。

上述研究針對煙幕初始云團爆炸分散形成過程建模與計算的不同方面進行了創(chuàng)新與改進，但仍存在以下問題：(1)根據(jù)數(shù)值計算結果分析，等熵膨脹階段的氣動阻力對模型精度有一定影響，建立理論模型時不能忽略；(2)理論模型采用四階龍格—庫塔法求解數(shù)值解的效率低，收斂性和平滑性存在不足。針對以上問題，為進一步提高模型精度及計算效率，本文對煙幕初始云團爆炸分散的理論模型和模型計算方法作進一步改進，將爆炸分散過程分為等熵膨脹和自由擴散兩個階段，考慮等熵膨脹階段的氣動阻力，推導并建立了煙幕初始云團爆炸分散理論模型，利用歐拉法求解模型的數(shù)值解，并與四階龍格—庫塔法進行對比分析。通過開展野外測試試驗，對改進后的模型和所提算法進行了分析與驗證。

1 理論部分

發(fā)煙裝置起爆后，產生爆轟波，并迅速形成大量高溫高壓氣體(中心黑索金藥柱爆炸形成)向四周膨脹做功，驅動周圍的碳基干擾劑向外擠壓，殼體內部壓力逐漸增大，當超過殼體承受極限后，發(fā)煙裝置發(fā)生爆炸，形成一個近似球體且不斷膨脹的高溫高壓云團。殼體膨脹到破裂過程時間極短，為使研究問題簡化，假設爆轟前后爆轟產物和云團體積無變化，粒子速度為零，即認為爆轟瞬時完成。

爆炸分散過程實際為爆轟產物在空氣介質中流動并驅動碳基干擾劑向外膨脹的過程，煙幕初始云團是爆轟產物與碳基干擾劑形成的混合物。根據(jù)等熵膨脹基本理論，同時，考慮便于計算爆轟產物壓力，假設煙幕初始云團在內部壓力衰減為零之前的過程視為爆轟波在空氣中的絕熱等熵膨脹過程，即爆炸分散過程的時間分界點為爆轟壓力衰減為零時，分為等熵膨脹和自由膨脹兩個階段。

1.1 等熵膨脹階段模型建立

取云團邊緣一個小微元體，其質量為dm，受爆轟產物壓力作用面積為dS(微元體體積遠小于云團體積)，受到的作用力主要是爆轟產物壓力和氣動阻力，則有：

(1)

式(1)中：a為微元體膨脹過程的加速度，可表達為a=d2R/dt2；P為爆轟產物壓強；k′為氣動阻力系數(shù)；ρ′為標準條件下大氣密度；v為微元體的速度。

微元體的狀態(tài)方程可表示為：

p=A(S′)·ρk

(2)

式(2)中，A(S′)為常數(shù)。

為求得云團運動參數(shù)的數(shù)值解，在等熵膨脹假設前提下，得到以上微元體壓力表達式，則可進一步得到微元體爆轟發(fā)生后任意時刻的狀態(tài)與初始狀態(tài)之間的關系如下：

(3)

式(3)中：P0為爆炸初始時刻形成的高溫高壓云團的壓強；ρ0為爆炸初始時刻高溫高壓云團的密度；ρ為爆炸后等熵膨脹階段任意時刻高溫高壓云團的密度；k為等熵指數(shù)。

假設爆炸后形成的云團為不斷膨脹的球體，則云團密度可用體積和質量代替，并代入式(3)中可得：

(4)

式(4)中：m為煙幕初始云團質量；R為煙幕初始云團的膨脹半徑。

將式(4)代入式(1)整理得到如下表達式：

(5)

將等熵膨脹階段模型整理，得到如下微分方程組：

(6)

1.2 自由膨脹階段模型建立

在等熵膨脹階段，爆轟作用力在云團向外膨脹過程中不斷衰減。當壓力等于標準大氣壓時，進入自由擴散階段。這一階段，微元體主要受氣動阻力作用，則：

(7)

由式(7)知，自由膨脹階段，微元體加速度為負值，速度隨時間不斷衰減，直到速度衰減到零，此階段結束。經整理化簡得到煙幕初始云團爆炸分散自由膨脹階段的微分方程組為：

(8)

2 實驗

根據(jù)特種彈藥戰(zhàn)斗部設計原理，建立了一種外形為圓柱體的發(fā)煙裝置裝藥模型，如圖1所示。殼體材料為合金鋼，厚度為2.50 mm，高為127.50 mm。中心裝藥為黑索金藥柱，裝填密度為1.70 g/cm3，裝藥半徑為9.00 mm，雷管位于中心裝藥上端，用于起爆黑索金藥柱。四周為碳基干擾劑，裝填密度為2.30 g/cm3，裝藥半徑為56.50 mm。

為有效獲取煙幕初始云團爆炸分散過程中的圖像數(shù)據(jù)，設計了發(fā)煙裝置野外煙幕測試試驗，圖2為試驗系統(tǒng)的示意圖。試驗主要利用高速攝像機對煙幕爆炸分散過程進行記錄，并通過測距儀測量吊架兩側距離，基本參數(shù)值見表1。

圖2 試驗系統(tǒng)示意圖

型號幀率/(幀/s)分辨率VEO4001 0001 280×800

利用獲取的圖像數(shù)據(jù)，采用圖像邊緣檢測法，獲取吊架兩側的像素點坐標和云團邊界兩側的像素點坐標，通過比例關系求解云團徑向膨脹的最大半徑，其表達式為：

(9)

式(9)中：L為吊架兩側的距離；Δx1為圖像中吊架兩側像素點距離；Δx2為圖像中云團兩側像素點的最大距離。

3 結果與討論

3.1 理論模型計算與分析

3.1.1基本原理

式(6)和式(8)為二階非線性微分方程組，可采用數(shù)值積分法求得數(shù)值解。常用的數(shù)值積分法包括歐拉法和四階龍格—庫塔法等。一般針對低階微分方程，為確保數(shù)值解的光滑性，采用歐拉法進行計算，計算公式如下：

(10)

為進一步計算模型的數(shù)值解，首先需確定初值。根據(jù)中心裝藥的體積和裝藥密度等參數(shù)，根據(jù)Liu等[10]提供的方法，通過計算不同氮當量計算爆壓的公式為：

P′=1.092(ρD∑N)2-5.74

(11)

式(11)中：P′為炸藥爆轟壓力；ρD為炸藥裝藥密度；∑N為炸藥氮當量。

通過式(10)計算得到瞬時爆轟時刻，爆轟壓力P′=27.76 GPa。圓柱殼體破裂的強度為1.20 GPa，因此，可確定求解理論模型爆炸瞬時的初值為：P0=1.20 GPa，R0=56.50 mm。

3.1.2數(shù)值計算與分析

根據(jù)初始參數(shù)，編寫了歐拉法的MATLAB程序，綜合考慮計算效率和精度，時間步長設為1×10-8s，根據(jù)圖3～圖5可知，第一階段發(fā)生在0～4.62×10-5s時間范圍內，由于破殼壓力遠小于中心炸藥爆炸所產生的爆壓，因此，壓力在極短時間達到最大值，之后，隨著爆轟能量的衰減，壓力逐漸衰減到標準大氣壓。煙幕初始云團加速度迅速減小，當膨脹力小于氣動阻力時，加速度成為負值，后緩慢增加；速度先增大后較小，最大半徑膨脹到初始半徑的2～3倍。在4.62×10-4s以后，煙幕初始云團進入爆炸分散的第二階段。在第二階段的短時間內，加速度和速度急劇衰減并趨于0為止。此時，膨脹半徑增大到約3 m。

3.2 實驗結果與討論

采用VEO400高速攝影機進行拍攝，截取0～1 s煙幕初始云團膨脹的圖像，如圖6所示。以1 s這一時刻的圖像為例，當背景容易區(qū)分情況下，采用圖像的二值化，將云團和背景區(qū)分開，通過局部灰度拉伸，獲取灰度突變的邊緣灰度閥值，得到煙幕初始云團邊緣形態(tài)圖像，并求解得到煙幕初始云團最大半徑。

圖3 云團膨脹加速度隨時間變化曲線

圖4 云團膨脹速度隨時間變化曲線

圖5 云團膨脹半徑隨時間變化曲線

圖6 發(fā)煙裝置靜爆試驗圖

通過上述算法處理得到的煙幕初始云團邊緣曲線較為清晰，如圖7所示，利用式(9)分別計算各時間節(jié)點煙幕初始云團最大半徑，得到試驗值。同時，分別計算了改進和未改進模型對應時刻云團最大半徑并進行了比較與分析，如圖8所示，改進模型求解得到的膨脹半徑隨時間變化曲線比未改進模型更接近試驗結果，為進一步證明了改進模型的適用性，計算了二種模型的相對誤差曲線，如圖9所示。誤差隨時間逐步趨于穩(wěn)定，未改進模型誤差在10%左右，改進模型的誤差約為5%，改進模型計算結果的相對誤差較小，相比而言，相對誤差降低了約5%～10%。煙幕初始云團膨脹半徑，在爆炸發(fā)生短時間內，急劇增大。之后，逐漸趨于平緩，其變化規(guī)律與試驗測試得到的變化規(guī)律基本一致。

云團爆炸分散的整個階段，其理論計算值大于試驗測試值，其主要原因在于：爆炸分散的整個階段，爆轟產物的能量，有一部分被殼體的塑性變形及其他因素所消耗，有效驅動的能量低于理論值，因此，基于一定假設的理論模型計算結果要大于試驗結果；試驗測試時，高速攝影機圖像捕獲精度為1×10-3s，第一階段的整個發(fā)生時間小于1×10-3s，捕獲精度存在誤差，因此，試驗測試結果低于理論值。為提高工程計算的適用性，在建立的爆炸分散模型基礎上，需增加修正項，對理論模型作進一步修正，降低初始時刻模型的計算誤差。

圖7 云團邊界形態(tài)圖像

圖8 云團膨脹半徑隨時間變化曲線

圖9 相對誤差隨時間變化曲線

4 結論

通過理論建模、數(shù)值計算與實驗分析相結合的方法，對圓柱體發(fā)煙裝置煙幕初始云團運動規(guī)律進行了分析。結果表明考慮第一階段氣動阻力的改進模型相比未改進模型，計算值誤差較小且更接近試驗值，能夠描述煙幕初始云團膨脹半徑的變化規(guī)律。