基于增強型離散相模型和化學燃燒的紅外誘餌彈建模方法

2020-12-29 02:32:44楊春玲張振東張潼奕聿

航空學報 2020年12期

楊春玲，張振東，張潼奕聿

哈爾濱工業大學電氣工程及自動化學院，哈爾濱 150001

紅外誘餌彈是一種使用率很高的紅外對抗武器。問世至今，它憑借成熟、可靠及價格低廉等優點，被各軍事大國廣泛應用于實戰[1-2]。當紅外目標被鎖定時，會發射紅外誘餌彈來誘騙進攻武器，達到保護自身的目的。因此，建立一種準確的紅外誘餌彈仿真模型，無論是針對紅外制導算法的抗干擾性研究還是紅外誘餌彈的誘騙效果研究，都將具有重大的意義。

目前針對紅外誘餌彈的主流建模方法有基于紋理特征的、基于仿真粒子算法的和基于計算機流體仿真(Computational Fluid Dynamics，CFD)的紅外誘餌彈建模法等。

基于紋理特征的建模方法[3-5]是一種相對較為傳統的方法，它的原理是先通過觀測和總結火焰、煙幕等對象的條紋信息，再利用數理統計的方法建立出相應的經驗模型。由于這種方法屬于早期的經驗公式建模法，缺少機理性的研究，因此在計算機運算能力明顯提升后很快就被其他數值方法取代了。

基于仿真粒子算法的紅外誘餌彈建模方法[6-7]思想如下：首先，假設紅外誘餌彈是由空間中初始位置隨機分布的大量紅外輻射粒子構成的；然后，以每一個粒子為研究對象，通過施加大小和方向都隨機的力驅動粒子做布朗運動；接下來，根據仿真要求添加“仿真風”，使粒子具有某些一致性的運動規律；最后，通過添加這兩種因素仿真得到每個顆粒的運動軌跡，進而仿真得到紅外誘餌彈的模型。這種方法屬于統計學建模方法，它運算量較小且生成圖像速度較快，目前在很多半實物仿真平臺上有廣泛應用[8]。然而，這種方法缺乏機理性研究，也并沒有過多地研究紅外誘餌彈輻射特征的產生機理。

隨著計算機運算能力的不斷增強，基于CFD的紅外誘餌彈仿真方法[9-10]憑借其仿真精度高、機理性研究能力強的優點，已逐漸成為主流研究方法。而基于CFD的仿真方法主要有連續相建模和離散相建模兩種方法。

其中，連續相建模方法將紅外誘餌彈噴射的物質簡化為熱氣流或可參加化學反應的氣體物質(如CO)，然后求解Navier-Stokes(N-S)方程，建立湍流模型，計算得到誘餌彈的輻射特征。這兩種連續相的假設中，簡化為熱氣流計算較為容易，但由于氣體比熱容較小，因此換熱速度較快，模型輻射特征持續時間較短；而簡化為可參加化學反應的氣體物質盡管可以通過燃燒獲得化學能量補償，保持更長的輻射時間，但這種仿真方法將化學反應簡化為均相反應，不符合化學反應動力學過程，因此其輻射光譜準確性較差。

相較連續相建模方法，離散相建模方法是近些年發展起來的一種新型建模手段。它首先以顆粒為研究目標，通過建立Lagrange坐標系分析流場中顆粒的受力情況，計算分子的運動軌跡；然后借助計算機仿真能力，通過追蹤大量顆粒建立宏觀物體的運動、輻射等模型。

離散相模型對于由大量顆粒(如煤、沙塵等)組成的宏觀現象研究效果較好，而紅外誘餌彈也有著上述相同的特征，它是一種由藥劑顆粒運動和燃燒產生輻射特征的紅外目標。因此，離散相建模方法將更加貼近真實物理過程，具有更高的建模精度。張振東等[11-12]建立了一種基于離散相模型(Discrete Phase Model，DPM)的紅外誘餌彈輻射特征模型，首先以顆粒作為研究對象，通過求解力學平衡方程和N-S方程計算每個顆粒的受力情況和運動軌跡；然后加入離散坐標(Discrete Ordinates,DO)輻射模型建立紅外誘餌彈的輻射特征模型。在此基礎上，張振東等[13]又建立了增強型的DPM模型，通過編寫用戶自定義編程文件(User Define Files，UDF)的方法引入了多種附加力，建立了復合顆粒噴口，并利用擴展型“等效黑體分子”輻射模型建立了紅外誘餌彈的輻射特征模型。

