鄒濤，王超哲，童中翔＊，賈林通，童奇

空軍工程大學(xué) 航空航天工程學(xué)院，西安 710038

箔片型紅外面源誘餌擴(kuò)散規(guī)律

鄒濤，王超哲，童中翔＊，賈林通，童奇

空軍工程大學(xué) 航空航天工程學(xué)院，西安 710038

目前紅外面源誘餌已成為對(duì)抗紅外成像制導(dǎo)導(dǎo)彈的重要方法。鑒于此，以箔片型紅外面源誘餌為研究對(duì)象，利用計(jì)算流體力學(xué)（CFD）方法得到單個(gè)箔片和兩平行箔片在不同迎角下的氣動(dòng)力系數(shù)。通過(guò)對(duì)箔片空中運(yùn)動(dòng)分析，特別是考慮到燃燒力作用的特殊性，將箔片的運(yùn)動(dòng)分為兩個(gè)階段：起燃階段和完全燃燒階段，并建立了兩個(gè)階段的運(yùn)動(dòng)模型。在單個(gè)箔片運(yùn)動(dòng)模型的基礎(chǔ)上，設(shè)置上千個(gè)箔片初始姿態(tài)和旋轉(zhuǎn)角速度的概率分布情況，同時(shí)對(duì)上千個(gè)箔片求解運(yùn)動(dòng)方程，即可得到整個(gè)紅外面源誘餌的空間運(yùn)動(dòng)擴(kuò)散規(guī)律。仿真結(jié)果表明：在高速運(yùn)動(dòng)平臺(tái)上垂直向上發(fā)射紅外面源誘餌，其擴(kuò)散形狀近似為前部稀疏后部密集的錐形云團(tuán)；仿真得到的面源誘餌擴(kuò)散尺寸和運(yùn)動(dòng)位置與實(shí)測(cè)數(shù)據(jù)吻合度較高，特別是面源誘餌在x軸方向擴(kuò)散尺寸誤差小于5%。

箔片；紅外面源誘餌；計(jì)算流體力學(xué)；氣動(dòng)力系數(shù)；燃燒力；擴(kuò)散規(guī)律

為進(jìn)一步提高抗干擾能力，紅外導(dǎo)彈將向紅外成像導(dǎo)彈發(fā)展，后者可將目標(biāo)紅外圖像特征與光譜特征、運(yùn)動(dòng)特征綜合起來(lái)進(jìn)行目標(biāo)識(shí)別，具有非常強(qiáng)的抗干擾能力，使絕大多數(shù)紅外干擾裝置失去作用。紅外成像導(dǎo)彈及其他紅外探測(cè)設(shè)備對(duì)飛機(jī)的生存力形成了巨大而緊迫的威脅，為此面源紅外誘餌應(yīng)運(yùn)而生［1－5］。

趙非玉等［6］通過(guò)理論推導(dǎo)，研究了面源紅外誘餌的運(yùn)動(dòng)特性和三維模型建模方法，基于VC軟件實(shí)現(xiàn)了面源紅外誘餌的運(yùn)動(dòng)學(xué)仿真，利用OpenGL技術(shù)實(shí)現(xiàn)了面源紅外誘餌三維模型和動(dòng)態(tài)紅外圖像生成。付曉紅和樊秋林［7］在面源紅外誘餌的建模與仿真研究中，分別建立了面源紅外誘餌的運(yùn)動(dòng)模型，其中運(yùn)動(dòng)模型包括運(yùn)動(dòng)方程模型和擴(kuò)散模型，具有一定的參考價(jià)值。陳乃光［8］通過(guò)引入勻降速度，推導(dǎo)了誘餌水平和垂直方向速度、位移的計(jì)算公式，結(jié)合圖象分析研究了面源誘餌箔片的運(yùn)動(dòng)擴(kuò)散規(guī)律。黃蓓等［9］通過(guò)試驗(yàn)研究了光電干擾箔片拋撒與云團(tuán)擴(kuò)散問(wèn)題，將過(guò)程劃分為出艙、分離和擴(kuò)散沉降三個(gè)階段進(jìn)行研究，分別提取單個(gè)箔片運(yùn)動(dòng)軌跡和散布云團(tuán)的結(jié)構(gòu)進(jìn)行了統(tǒng)計(jì)分析。

目前國(guó)內(nèi)對(duì)機(jī)載面源紅外誘餌擴(kuò)散規(guī)律的研究才剛剛起步，雖然已經(jīng)成為研究熱點(diǎn)，但和國(guó)外相比還存在不小差距，比如對(duì)機(jī)載面源紅外誘餌擴(kuò)散規(guī)律的研究不夠深入，既沒(méi)有有效的試驗(yàn)測(cè)試，又缺乏可信通用的仿真模型，特別是高空高速條件下面源誘餌運(yùn)動(dòng)特性的研究較少。

針對(duì)這一情況，本文首先，利用計(jì)算流體力學(xué)（CFD）軟件計(jì)算單個(gè)箔片和兩平行箔片在不同迎角下的氣動(dòng)力系數(shù)；然后，通過(guò)對(duì)箔片空中的運(yùn)動(dòng)分析，建立單個(gè)箔片運(yùn)動(dòng)模型；最后設(shè)置上千個(gè)箔片初始姿態(tài)和旋轉(zhuǎn)角速度的概率分布情況，同時(shí)對(duì)上千個(gè)箔片求解運(yùn)動(dòng)方程，得出整個(gè)面源誘餌的仿真結(jié)果。

與點(diǎn)源紅外誘餌不同，面源紅外誘餌發(fā)射出去后，有效載荷擴(kuò)散開(kāi)來(lái)并燃燒，形成具有一定空間體積（形狀）的輻射團(tuán)，因此，每個(gè)載荷單元（即每一個(gè)箔片）的運(yùn)動(dòng)特性對(duì)整個(gè)面源誘餌的擴(kuò)散形狀和運(yùn)動(dòng)特性有決定性影響。面源誘餌箔片是輕、薄、小的物體，其在大氣中的運(yùn)動(dòng)擴(kuò)散受外界條件的影響較大，并且?guī)в幸欢ǔ潭鹊碾S機(jī)性，建模難度較大［10］。鑒于此，本文研究面源誘餌的擴(kuò)散規(guī)律，是在以下幾個(gè)前提假設(shè)下進(jìn)行的：

1）空氣相對(duì)于地面靜止，即沒(méi)有風(fēng)的影響。這樣大氣環(huán)境按國(guó)際標(biāo)準(zhǔn)大氣的參數(shù)來(lái)設(shè)置。

2）箔片質(zhì)量分布均勻，且為剛體，不發(fā)生形變。面源誘餌的箔片是在金屬材質(zhì)的基底上涂敷可遇氧自燃的反應(yīng)物質(zhì)制成的，具有一定的強(qiáng)度和韌度。由于箔片在空中的飛行是自由散落的，且箔片的尺寸較?。ㄖ睆綆桌迕祝饔迷诓系臍鈩?dòng)力通常只能使箔片旋轉(zhuǎn)，而很難使其扭曲變形，因此研究箔片運(yùn)動(dòng)時(shí)可以將其看做剛體。

3）箔片空中飛行時(shí)，不考慮相互之間碰撞的作用效果。這是由于除了面源誘餌剛發(fā)射出的極短時(shí)間，在箔片散開(kāi)的大多數(shù)時(shí)間內(nèi)，箔片之間的間隔均大于箔片自身的尺寸?？梢哉J(rèn)為各個(gè)箔片的運(yùn)動(dòng)是相互獨(dú)立的。

1 箔片氣動(dòng)力數(shù)值計(jì)算

對(duì)面源誘餌運(yùn)動(dòng)擴(kuò)散起關(guān)鍵作用的力是箔片所受的空氣動(dòng)力，研究物體的氣動(dòng)力主要有試驗(yàn)研究、理論分析和數(shù)值模擬3種方法。試驗(yàn)研究真實(shí)可信，但耗時(shí)長(zhǎng)，成本高；理論分析可以得到一些普遍性的規(guī)律，但對(duì)于一些復(fù)雜問(wèn)題，理論分析無(wú)法解決；CFD方法利用計(jì)算機(jī)進(jìn)行數(shù)值模擬，速度快、費(fèi)用低并且用途廣泛，在流場(chǎng)分析中發(fā)揮著越來(lái)越重要的作用［11］。本文研究的箔片形狀薄小，氣流條件復(fù)雜多變，其氣動(dòng)力適宜采用CFD方法進(jìn)行數(shù)值計(jì)算［12］。

1．1 單個(gè)箔片氣動(dòng)力分析

1．1．1 幾何模型與網(wǎng)格

所研究的面源誘餌，其有效紅外載荷是上千片涂敷了可燃物質(zhì)的箔片，單個(gè)箔片呈圓形，直徑為2．54cm、厚度為0．02cm［13－14］。箔片是中心軸對(duì)稱的，因此具有無(wú)數(shù)個(gè)對(duì)稱面，定義箔片速度方向與箔片中軸所成的面為氣流對(duì)稱面，氣流對(duì)稱面與箔片的交線即為箔片的弦線，如圖1所示。

圖1 箔片幾何形狀Fig．1 Foil geometry

這樣箔片所受的氣動(dòng)力始終在氣流對(duì)稱面內(nèi)，可以分解為沿相對(duì)氣流方向（即與速度方向相反）的氣動(dòng)阻力D和垂直相對(duì)氣流方向的氣動(dòng)升力L。速度方向與箔片弦線的夾角定義為箔片迎角α，相對(duì)氣流從箔片下方流過(guò)時(shí)，α為正；從箔片上方流過(guò)時(shí)，α為負(fù)。則箔片的升力與阻力的大小僅與動(dòng)壓和迎角有關(guān)，即

式中：ρ為箔片周圍空氣的密度；V為箔片的速度大?。籗為箔片圓面的面積；CD為阻力系數(shù)；CL為升力系數(shù)；CD和CL的大小由迎角α決定。

利用ICEM建立箔片CFD計(jì)算模型的三維幾何區(qū)域，如圖2所示。

圖2 箔片計(jì)算區(qū)域Fig．2 Foil computation region

在來(lái)流方向上，箔片之前的流場(chǎng)對(duì)箔片周圍流場(chǎng)影響相對(duì)較小，因此為了提高網(wǎng)格質(zhì)量，又兼顧計(jì)算性能，計(jì)算域的前端距箔片中心為10cm，后端距箔片中心為16cm，左右兩端距箔片中心為13 cm。對(duì)計(jì)算域進(jìn)行網(wǎng)格劃分，生成四邊形結(jié)構(gòu)化網(wǎng)格，為提高計(jì)算精度，箔片附近進(jìn)行網(wǎng)格加密處理。箔片附近網(wǎng)格及其放大網(wǎng)格如圖3所示。計(jì)算區(qū)域劃分網(wǎng)格數(shù)為892 652，節(jié)點(diǎn)數(shù)為920 502。

圖3 箔片計(jì)算域網(wǎng)格劃分Fig．3 Mesh generation of foil computation region

1．1．2 氣動(dòng)力計(jì)算結(jié)果