雖然這種離散相模型確實具有較高的仿真精度，但以上研究[11-13]并沒有考慮顆粒燃燒對紅外誘餌彈輻射特征產生的影響。而通過分析紅外誘餌彈的工作機理[14-18]可以發現，誘餌彈的輻射特征主要由兩部分構成：① 由氮氧化合物和碳氧化合物構成的高溫氣體；② 由燃燒的紅外誘餌藥劑(如Mg、聚四氟乙烯(PTFE)和Viton混合物(MTV)材料)構成的高溫顆粒。由此可見，僅以氣體相或僅以固體相建模，而不考慮化學燃燒對紅外誘餌彈的影響，得到的結論無疑會在一定程度上存在誤差。

綜上，本文將針對化學組分燃燒對輻射特征產生的影響做進一步研究。首先，建立一種層流場中的MTV燃燒數值仿真模型，該模型以MTV顆粒為研究對象，通過分別求解流體力學基本方程組和化學組分守恒定律，先計算MTV顆粒的空間運動軌跡、與周圍環境在質量和能量上的交換情況，再引入增強型DPM模型，計算得到顆粒的運動狀態，并借助CFD仿真軟件計算不同時刻的MTV顆粒燃燒時的物理狀態；然后，建立擴展型“等效分子黑體”模型，計算紅外誘餌彈的輻射特征；最后，設計實驗進行驗證。

1 數學模型

實戰中，紅外誘餌彈會噴射大量燃燒顆粒在作戰空域形成紅外干擾。因此在理論研究紅外誘餌彈的數學模型時，需要分別從氣體相、固體離散相和化學燃燒相的角度進行建模研究。

1.1 氣體相的控制方程

建立層流低馬赫數下的氣相模型，其質量守恒定律、化學組分守恒定律、動量守恒定律和能量守恒定律可以變形為

(1)

(2)

(3)

(4)

(5)

(6)

假定統一Lewis數，化學組分的擴散系數可以根據熱系數計算得到。基于Stokes假說的牛頓流體假設，黏性張力可以通過計算得到：

(7)

式中：μg為混合氣體的動態速度；U為流體速度。

氣體相均勻反應釋放的熱量可以通過式(8)計算得到：

(8)

式中：hα為第α相的反應焓變量。

理想氣體的狀態方程用于使方程組收斂。組分傳輸系數、反應速率和混合物的熱力學關系可以通過Blanquart等[19]、Cai和Pitsch[20]提出的方法計算得到。

1.2 固體相的控制方程

對于固體相而言，能量方程可以通過式(9)進行求解：

(9)

式中：ρs、cp,s、Ts和λs分別為固體相的密度、定壓比熱、溫度和導熱系數。混合物中，Mg作為分散相球形顆粒散布在PTFE連續相基體中。二元混合物等效導熱系數λs可以通過Maxwell模型計算得到：

(10)

式中：λm為PTFE的導熱系數，λm=0.244 W·m-1·K-1；λd為Mg的導熱系數，λd=156 W·m-1·K-1；Vd為Mg的體積分數。

MTV顆粒的密度主要取決于MTV顆粒在高溫情況下的揮發速率和顆粒表面的氧氣分布程度。MTV煙火劑的密度為1.79 g/cm3，Mg、PTFE和Vtion的質量分數分別為60wt%、35wt%和5wt%。

1.3 DPM模型的數學模型

DPM模型是通過對Lagrangian參考系下顆粒的運動方程求積分計算運動軌跡的。考慮顆粒的慣性與受力平衡，分散相粒子運動方程(以直角坐標系內x方向為例)為

(11)

(12)

(13)

(14)

式中：μ為連續相黏度；CD為顆粒的比熱容；Re為雷諾數；dp為顆粒直徑；a1和a2為系數，a1=0.364 4，a2=98.33。