計(jì)算圖2所示箔片外流場(chǎng)，設(shè)外流場(chǎng)為地面靜止大氣場(chǎng)，溫度為293K，壓力為101 325Pa，箔片飛行馬赫數(shù)Ma＝0．7。采用CFD專業(yè)求解器FLUENT進(jìn)行流場(chǎng)計(jì)算。選擇Spalart－Allmaras湍流模型，基于密度的隱式求解器，其中的經(jīng)驗(yàn)常數(shù)取FLUENT中的默認(rèn)值［15－17］。計(jì)算區(qū)域外圍設(shè)置為壓力遠(yuǎn)場(chǎng)邊界，箔片的外表面設(shè)置為壁面邊界。計(jì)算時(shí)通過(guò)設(shè)置來(lái)流方向控制箔片的迎角變化，來(lái)流方向平行于zOy平面，采用分離隱式求解器，各物理量的離散采用二階迎風(fēng)格式［18－19］，迭代收斂準(zhǔn)則為殘差小于1×10－6。按照上述算例中設(shè)置，計(jì)算箔片迎角在0°～90°之間的阻力系數(shù)與升力系數(shù)，為了提高效率，迎角每隔10°計(jì)算一次。在應(yīng)用氣動(dòng)力系數(shù)進(jìn)行計(jì)算時(shí)，采用線性插值的方法得到各個(gè)角度的氣動(dòng)力系數(shù)。具體計(jì)算結(jié)果如表1所示。

得到的阻力系數(shù)和升力系數(shù)隨迎角的變化規(guī)律如圖4所示，可以看出單個(gè)箔片的阻力系數(shù)在迎角0°～90°之間，隨迎角的增大而增大，升力系數(shù)隨迎角的增大先增大后減小，當(dāng)箔片的迎角為45°時(shí)，阻力系數(shù)和升力系數(shù)基本相等。由理論分析可得箔片在迎角為0°和90°時(shí)，升力系數(shù)應(yīng)該為0，表1中的計(jì)算值是由于誤差造成的，因而在下文計(jì)算時(shí)將其修正為0。

表1 箔片在不同迎角下的阻力系數(shù)和升力系數(shù)Table 1 Foil’s drag and lift coefficients at different angles of attack

圖4 箔片阻力系數(shù)和升力系數(shù)隨迎角的變化規(guī)律Fig．4 Change of foil’s drag and lift coefficients with angle of attack

由圖5可以看出，隨著迎角的增大，氣流流經(jīng)箔片上表面時(shí)，前緣區(qū)的負(fù)壓增大，后緣附近的附面層因受較強(qiáng)逆壓梯度作用，分離點(diǎn)前移，渦流區(qū)擴(kuò)大，使得上下表面壓力差和升力系數(shù)隨迎角增大的增長(zhǎng)率下降，阻力系數(shù)不斷增大；當(dāng)迎角達(dá)到30°左右（即臨界迎角）后，繼續(xù)增大迎角，分離點(diǎn)迅速前移，渦流區(qū)迅速擴(kuò)大，致使箔片上表面前段流速減慢，吸力降低。圖5為在不同迎角下，單個(gè)箔片計(jì)算域氣流對(duì)稱面的壓力場(chǎng)分布，圖中p為壓力。

圖5 不同迎角下單個(gè)箔片計(jì)算域的壓力場(chǎng)分布Fig．5 Pressure field distribution of single foil’s computation region at different angles of attack

在靠近后緣的一段范圍內(nèi)，吸力雖然稍有增高，但很有限，補(bǔ)償不了前段的吸力降低，所以，升力系數(shù)減小，阻力系數(shù)增大。

箔片在空中飛行的時(shí)候由于箔片的對(duì)稱性可得，箔片迎角的變化范圍為－90°～90°之間。當(dāng)迎角在－90°～0°之間變化時(shí)，阻力系數(shù)曲線與迎角在0°～90°之間的阻力系數(shù)曲線關(guān)于阻力系數(shù)軸對(duì)稱，即互為相反數(shù)的兩角阻力系數(shù)相等；而升力系數(shù)曲線與迎角在0°～90°之間的升力系數(shù)曲線關(guān)于原點(diǎn)中心對(duì)稱，即互為相反數(shù)的兩角的升力系數(shù)互為相反數(shù)。

1．2 兩個(gè)箔片氣動(dòng)力計(jì)算

面源誘餌剛從發(fā)射筒內(nèi)發(fā)射到有效燃燒前，箔片距離很近，并且在此過(guò)程中箔片還未發(fā)生翻轉(zhuǎn)，箔片之間距離比較近，可近似認(rèn)為彼此之間是相互平行的。因而相互之間氣動(dòng)力是有影響的，為了確定相對(duì)影響的大小，研究了兩個(gè)箔片正平行和交錯(cuò)平行兩種狀態(tài)下的氣動(dòng)力變化情況。

圖6為兩個(gè)正平行箔片計(jì)算域，計(jì)算域邊界到下箔片中心的距離與單片計(jì)算域設(shè)置相同。兩箔片平行于xOy面，定義箔片中心點(diǎn)z坐標(biāo)值較大的為上箔片，上下兩箔片之間的距離為2．54cm。

圖7為兩交錯(cuò)平行箔片計(jì)算域，計(jì)算域邊界到下箔片中心的距離與單片計(jì)算域設(shè)置相同。兩箔片平行于xOy面，定義箔片中心點(diǎn)z坐標(biāo)值較大的為上箔片，兩箔片直接垂直距離為1．27cm，將兩箔片沿y軸方向交錯(cuò)，距離為1．27cm。運(yùn)用CFD計(jì)算軟件對(duì)兩種平行狀態(tài)下箔片進(jìn)行計(jì)算時(shí)，參數(shù)設(shè)置與單個(gè)箔片計(jì)算時(shí)相同。

圖6 兩正平行箔片計(jì)算域Fig．6 Computation region of two orthographic parallel foils

圖7 兩交錯(cuò)平行箔片計(jì)算域Fig．7 Computation region of two offset parallel foils

表2給出了兩正平行箔片在不同迎角下的阻力系數(shù)和升力系數(shù)，圖8為兩箔片正平行時(shí)升力系數(shù)和阻力系數(shù)隨迎角的變化規(guī)律。表2和圖8中：CD1、CL1分別為上箔片的阻力系數(shù)、升力系數(shù)；CD2、CL2分別為下箔片的阻力系數(shù)、升力系數(shù)；CD、CL分別為單片計(jì)算時(shí)箔片的阻力系數(shù)、升力系數(shù)?？梢钥闯?，當(dāng)兩個(gè)箔片相距比較近時(shí)，彼此之間的流場(chǎng)互相影響，對(duì)氣動(dòng)力系數(shù)也產(chǎn)生了影響。位于來(lái)流方向前端的箔片（即下箔片）氣動(dòng)力系數(shù)與單個(gè)箔片的氣動(dòng)力系數(shù)基本相同，而位于來(lái)流方向后端的箔片（即上箔片）氣動(dòng)力系數(shù)在迎角較大時(shí)，氣動(dòng)力系數(shù)明顯小于單個(gè)箔片的氣動(dòng)力系數(shù)。

表2 兩正平行箔片在不同迎角下的阻力系數(shù)和升力系數(shù)Table 2 Drag and lift coefficients of two orthographic parallel foils at different angles of attack

表3給出了兩交錯(cuò)平行箔片在不同迎角下的阻力系數(shù)和升力系數(shù)，圖9為兩箔片交錯(cuò)平行時(shí)阻力系數(shù)和升力系數(shù)隨迎角的變化規(guī)律。表3和圖9中：CD3、CL3分別為下箔片的阻力系數(shù)、升力系數(shù)；CD4、CL4分別為上箔片的阻力系數(shù)、升力系數(shù)?？梢钥闯觯?dāng)兩個(gè)箔片交錯(cuò)平行時(shí)，位于來(lái)流方向前端的箔片（即下箔片）對(duì)位于來(lái)流方向后端的箔片（即上箔片）影響更大。隨著迎角的增大，影響越劇烈，特別是迎角在40°～60°之間時(shí)，由于下箔片的影響，上箔片的氣動(dòng)力系數(shù)為負(fù)值。

表3 兩交錯(cuò)平行箔片在不同迎角下的阻力系數(shù)和升力系數(shù)Table 3 Drag and lift coefficients of two offset parallel foils at different angles of attack

圖8 兩箔片正平行時(shí)升力和阻力系數(shù)隨迎角的變化Fig．8 Change of drag and lift coefficients of two orthographic parallel foils with angle of attack

圖9 兩箔片交錯(cuò)平行時(shí)升力和阻力系數(shù)隨迎角的變化Fig．9 Change of drag and lift coefficients of two offset parallel foils with angle of attack

圖10 不同迎角下兩正平行箔片計(jì)算域壓力場(chǎng)分布Fig．10 Pressure field distribution of computation region of two orthographic parallel foils at different angles of attack

圖10為在不同迎角下兩正平行箔片計(jì)算域氣流對(duì)稱面的壓力場(chǎng)分布。由圖10可以看出，當(dāng)兩箔片正平行時(shí)，在迎角比較小時(shí)，兩個(gè)箔片相互影響較小，氣動(dòng)力系數(shù)與單片計(jì)算時(shí)基本相同；隨著迎角的增大，由于上面箔片的影響，下面箔片后方的負(fù)壓區(qū)明顯大于單個(gè)箔片計(jì)算時(shí)的負(fù)壓區(qū)，前后壓力差增大，因此對(duì)于下面箔片而言其氣動(dòng)力系數(shù)要比單個(gè)箔片計(jì)算時(shí)的稍大一點(diǎn)；而對(duì)于上面箔片而言，來(lái)流方向前端負(fù)壓區(qū)的增大，使得其前后壓力差不斷減小，升力系數(shù)、阻力系數(shù)明顯減小，并且當(dāng)迎角大于80°時(shí)，前端的壓力小于后端的壓力，使得阻力系數(shù)為負(fù)。圖11為在不同迎角下兩交錯(cuò)平行箔片計(jì)算域氣流對(duì)稱面的壓力場(chǎng)分布。從圖11可以看出，由于兩箔片交錯(cuò)平行時(shí)，上箔片沿y軸方向與下箔片交錯(cuò)了1．27cm且兩箔片之間距離比正平行狀態(tài)下兩箔片之間的距離還小，因而上箔片在迎角為20°時(shí)，相較于正平行狀態(tài)下，前段負(fù)壓區(qū)明顯增大，因而氣動(dòng)力系數(shù)明顯小于上箔片氣動(dòng)力系數(shù)；特別是迎角在40°到60°之間變化時(shí)，上箔片完全置于下箔片的負(fù)壓區(qū)，在來(lái)流方向上，下箔片前段為負(fù)壓，后端為正壓，因而氣動(dòng)力系數(shù)為負(fù)值；隨著迎角的增大，上箔片前段正壓區(qū)開(kāi)始變大，因而氣動(dòng)力系數(shù)又變?yōu)檎怠?/p>

圖11 不同迎角下兩交錯(cuò)平行箔片計(jì)算域壓力場(chǎng)分布Fig．11 Pressure field distribution of computation region of two offset parallel foils at different angles of attack

2 箔片運(yùn)動(dòng)方程建立

由于實(shí)際使用中面源誘餌在高速運(yùn)動(dòng)的飛機(jī)上投放［20］，箔片初始速度Vf0為載機(jī)運(yùn)動(dòng)速度Vp與誘餌彈拋射速度Vf－p的矢量疊加，如圖12所示。

圖12 箔片初始運(yùn)動(dòng)的速度和方向Fig．12 Velocity and direction of foils initial motion

箔片在空中運(yùn)動(dòng)的整個(gè)過(guò)程中，始終受重力和空氣動(dòng)力的影響，但箔片表面同時(shí)還進(jìn)行著化學(xué)反應(yīng)，這使得箔片運(yùn)動(dòng)具有特殊性。在箔片剛被拋撒到空氣中時(shí)，各個(gè)箔片相互距離很近，與氧氣接觸也不充分，氧化反應(yīng)不劇烈。之后箔片與氧氣充分接觸，氧化反應(yīng)逐漸劇烈，釋放出大量的熱量，致使箔片與箔片之間的空氣受熱膨脹，箔片彼此之間產(chǎn)生一種熱力，稱之為燃燒力。