在固體相中，質量傳輸和組分傳輸被忽略了。因此，對于固體相而言，只需通過式(9)計算它的能量方程。

1.4 氣-固交界面處的控制方程

根據組分、質量和能量源與氣體和固體通量之間的平衡，可以得到組分、質量和能量交界面方程為

(15)

(16)

(17)

圖1 化學反應時組分微觀變化示意圖Fig.1 Schematic diagram of microscopic changes of components in chemical reactions

(18)

(19)

式中：ε為發射率，ε=0.73；σ為Stefan-Boltzmann常數；TW為壁面溫度。

固體相表面非均相反應決定的每一組分的產生和消失速率可以通過如式(20)和式(21)所示的化學反應進行計算。由于交界面處的非均勻性反應，Waite等[21]提出了一種用于計算每一項產生和消耗速率的化學反應機理：

(m-2n)Mg(g)+nC(s)+qan

(20)

2CO2+qae

(21)

式中：qan為厭氧環境中化學反應焓；qae為空氣中化學反應焓。

化學組分在燃燒時的燃燒特性如文獻[15-17]所示，紅外煙幕的光照強度Lλ(W·sr-1·cm-2)由化學組分的消耗速率?m/?t(g·s-1·cm-2)和光譜效率Eλ(J·g-1·sr-1)共同決定，他們之間的數值關系為

Lλ=Eλ?m/?t

(22)

光譜效率Eλ與整體化學反應的焓變(ΔrH=qan+qae)、發射率ε和化學反應生成物的溫度之間的數值關系為

Eλ=(4π)-1ΔrHF(λ,T)

(23)

式中：F(λ,T)為不同溫度下光譜輻射率，可根據普朗克公式計算；λ為光的波長。

Koch等[22-25]給出了MTV煙火劑中各組分的物理特性和化學特性，如表1所示。

表1 MTV煙火劑中各組分的物理特性[22-25]

2 基于動態計算網格的紅外誘餌彈仿真模型

在分別從氣體相、固體離散相和化學燃燒時的氣固耦合等方面對紅外誘餌彈進行理論建模后，進一步建立動計算網格研究紅外誘餌彈在旋轉時的紅外特征。

2.1 “復合顆粒噴口”邊界條件設計

在連續相入口處，流場的速度和湍流強度等邊界條件被分別設置為常數。由于離散相的空間分布存在明顯的隨機性，通過編寫用戶自定義函數(UDF)設計了復合顆粒邊界條件。復合顆粒噴口邊界條件是將顆粒按照直徑分成n組，并保證每一組的顆粒均滿足高斯分布。因此，當n的數值很大時，數值模型就將能夠更加準確地模擬實驗模型。為了計算方便，令n′=6，直徑分別為0.50、0.75、1.00、1.25、1.50、1.75 mm，每一組均占顆粒總量的16.7%，噴口的截面圖如圖2所示。

在出口和遠場處，為防止由氣壓差造成的空氣回流，出口設置為數值是1 atm的壓力出口。同時，當顆粒經過出口處后，顆粒將會被從計算域中刪除。

彈體表面的壁面條件選擇無滑移標準壁面函數，同時假設當顆粒碰撞彈體時會產生反彈。為研究標準壁面函數中顆粒反彈因子對仿真結果的影響，分別設計了反彈參數為0.5、0.7和1.0的3個實驗，結果如圖3所示。結果表明，反彈因子對顆粒在彈體表面處的分布影響并不明顯，因此反彈參數選擇1.0。

圖2 噴口處個質量顆粒空間分布示意圖Fig.2 Spatial distribution diagram of mass particles at nozzle

圖3 壁面反彈因素對顆粒空間部分的影響Fig.3 Influence of wall rebound on spatial part of particles

2.2 動網格計算域劃分

在借助有限元仿真工具對增強型DPM模型的數值模型進行仿真驗證時，需要根據研究對象的結構特征和運動特征等設計適當的流場計算網格。由于建立的模型考慮了彈體多自由度運動對仿真結果的影響，因此設計了一種基于動態網格的流場計算網格。大量仿真實驗表明，動態網格的網格形式、參數和結構等因素將直接影響仿真的時間長度和最終的仿真結果。所以，在設計多組對比實驗后，最終建立了一套恰當的流場計算網格。