由于燃燒力作用在箔片表面并不是均勻的，因此燃燒力的作用點(diǎn)與箔片的幾何中心并不重合，對(duì)箔片產(chǎn)生一個(gè)翻轉(zhuǎn)力矩，從而使箔片質(zhì)心運(yùn)動(dòng)的同時(shí)姿態(tài)發(fā)生翻轉(zhuǎn)。可見(jiàn)燃燒作用使得箔片不同階段的運(yùn)動(dòng)具有不同受力特征。從發(fā)射到充分燃燒時(shí)間稱之為起燃時(shí)間ts，并根據(jù)起燃時(shí)間將箔片的運(yùn)動(dòng)分為起燃階段和充分燃燒兩個(gè)階段。

2．1 箔片起燃階段運(yùn)動(dòng)方程

起燃階段，各個(gè)箔片之間距離較近，不考慮相互碰撞，箔片所受的氣動(dòng)力會(huì)相互影響。根據(jù)氣動(dòng)力系數(shù)的計(jì)算結(jié)果，在單個(gè)箔片氣動(dòng)力系數(shù)的基礎(chǔ)上乘上一個(gè)修正系數(shù)即為相互干擾情況下箔片的平均氣動(dòng)力系數(shù)：

式中：aD和aL分別為阻力和升力系數(shù)的修正系數(shù)。根據(jù)牛頓第二定律，箔片在任意活動(dòng)坐標(biāo)系中的質(zhì)心動(dòng)力學(xué)方程為

式中：m為箔片的質(zhì)量；Vx、Vy和Vz分別為速度在x、y和z 3個(gè)坐標(biāo)軸上的分量；ωx、ωy和ωz為角速度在3個(gè)坐標(biāo)軸上的分量；Fx、Fy和Fz為合外力在3個(gè)坐標(biāo)軸上的分量。

圖13給出了本文所建立的地面坐標(biāo)系（Oxgygzg）、箔片坐標(biāo)系（O′xbybzb）和航跡坐標(biāo)系（O′xhyhzh）的關(guān)系。

圖13 地面坐標(biāo)系、箔片坐標(biāo)系和航跡坐標(biāo)系Fig．13 Ground coordinate system，foil coordinate system and flight－path coordinate system

為使箔片質(zhì)心動(dòng)力學(xué)方程具有最簡(jiǎn)單的形式，選用航跡坐標(biāo)系。將箔片所受的外力，即重力G和氣動(dòng)力R在航跡坐標(biāo)系中投影：

根據(jù)航跡坐標(biāo)系的定義，速度V在航跡坐標(biāo)系中只有沿Oxh方向的分量，且Vxh＝V，而Vyh＝Vzh＝0。

航跡坐標(biāo)系的轉(zhuǎn)動(dòng)角速度ωh可看做是沿Oyg軸方向的角速度 ψ和沿Ozh軸方向的角速度 θ的合矢量：

將上述諸力、速度和角速度在航跡軸系中的分量表達(dá)式代入式（5），并經(jīng)整理得到

這就是在航跡坐標(biāo)系中的箔片質(zhì)心動(dòng)力學(xué)方程。方程式中還隱含一個(gè)迎角α，它決定了升力和阻力的大小。

根據(jù)空間幾何的知識(shí)，可由箔片速度方向和箔片姿態(tài)計(jì)算出速度滾轉(zhuǎn)角γ和迎角α。箔片速度方向由θ和ψ確定，箔片姿態(tài)由θs和ψs確定，則由θ、ψ和θs、ψs可以求出α和γ：

箔片質(zhì)心的運(yùn)動(dòng)學(xué)方程為

給出箔片初始位置，初始速度大小V，初始速度方向θ0和ψ0，以及箔片初始姿態(tài)θs0和ψs0之后，積分方程式（9）和式（14）即可得到充分燃燒前箔片的飛行軌跡，即箔片空間位置（xg，yg，zg）隨時(shí)間的變化規(guī)律。

2．2 箔片充分燃燒階段運(yùn)動(dòng)方程

在箔片充分燃燒階段，各個(gè)箔片之間距離較遠(yuǎn)，箔片所受的氣動(dòng)力互不影響，因此單個(gè)箔片的氣動(dòng)力系數(shù)直接采用1．1．2節(jié)的數(shù)值計(jì)算結(jié)果。而箔片的動(dòng)力學(xué)方程和運(yùn)動(dòng)學(xué)方程仍是方程式（9）和式（14）的形式。

同時(shí)，由于氣動(dòng)力和燃燒力作用在箔片表面的非均勻性，箔片會(huì)發(fā)生翻轉(zhuǎn)，因此具有一個(gè)旋轉(zhuǎn)角速度。箔片的轉(zhuǎn)動(dòng)可分解為中心軸兩個(gè)姿態(tài)角的轉(zhuǎn)動(dòng)：

式中：ωθs為中心軸俯仰角的旋轉(zhuǎn)角速度；ωψs為中心軸偏轉(zhuǎn)角的旋轉(zhuǎn)角速度。為簡(jiǎn)化計(jì)算和分析，假設(shè)ωθs和ωψs不隨時(shí)間變化。

已知充分燃燒前階段的最后時(shí)刻（也就是ts時(shí)刻）箔片的位置、速度大小和方向以及箔片的姿態(tài)之后，積分方程式（9）、式（14）和式（15）即可得到充分燃燒后箔片的飛行軌跡。

設(shè)置上千個(gè)箔片初始姿態(tài)、速度和旋轉(zhuǎn)角速度的概率分布情況，同時(shí)對(duì)上千個(gè)箔片求解運(yùn)動(dòng)方程，即可得到整個(gè)面源干擾彈的空間運(yùn)動(dòng)擴(kuò)散規(guī)律。

3 面源誘餌運(yùn)動(dòng)擴(kuò)散仿真