圖4(a)為“∏P∏-50”紅外誘餌彈的彈體結構。可以看出，噴口平均分布在彈體的表面，使紅外誘餌彈的彈體表面呈現一種“玉米棒”形。根據這種彈體結構模型建立如圖4(b)所示的3D模型。

圖4 紅外誘餌彈結構Fig.4 Structure of infrared decoy bomb

由于實戰中誘餌彈出膛瞬間附帶的切向力會導致彈體旋轉，從而影響整個紅外誘餌彈在空間中的分布。因此，為了能夠仿真誘餌彈旋轉時紅外誘餌彈的空間分布和物理特征，將流場分割為如圖5所示的形式，即靠近彈體的流場區域被設計為半球形，其近場網格均為結構化網格，如圖5(b)所示；遠場部分的網格同樣是由尺寸更大的結構化網格構成，如圖5(c)所示。

在將流場整體劃分為遠近兩個子計算域后，下一步即分別研究它們的網格形式和網格數量等參數的劃分方案。在設計網格時常常需要面對一些兩難的選擇，例如：盡管網格數量的增加可以使仿真結果更加接近真實結果，卻會帶來計算時間成倍增長的問題；結構化網格雖然會使計算時間縮短，使計算結果呈現更好的收斂性，但同時也為網格設計工作帶來很大的不便。

針對這些網格參數的選擇設計了對比實驗，計算網格參數選擇對比試驗結果如表2所示。對比實驗的判斷依據綜合考慮了網格質量和尾焰仿真結果溫度變化率。其中，計算溫度變化率的方法是首先計算相同彈體和流場尺寸時，使用結構化靜態網格仿真得到的最高溫度T0，然后分別計算各個網格劃分方案中仿真結果的溫度Tn′，最后計算1-(Tn′/T0)。

圖5 流場的空間分布和計算網格Fig.5 Spatial distribution of flow field and computational mesh

2.3 數值仿真結果及分析

在確定網格形式和網格參數后，使用該計算網格進行數值仿真，分別得到了彈體的組分分布圖和彈體旋轉時的溫度云圖，如圖6和圖7所示。

從圖6中可以看出組分分布具有如下共同特點：MTV藥劑剛被噴射出來時，C含量和Mg蒸汽含量比較高；高濃度MTV藥劑的劇烈化學反應迅速消耗了氧氣，在彈體近處形成厭氧區；此時的化學反應過程如式(20)所示。

表2 計算網格參數選擇對比試驗結果Table 2 Test results comparison of parameter selection of computational grid

圖6 紅外誘餌彈尾焰組分分布Fig.6 Component distribution at tail flame of IR decoy

圖7 紅外誘餌彈溫度分布云圖的CFD仿真結果Fig.7 CFD simulation results of infrared decoy temperature distribution nephograms

同時，根據圖6也可以看出，風速越快，單位時間內燃燒區域會加入越多的氧氣，因此彈體周圍厭氧區域會縮小，如圖6(a)～圖6(d)和圖6(g)～圖6(j)所示；盡管厭氧區域被壓縮了，但由于化學反應物總量并沒有發生變化，所以在有氧區的化學反應更加劇烈了，因此最終的化學反應產物CO2和MgO生成速度更快，擴散區域更大，如圖6(e)、圖6(f)、圖6(k)和圖6(l)所示。

總的來講，風速越快，組分云圖橢度e(即長寬比)越大；反之，圖像越趨近于圓形。

隨著C顆粒和氣態Mg的擴散，它們在經過了“厭氧、有氧交界面”后逐漸進入“有氧區域”。化學過程中的氧化物也逐漸由F變為O。因此在交界面附近，化學反應的中間產物CO和MgF2含量較高，這一特點如圖6(c)、圖6(d)、圖6(i)和圖6(j)所示；而在富氧區中，CO進一步與O2反應，全部變成CO2。Mg的所有中間產物與O2反應，最終全變為MgO，這一特點如圖6(e)、圖6(f)、圖6(k)和圖6(l)所示。

圖7為當誘餌彈的飛行速度為1Ma、彈體旋轉速度為1 rad/s時加入燃燒相和未加燃燒相的對比溫度云圖。從圖7可以看出，燃燒相的引入主要導致溫度云圖出現了一些明顯的變化。首先，燃燒相使計算域內高溫區域擴大，并且溫度的最高值從約3 500 K提高至約5 000 K。這是由于如果不考慮化學燃燒，DPM模型會將噴射顆粒等效為發熱顆粒，這時流場的能量來源主要是與顆粒熱交換；而燃燒相的加入會使顆粒釋放大量的化學能，并使溫度進一步上升。因此，溫度云圖中心溫度必然升高，也必將導致作用空域面積增大。

另外，燃燒相的引入使流場中各局部區域因為氧含量不同導致化學反應熱略有不同，因此化學反應焓也不同。體現在溫度云圖上的變化主要就是波動性增強。

3 基于“復合顆粒等效分子黑體模型”的紅外輻射特征計算方法

由于建立的仿真模型將用于紅外實景仿真系統中，而該系統顯示的內容均為紅外目標的紅外云圖。因此，進一步建立非均勻輻射模型計算紅外誘餌彈的輻射特征。

紅外煙幕內部燃燒顆粒的空間分布具有非均勻性，這將導致煙幕中各處紅外光的透射率存在差異。所以空間中的紅外煙幕事實上是一種非均勻性半透明紅外輻射體。針對這種輻射體，如直接使用Planck定律計算其輻射強度，勢必會帶來計算上的復雜化，并影響計算精度。

根據上述原因，選用張文華等[26]提出的“等效分子黑體”模型。由于這種模型僅適用于單一直徑顆粒的離散相仿真，無法滿足復合直徑顆粒的計算場景，本文在原有模型的基礎上進行了復合顆粒噴口型擴展，從紅外煙幕計算紅外煙幕的輻射特征，并成功獲得紅外煙幕的紅外分布云圖。然而，張文華等[26]規定了所有顆粒的物理特性均相同，因此為了滿足計算需求，擴展并建立了“復合顆粒等效分子黑體”算法。

將流場在空間上被劃分為大量邊長為D的“微元”，如圖8所示，每個微元中都包含不確定數量的輻射顆粒。流場中每一個輻射體都存在它對應的“等效輻射面積”，例如可以假設某一個顆粒的輻射面積為A，某一個“微元”的輻射面積為S。

圖8 微元結構Fig.8 Structure of micro unit

基于這種方法的“遮擋原理”，每個顆粒對應的輻射面積均被假設為非透明體。因此，當空間中的多個顆粒互相發生遮擋(如圖9所示的3個顆粒發生遮擋)時，它們對應的“等效輻射面積”總和Athree=A1+A2+A3。

圖9 發生遮擋時的等效輻射面圖Fig.9 Equivalent radiation surface map in occlusion

由于增強型DPM將顆粒分為了多組，每一組中的所有顆粒均具有完全相同的物理特征(ρn,Tn,pn)，因此第n組顆粒對應的微元的密度就可以被表示為

(24)

式中：Mn為該微元中第n組顆粒的數量。

根據“復合顆粒等效分子黑體”法的“遮擋原理”，由第n組顆粒產生的等效黑體輻射面積總和可以被表示為

(25)

(26)

在計算了一個“微元”中第n組顆粒產生的等效輻射面積后，進一步計算紅外煙幕的整體輻射特征。首先，整個計算域被劃分成I×J×K個微元，如圖10(a)所示，其中每一個微元的邊長仍為D；然后，取其中一列作為研究對象，如圖10(b)所示，此列中第k個微元的等效輻射面積可以根據式(27)計算求得：

(27)

根據式(27)，第k個微元的輻射面積通量可以表示為

(28)

式中：lk為第k個微元的輻射平面邊長；d為等效輻射體邊長。

根據式(28)，一列微元輻射量可以通過式(29)計算求得：