在第3節(jié)單個(gè)箔片運(yùn)動(dòng)模型的基礎(chǔ)上，設(shè)置上千個(gè)箔片初始姿態(tài)和旋轉(zhuǎn)角速度的概率分布情況，同時(shí)對(duì)上千個(gè)箔片求解運(yùn)動(dòng)方程，即可得到整個(gè)面源干擾彈的空間運(yùn)動(dòng)擴(kuò)散規(guī)律。

一枚面源紅外誘餌的有效載荷通常包含一千余個(gè)箔片，雖然疊壓緊密，但根據(jù)文獻(xiàn)［21］中的內(nèi)彈道特性分析，各個(gè)箔片出口速度是有細(xì)微差別的，箔片束出口是一個(gè)加速過(guò)程，第一片的出口速度最小，最后一片的出口速度最大，本文經(jīng)過(guò)與實(shí)測(cè)數(shù)據(jù)的對(duì)比分析，假設(shè)一千個(gè)箔片的出口速度Vf0服從U（25，35）的均勻分布（單位 m／s）。同時(shí)由于大氣擾動(dòng)等隨機(jī)因素的影響，一千個(gè)箔片的初始姿態(tài)也有一定的差別，經(jīng)過(guò)與實(shí)測(cè)數(shù)據(jù)的對(duì)比分析，近似認(rèn)為箔片中心軸初始俯仰角θd0服從N（π／2，π2／300）的正態(tài)分布，初始偏轉(zhuǎn)角ψd0服從N（π／2，π2／500）的正態(tài)分布（單位rad）。

箔片充分燃燒后，由于氣動(dòng)力和燃燒力作用在箔片表面的非均勻性，箔片開(kāi)始旋轉(zhuǎn)，由于氣動(dòng)作用和燃燒情況具有一定的隨機(jī)性，一千余個(gè)箔片的旋轉(zhuǎn)角速度也不同，經(jīng)過(guò)與實(shí)測(cè)數(shù)據(jù)的對(duì)比分析，近似認(rèn)為箔片中心軸俯仰角的旋轉(zhuǎn)角速度ωθd和中心軸偏轉(zhuǎn)角的旋轉(zhuǎn)角速度ωψd均服從U（－16π，16π）的均勻分布（單位rad）。

由于高空高速發(fā)射平臺(tái)上，面源誘餌運(yùn)動(dòng)擴(kuò)散情況很難得到實(shí)測(cè)數(shù)據(jù)或圖像，因此為驗(yàn)證仿真模型的正確性，先計(jì)算一個(gè)在地面高速運(yùn)動(dòng)平臺(tái)上發(fā)射面源誘餌的擴(kuò)散算例。

設(shè)一枚面源紅外誘餌有1 000個(gè)箔片，單個(gè)箔片質(zhì)量為0．65g，直徑為25．4mm，起燃時(shí)間為0．08s。箔片垂直向上發(fā)射，發(fā)射平臺(tái)的運(yùn)動(dòng)馬赫數(shù)為0．7，發(fā)射點(diǎn)高度為3m。同時(shí)求解1 000個(gè)箔片的運(yùn)動(dòng)方程，得到正側(cè)方面源誘餌隨時(shí)間的運(yùn)動(dòng)擴(kuò)散情況如圖14所示，正后方面源誘餌隨時(shí)間的運(yùn)動(dòng)擴(kuò)散情況如圖15所示，1s時(shí)刻箔片三維分布散點(diǎn)圖如圖16所示。仿真圖中，y坐標(biāo)為負(fù)值時(shí)，表示箔片已經(jīng)落地。

圖14 xOy面內(nèi)面源誘餌隨時(shí)間的運(yùn)動(dòng)擴(kuò)散圖像Fig．14 Motion diffusion image of surface－type infrareddecoy in xOy

圖15 zOy面內(nèi)面源誘餌隨時(shí)間的運(yùn)動(dòng)擴(kuò)散圖像Fig．15 Motion diffusion image of surface－type infrared decoy in zOy

圖16 t＝1s箔片三維分布散點(diǎn)圖Fig．16 Three－dimensional scatter diagram of foils at t＝1s

由圖14～圖16仿真的面源誘餌擴(kuò)散圖像可以看出，在高速運(yùn)動(dòng)平臺(tái)上垂直向上發(fā)射面源紅外誘餌，其擴(kuò)散形狀近似為前部稀疏后部密集的錐形云團(tuán)。面源誘餌發(fā)射出去后迅速擴(kuò)散，在0．8s之內(nèi)就能擴(kuò)散成較大的形狀，形成面源假目標(biāo)。誘餌所受的氣動(dòng)阻力很大，因此沿初始運(yùn)動(dòng)方向向前運(yùn)動(dòng)較短的距離（大多數(shù)箔片均不超過(guò)100 m）就開(kāi)始下落，約2s之后就開(kāi)始垂直下落。可見(jiàn)在高速運(yùn)動(dòng)的平臺(tái)上發(fā)射面源誘餌，形成的擴(kuò)散假目標(biāo)基本上滯留在發(fā)射點(diǎn)附近，而發(fā)射平臺(tái)會(huì)迅速遠(yuǎn)離。因此只有在探測(cè)方向與平臺(tái)運(yùn)動(dòng)方向的夾角較小或探測(cè)器距離較遠(yuǎn)時(shí)，面源假目標(biāo)才有可能落入探測(cè)器視場(chǎng)較長(zhǎng)時(shí)間，形成有效干擾。