(29)

(30)

圖10 擴展型“等效分子黑體”模型結構示意圖Fig.10 Structural schematic diagram of extended “equivalent molecular blackbody” model

因此，紅外煙幕的紅外特征就可以通過式(31)計算求得：

(31)

彈體參數、計算域條件和初始化條件如下：彈體尺寸為?50 mm×300 mm，彈體飛行速度為1.2Ma，海拔高度為5 000 m，組分為Mg和Al。在建立了擴展性“等效分子黑體”模型后，根據彈體參數仿真得到了4種不同彈體旋轉速度對應的紅外誘餌彈的紅外云圖，如圖11所示。

圖12所示為轉速為0時不同時刻誘餌彈輻射特征光譜圖。可知，4 s時誘餌彈已經處于穩定燃燒狀態，燃燒可以一直持續至14 s。

另外，誘餌彈分別在波數約為1 000 cm-1和2 300 cm-1處輻射特征較為明顯，這與林長津[27]的實驗結果變化趨勢相符。

4 實驗驗證

為驗證所提模型的仿真精度，設計了兩組實驗，從圖形特征角度驗證紅外誘餌彈的仿真精度。

Labonté和Deck[28]提出了一種徑向基(RBF)神經網絡結構以及12種圖形特征用于辨識同一幅圖像中的紅外目標與紅外誘餌彈。仿真結果表明其中6種圖形特性的辨識率達到94%。因此，分別使用這6種圖形特征對比實驗圖像與仿真圖像的相似度。這些參數的定義詳見表3。

設計兩組對比實驗：實驗1的實驗環境為多云，高空風速36 m/s，誘餌彈飛行速度0.97Ma；實驗2的實驗環境為晴天，高空風速19 m/s，誘餌彈飛行速度0.87Ma。

為驗證燃燒對紅外圖形特征的影響，分別針對實驗1和實驗2設計了兩組對比實驗。實驗結果如表4所示。可以發現由于紅外藥劑的燃燒產生大量高溫燃氣，這些高溫燃氣將直接增大紅外目標在整體圖像中的面積，這直接影響了離心率e和半徑R兩個特征，從而明顯提高了紅外模型的仿真精度。

針對這種兩組實驗條件的仿真結果如圖13所示。可以看出，風速對于誘餌彈成像是有影響的。風速越快時，彈體尾焰拖尾越長，這是由于更快的風速在單位時間內引入了更多的氧氣，導致化學反應的厭氧區被壓縮；同時更快的風速也會促使有氧區與周邊熱交換速度更快。

圖11 不同初始條件下的紅外誘餌彈仿真結果Fig.11 Infrared decoy simulation results under different initial conditions

圖12 MTV煙火劑顆粒在不同時刻燃燒的光譜變化情況Fig.12 Spectral changes of MTV pyrotechnic particles during combustion at different times

在對比風速對仿真結果的影響后，繼續跟進實驗1和實驗2的工況，對比研究了化學燃燒相加入前后對仿真結果的影響。對比實驗分別選取起燃時刻和穩定燃燒時刻作為研究對象，結果如圖14所示。

通過觀察對比實驗結果，可以發現如下3個特點：

1) 在起燃時刻，未加入燃燒相的誘餌彈圖像特征更為明顯。這是由于未加入燃燒相時，DPM模型將顆粒初始狀態假設為高溫顆粒。這些高溫顆粒在與周圍流場換熱后逐漸冷卻，最終失去紅外特征；而考慮燃燒相的DPM模型將顆粒的初始狀態假設為常溫狀態，通過劇烈的化學反應釋放化學能，從而產生紅外特征。這就導致起燃時刻未加入燃燒相的仿真結果更為明顯，這也更符合真實物理過程。

2) 穩定時，加入燃燒相的仿真結果紅外特征更為明顯。這時由于燃燒相對附加化學能，導致反應焓釋放更多熱量。因此其仿真結果中，反應持續時間更長，能量更高。

3) 整個仿真過程中，加入燃燒相后，紅外目標邊緣更加不規則。這是因為仿真流場中出現了大量未完全反應的中間產物，這些中間產物對氧