表4給出了面源誘餌運(yùn)動(dòng)擴(kuò)散仿真值與實(shí)測(cè)值的對(duì)比，其中中心點(diǎn)位置是相對(duì)于發(fā)射點(diǎn)的位置。由于1s之后大量箔片已經(jīng)落地，在y方向上仿真值與實(shí)測(cè)值出現(xiàn)較大偏差。但總體上看，仿真得到的面源誘餌擴(kuò)散尺寸和運(yùn)動(dòng)位置也是較為符合實(shí)際。特別是x方向擴(kuò)散尺寸誤差小于5%，說(shuō)明上述所建立的模型以及對(duì)箔片氣動(dòng)力系數(shù)的計(jì)算是正確的。誤差的存在主要是在實(shí)際測(cè)量中還有一些不可控因素的存在，例如風(fēng)的因素，儀器的測(cè)量誤差等。

表4 面源誘餌運(yùn)動(dòng)擴(kuò)散仿真值與實(shí)測(cè)值的對(duì)比Table 4 Comparison of surface－type infrared decoy motion diffusion simulation and measured values

4 結(jié) 論

1）通過(guò)計(jì)算兩個(gè)箔片正平行和交錯(cuò)平行兩種狀態(tài)下的氣動(dòng)力系數(shù)。發(fā)現(xiàn)來(lái)流后端的箔片處于來(lái)流前端箔片的負(fù)壓區(qū)，因而處于來(lái)流后端的箔片氣動(dòng)力系數(shù)變化很大。

2）在建立箔片運(yùn)動(dòng)方程時(shí)，考慮到燃燒力的存在，因而將箔片的運(yùn)動(dòng)分為充分燃燒前和充分燃燒兩個(gè)階段，從而使仿真結(jié)果與試驗(yàn)結(jié)果更吻合。

3）由仿真的面源誘餌擴(kuò)散圖像可以看出，在高速運(yùn)動(dòng)平臺(tái)上垂直向上發(fā)射面源紅外誘餌，其擴(kuò)散形狀近似為前部稀疏后部密集的長(zhǎng)扁形云團(tuán)。將仿真數(shù)據(jù)與實(shí)測(cè)數(shù)據(jù)對(duì)比可知，仿真結(jié)果與試驗(yàn)結(jié)果非常接近，說(shuō)明所建模型的正確性。

4）本文所建立的面源誘餌擴(kuò)散模型對(duì)機(jī)載面源誘餌擴(kuò)散規(guī)律進(jìn)行了更深入的研究，通用性好；可服務(wù)于機(jī)載面源紅外誘餌的設(shè)計(jì)，為誘餌技術(shù)指標(biāo)的制定提供理論依據(jù)；還可服務(wù)于紅外成像探測(cè)器的改進(jìn)，為成像探測(cè)器的抗干擾算法設(shè)計(jì)提供理論參考和干擾模型。

［1］ GOPALAPILLAI R，KIM H D，SETOGUCHI T，et al．On the near－field aerodynamics of a projectile launched from a ballistic range［J］．Journal of Mechanical Science and Technology，2007，21（7）：1129－1138．

［2］ FEGG M，BANNASH H．Method for offering aphantom target and decoy：America，US6513438B1［P］．2003－02－03．

［3］ SCHLEIJPEN R H M A．Toolset for evaluating infrared countermeasures and signature reduction for ships［C］／Proceedings of SPIE－The Intenational Society for Optical Engineering．Toulouse，F(xiàn)rance：SPIE，2010，7836：78360F－1－79360F－12．

［4］ KOCH E C．Review on pyrotechnic aerial infrared decoys［J］． Propellants， Explosives， Pyrotechnics， 2001，26（1）：3－11．

［5］ BUTTERS B，NICHOLLS E，WALMSLEY R，et al．Infrared decoy and obscurant modeling and simulation for ship protection［C］／Proceedings of SPIE－The Intenational Society for Optical Engineering．Toulouse，F(xiàn)rance：SPIE，2011，8187（1）：75－78．

［6］ 趙非玉，盧山，蔣沖，等．面源紅外誘餌仿真建模方法研究［J］．光電技術(shù)應(yīng)用，2012，27（2）：66－69．ZHAO F Y，LU S，JIANG C，et al．Research on the modeling method of surface－type infrared decoy［J］．Electro－Optic Technology Application，2012，27（2）：66－69（in Chinese）．

［7］ 付曉紅，樊秋林．面源紅外誘餌的建模與仿真研究［J］．光電技術(shù)應(yīng)用．2013，28（6）：81－86．FU X H，F(xiàn)AN Q L．Research on modeling and simulation of special material decoy ［J］．Electro－Optic Technology Application，2013，28（6）：81－86（in Chinese）．

［8］ 陳乃光．光電干擾箔片云團(tuán)運(yùn)動(dòng)規(guī)律的仿真研究［J］．航天電子對(duì)抗，2008，24（5）：20－23．CHEN N G．Moving characters simulation of foils cloud used in optics－electronic countermeasure［J］．Aerospace Electronic Warfare，2008，24（5）：20－23（in Chinese）．

［9］ 黃蓓，王浩，王帥，等．光電干擾箔片拋撤與擴(kuò)散試驗(yàn)研究［J］．含能材料，2012（5）：605－609．HUANG B，WANG H，WANG S，et al．Dispersing experiment for separation and dispersion of multi－plates［J］．Chinese Journal of Energetic Materials，2012（5）：605－609（in Chinese）．

［10］ BADURA W．Device for producing a decoy cloud，in particular an infrared decoy cloud：America，US4621579A［P］．1986－11－11．

［11］ 郭江濤，邢存震，譚獻(xiàn)忠，等．基于Fluent的亞音速翼型氣動(dòng)特性數(shù)值研究［J］．中原工學(xué)院學(xué)報(bào)，2012，23（1）：36－40．GUO J T，XING C Z，TAN X Z，et al．Numerical study on subsonic airfoil aerodyamic characteristics based on the Fluent［J］．Journal of Zhongyuan University of Technology，2012，23（1）：36－40（in Chinese）．

［12］ HU K，DING Y，ZHANG X Y，et al．A scalable infrastructure for online performance analysis on CFD application［J］．Chinese Journal of Aeronautics，2012，25（4）：546－558．

［13］ OLLIE H．Technology survey：Airborne dispensers and ir expendables［J］．Journal of Electronic Defense，2011，34（6）：59－66．

［14］ KOCH E C．2006－2008Annual review on aerial infrared decoy flares［J］．Propellants Explos Pyrotech，2009，34（1）：6－12．

［15］ FENG X Q，LI Z K，SONG B F．Research of low boom and low drag supersonic aircraft design［J］．Chinese Journal of Aeronautics，2014，27（3）：531－541．

［16］ 李鵬，招啟軍．懸停狀態(tài)傾轉(zhuǎn)旋翼／機(jī)翼干擾流場(chǎng)及氣動(dòng)力的CFD計(jì)算［J］．航空學(xué)報(bào)，2014，35（2）：361－371．LI P，ZHAO Q J．CFD calculations on the interaction flow field and aerodynamic force of tiltrotor／wing in hover［J］．Acta Aeronautica et Astronautica Sinica，2014，35（2）：361－371（in Chinese）．

［17］ 劉雄，王翼，梁劍寒．方轉(zhuǎn)圓對(duì)三維側(cè)壓進(jìn)氣道的流動(dòng)特性影響［J］．航空學(xué)報(bào)，2014，35（11）：2939－2948．LIU X，WANG Y，LIANG J H．Effect of rectangular－tocircular on flow characteristics of three－dimensional sidewall compression inlet［J］．Acta Aeronautica et Astronautica Sinica，2014，35（11）：2939－2948（in Chinese）．

［18］ 郭則慶，王楊，姜孝海，等．膛口初始流場(chǎng)對(duì)火藥燃?xì)饬鲌?chǎng)影響的數(shù)值研究［J］．兵工學(xué)報(bào)，2012，33（6）：663－668．GUO Z Q，WANG Y，JIANG X H，et al．Numerical study on effects of precursor flow on muzzle proprellant flow field［J］．Acta Armamen Taii，2012，33（6）：663－668（in Chinese）．

［19］ WU Z J，XU Y Z，WANG W B，et al．Re view of shock wave detection method in CFD post－processing［J］．Chinese Journal of Aeronautics，2013，26（3）：501－513．

［20］ 馬東立，劉亞楓，林鵬．航空拖曳誘餌系統(tǒng)的動(dòng)態(tài)特性研究［J］．航空學(xué)報(bào)，2014，35（1）：161－170．MA D L，LIU Y F，LIN P．Study of dynamic characteristics of aeronautic toeed decoy system［J］．Acta Aeronautica et Astronautica Sinica，2014，35（1）：161－170（in Chinese）．

［21］ 黃長(zhǎng)強(qiáng)，趙輝，杜海文，等．機(jī)載彈藥精確制導(dǎo)原理［M］．北京：國(guó)防工業(yè)出版社，2011：164－169．HUANG C Q，ZHAO H，DU H W，et al．The precision guidance theory of airborne mention［M］．Beijing：National Defense Industry Press，2011：164－169（in Chinese）．

Diffusion rule of foil－surface－type infrared decoy

ZOU Tao，WANG Chaozhe，TONG Zhongxiang＊，JIA Lintong，TONG Qi
Aeronautics and Astronautics Engineering Institute，Air Force Engineering University，Xi’an 710038，China

The surface－type infrared decoy has become an important method against infrared imaging guided missile．In this paper，foil－surface－type infrared decoy is regarded as the research object and the computational fluid dynamics（CFD）is used to obtain the aerodynamic coefficients of a single foil and two parallel foils at different angles of attack．With the analysis of foil movement in the air，the movement of the foil is divided into two phases，i．e．，the light－off and complete combustion considering the particularity of combustion pressure．And the motion models of the two stages are established．Based on the motion model of the single foil，thousands of foils’initial probability distribution of posture and rotating angular velocity are set and then the motion equations of thousands of foils are solved at the same time．In this way，the whole space motion diffusion rule of foil－surface－type infrared decoy is obtained．The simulation results show that the diffusion shape of surfacetype infrared decoy，which is launched vertically from high－speed motion platform，is the taper cloud with sparse front and dense rear approximately．The diffusion size and motion position of simulation are highly identical with the measured value，and the error is less than 5%of the diffusion size in x direction in particular．

foil；surface－type infrared decoy；computational fluid dynamics；aerodynamic coefficients；combustion force；diffusion rule

2015－10－26；Revised：2015－11－29；Accepted：2016－01－13；Published online：2016－03－02 14：49

URL：www．cnki．net／kcms／detail／11．1929．V．20160302．1449．006．html

National Natural Science Foundation of China（61471390）

V211．3；TN972

A

1000－6893（2016）09－2634－12

10．7527／S1000－6893．2016．0023

2015－10－26；退修日期：2015－11－29；錄用日期：2016－01－13；網(wǎng)絡(luò)出版時(shí)間：2016－03－02 14：49

www．cnki．net／kcms／detail／11．1929．V．20160302．1449．006．html

國(guó)家自然科學(xué)基金 （61471390）

＊通訊作者．Tel．：029－84787505 E－mail：tzxkgy＠163．com

鄒濤，王超哲，童中翔，等．箔片型紅外面源誘餌擴(kuò)散規(guī)律［J］．航空學(xué)報(bào)，2016，37（9）：26342－645．ZOU T，WANG C Z，TONG Z X，et al．Diffusion rule of foils－urfacet－ype infrared decoy［J］．Acta Aeronautica et Astronautica Sinica，2016，37（9）：26342－645．

鄒濤 男，碩士研究生。主要研究方向：飛行器紅外對(duì)抗仿真。E－mail：zoutao＿af＠163．com

童中翔 男，博士，教授，博士生導(dǎo)師。主要研究方向：飛行性能和飛行仿真，飛行器紅外對(duì)抗技術(shù)。

Tel．：029－84787505

E－mail：tzxkgy＠163．com

＊Corresponding author．Tel．：029－84787505 E－mail：tzxkgy＠163．com