陳 昊,孫雪明,劉建功,阮文俊

(1.南京理工大學(xué) 能源與動(dòng)力工程學(xué)院,江蘇 南京 210094;2.江西星火軍工工業(yè)有限公司,江西 南昌 331709;3.河北太行機(jī)械工業(yè)有限公司,河北 石家莊 052160)

箔條能在一定的空間范圍產(chǎn)生干擾回波[1],使得雷達(dá)無(wú)法準(zhǔn)確識(shí)別出目標(biāo)而達(dá)不到預(yù)期攻擊效果[2]。研究箔條云從產(chǎn)生到擴(kuò)散過(guò)程中的形態(tài)大小和取向分布等性質(zhì),對(duì)探究其電磁散射特性的變化極為重要[3]。

ALVAREZ等[4]提出了一種用于一般箔條云電磁模擬的離散箔條云模型,可以用來(lái)計(jì)算諧波或?qū)拵}沖散射響應(yīng)。該技術(shù)考慮了箔條云的幾何形狀和箔條在密度和方向上的統(tǒng)計(jì)分布。SEO等[5]運(yùn)用六自由度擴(kuò)散模型得到了低雷諾數(shù)下箔條位置和取向隨時(shí)間的變化情況。楊學(xué)斌等[6]在忽略箔條質(zhì)量的前提下建立了箔條云團(tuán)的布朗運(yùn)動(dòng)擴(kuò)散模型。李志輝等[7,8]應(yīng)用隨機(jī)動(dòng)力學(xué)加權(quán)技術(shù)及概率統(tǒng)計(jì)分布原理對(duì)箔條云整體性能進(jìn)行了分析。曾濤等[9]基于對(duì)箔條云邊界處若干高速箔條位置的計(jì)算提出了一種箔條云整體運(yùn)動(dòng)模型。TONG等[10]為了實(shí)現(xiàn)箔條快速散開(kāi)提出將一種高速旋轉(zhuǎn)的飛行器作為箔條投放器。

現(xiàn)有的研究為模擬箔條云飛行覆蓋范圍及整體運(yùn)動(dòng)性能提供了方法,但存在以下不足:①部分學(xué)者提出的模型只適用于成熟箔條云的研究[6],無(wú)法對(duì)處于發(fā)展過(guò)程中的箔條云及其性能進(jìn)行分析;②箔條云模型大多由球形發(fā)展至橢球型[7-11],與實(shí)戰(zhàn)過(guò)程中箔條云形態(tài)發(fā)展的有一定差距;③隨著近年來(lái)反箔條干擾工作的進(jìn)展[12,13],僅研究云團(tuán)的宏觀(guān)外形發(fā)展已不再能滿(mǎn)足有效干擾的需求,云團(tuán)內(nèi)部箔條如何分布已成為高成功率干擾的關(guān)鍵。

鑒于以上問(wèn)題,本文運(yùn)用剛體動(dòng)力學(xué)理論和碰撞概率模型,在單根箔條運(yùn)動(dòng)模型的基礎(chǔ)上,研究具有高牽連速度的箔條集群動(dòng)態(tài)拋撒后的多體分離過(guò)程,得到箔條云團(tuán)軸向截面由錐形發(fā)展至矩形的分布模型,并與試驗(yàn)結(jié)果對(duì)比,驗(yàn)證計(jì)算模型的準(zhǔn)確性和有效性,在此基礎(chǔ)上分析處于擴(kuò)散階段及沉降階段云團(tuán)內(nèi)部箔條的分布規(guī)律。為后續(xù)研究不同牽連速度下箔條的動(dòng)態(tài)拋撒提供參考。

1 計(jì)算模型

1.1 剛體動(dòng)力學(xué)及運(yùn)動(dòng)模型

由于在模擬箔條拋撒的計(jì)算空間范圍內(nèi)箔條所占體積比例遠(yuǎn)低于氣相場(chǎng),箔條對(duì)氣相場(chǎng)的作用可以忽略不計(jì),因而本文只考慮氣相場(chǎng)對(duì)固相的作用,而忽略固相場(chǎng)對(duì)氣相場(chǎng)的反作用。

為避免過(guò)大的計(jì)算量,本文以一定量的取樣箔條群代替數(shù)目龐大的真實(shí)箔條群,通過(guò)跟蹤計(jì)算每一個(gè)箔條求解離散場(chǎng)。計(jì)算時(shí)將每一個(gè)箔條離散成若干等份的箔條小段,對(duì)箔條小段進(jìn)行受力分析,將箔條小段受力合成求得箔條所受合力,并求得箔條所受力矩;根據(jù)箔條所受合力及力矩,求其在計(jì)算空間范圍內(nèi)的運(yùn)動(dòng)。

如圖1所示,Oxyz為固定參考坐標(biāo)系的平移坐標(biāo)系,Oξηζ為體坐標(biāo)系。描述箔條姿態(tài)的3個(gè)角為進(jìn)動(dòng)角ψ、章動(dòng)角θ、自轉(zhuǎn)角φ。上述3個(gè)角即稱(chēng)為歐拉角,在一定的條件下,剛體的任一方位可用一組歐拉角唯一地表示。

圖1 歐拉角示意圖Fig.1 Diagram of Euler angles

此模型中沒(méi)有考慮其他性質(zhì)力對(duì)箔條運(yùn)動(dòng)的影響,而只考慮流場(chǎng)對(duì)箔條的曳力和箔條所受的重力,相關(guān)計(jì)算公式為

ua=u0+(I·ω)×ra

(1)

式中:ua為箔條離散分段a的速度;u0為箔條的質(zhì)心速度;I為體坐標(biāo)系到定坐標(biāo)系的轉(zhuǎn)換矩陣;ω為箔條在體坐標(biāo)系中的角速度;ra為定坐標(biāo)系中箔條離散分段相對(duì)于基準(zhǔn)點(diǎn)的位置矢量。

箔條受到的流體曳力計(jì)算公式[14]為

(2)

式中:FD為箔條所受的流體曳力;FD,a為箔條離散分段所受的流體曳力;n為箔條離散段數(shù);CD為流體曳力系數(shù);ρc為流場(chǎng)密度;Ap為曳力作用的有效截面積;uc為流場(chǎng)速度;ua為離散分段的速度。

Ap=πdlasinβ

(3)

式中:d為箔條直徑,la為箔條離散分段的長(zhǎng)度,β為箔條軸向與矢量uc-ua的夾角。

曳力系數(shù)CD的計(jì)算公式[15]為

CD=0.71+0.07RN

(4)

(5)

式中:Re為箔條的雷諾數(shù);l為箔條的特征長(zhǎng)度,對(duì)于圓柱形箔條而言為箔條的直徑;μ為空氣的動(dòng)力黏度。對(duì)于單根箔條,平動(dòng)方程為

(6)

式中:m為細(xì)長(zhǎng)箔條的質(zhì)量;g為重力加速度;s為箔條質(zhì)心的位移。轉(zhuǎn)動(dòng)方程為[16]:

(7)

(8)

在運(yùn)動(dòng)學(xué)計(jì)算過(guò)程中引入歐拉參數(shù)。歐拉參數(shù)實(shí)際上是一組四元素,由λ0,λ1,λ2和λ34個(gè)元素組成。并且這4個(gè)元素滿(mǎn)足以下關(guān)系[17-19]:

(9)

歐拉參數(shù)與歐拉角關(guān)系[20]為

(10)

歐拉參數(shù)表示的定點(diǎn)轉(zhuǎn)動(dòng)的運(yùn)動(dòng)方程為

(11)

再根據(jù)歐拉角和歐拉參數(shù)之間的關(guān)系求出一個(gè)時(shí)間步之后的歐拉角。

(12)

式中:ψ′,θ′,φ′為一個(gè)計(jì)算時(shí)間步后的進(jìn)動(dòng)角、章動(dòng)角和自轉(zhuǎn)角。

1.2 箔條間碰撞模型

為考慮箔條間碰撞,對(duì)計(jì)算空間進(jìn)行網(wǎng)格劃分,箔條間的碰撞在這些網(wǎng)格內(nèi)進(jìn)行。為了找到不同箔條發(fā)生碰撞的位置,本文在箔條離散分段的基礎(chǔ)上,將這些離散分段轉(zhuǎn)化為具有等效直徑的小球,運(yùn)用修正后的Nanbu模型,得箔條分段a和同一網(wǎng)格內(nèi)其他所有箔條分段(除與分段a處在同一根箔條上的其他分段外)的碰撞概率[21]為

(13)

式中:n和N分別為箔條分段的真實(shí)數(shù)目和取樣數(shù)目;Dp為箔條分段的等效直徑;Gab為箔條分段a和b的相對(duì)速度;go為徑向分布函數(shù);Δt為時(shí)間步長(zhǎng)。

go=1/[1-εs/εmax]1/3

(14)

式中:εmax為填充狀態(tài)時(shí)的箔條分段密度,εs為不同時(shí)刻各網(wǎng)格內(nèi)的箔條分段密度。

在箔條分段a總碰撞概率Pa<1時(shí),利用隨機(jī)數(shù)R(0(b/N-Pab)則發(fā)生碰撞。

假設(shè)箔條分段a、b分別在箔條i、j上,通過(guò)上述方法就找到了箔條i、j發(fā)生碰撞的位置,根據(jù)動(dòng)量守恒定律可以獲得箔條i和j發(fā)生碰撞后的速度和角速度。具體的碰撞模型[22]為

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

式中:下標(biāo)P代表平行,V代表垂直;ui,uj為箔條碰撞前碰點(diǎn)速度;vi0,vj0為箔條碰撞前質(zhì)心速度;Ii,Ij為坐標(biāo)轉(zhuǎn)換余弦矩陣;ωi0,ωj0為碰撞前箔條在體坐標(biāo)系中的角速度;ri,rj為碰點(diǎn)到箔條質(zhì)心的位置矢量;v為碰點(diǎn)最大壓縮形變時(shí)速度;ri0,rj0分別為ri,rj的單位矢量;Pi,Pj為ri0,rj0的并矢張量;viV,vjV為在ri0,rj0垂直方向的速度變化;qi,qj為沖量;e為恢復(fù)系數(shù);vi,vj為碰撞后箔條質(zhì)心的速度;ωi,ωj為碰撞后箔條在體坐標(biāo)系中角速度;Π為定坐標(biāo)系中的單位張量;J′i,J″j為箔條在體坐標(biāo)系中的轉(zhuǎn)動(dòng)慣量;Ji,Jj為箔條在定坐標(biāo)系中的轉(zhuǎn)動(dòng)慣量;e為碰撞過(guò)程中的恢復(fù)系數(shù),當(dāng)彈性碰撞時(shí)e=1,當(dāng)完全非彈性碰撞時(shí)e=0。

1.3 計(jì)算過(guò)程

由式(1)～式(22)即可得涉及相互間碰撞的箔條動(dòng)力學(xué)及運(yùn)動(dòng)學(xué)特征,其中微分方程組由四階龍格-庫(kù)塔法求解。具體過(guò)程由C++實(shí)現(xiàn)。運(yùn)動(dòng)參數(shù)計(jì)算流程圖如圖2所示。

圖2 運(yùn)動(dòng)參數(shù)計(jì)算流程圖Fig.2 Flow chart of motion parameter calculation

2 數(shù)值計(jì)算和試驗(yàn)結(jié)果對(duì)比

設(shè)定載體火箭的飛行方向?yàn)閥軸正向,飛行馬赫數(shù)為0.8;箔條集束以35 m/s的速度垂直于載體火箭向上出艙,開(kāi)艙點(diǎn)的坐標(biāo)為(2,1,1),設(shè)定開(kāi)艙時(shí)刻為零時(shí)刻,開(kāi)艙后箔條由集束表面分離。假設(shè)不同時(shí)刻由箔條集束表面上分離的箔條的進(jìn)動(dòng)角ψ、章動(dòng)角θ、自轉(zhuǎn)角φ分別在0～2π,0～π,0～2π上呈隨機(jī)分布,角速度ω在體坐標(biāo)系3個(gè)軸上的投影ωξξ,ωηη,ωζζ在0～2 rad/s上隨機(jī)分布。仿真計(jì)算中箔條總數(shù)為2 000 根,時(shí)間步長(zhǎng)為1 ms,計(jì)算總時(shí)長(zhǎng)為1.5 s。

為了驗(yàn)證模型的合理性,進(jìn)行了近地箔條拋撒試驗(yàn)。針對(duì)本文試驗(yàn)背景需求,設(shè)計(jì)了固體火箭發(fā)動(dòng)機(jī)攜帶干擾彈高速拋撒的箭載試驗(yàn)系統(tǒng)。點(diǎn)火后火箭加速至干擾彈發(fā)射所需的模擬速度時(shí),采用延時(shí)點(diǎn)火炬點(diǎn)燃干擾彈發(fā)射藥,箔條集束被拋出。試驗(yàn)時(shí)高速攝像機(jī)布置在彈道側(cè)方,用于記錄箔條集束分離、云團(tuán)擴(kuò)散的全過(guò)程。試驗(yàn)中通過(guò)六通道同步觸發(fā)器控制火箭與高速攝像機(jī)的同步點(diǎn)火與觸發(fā)。圖3為箔條拋撒后不同時(shí)刻yz截面云團(tuán)分布的試驗(yàn)結(jié)果和仿真結(jié)果。圖4為箔條拋撒后xz截面云團(tuán)仿真結(jié)果。

圖3 yz截面云團(tuán)分布圖Fig.3 Chaff cloud distribution in the yz section

圖4 xz截面云團(tuán)分布圖Fig.4 Chaff cloud distribution in the xz section

由圖3和圖4中不同時(shí)刻云團(tuán)散布對(duì)比可知,當(dāng)所有箔條拋撒完畢,即t=50 ms時(shí),箔條云團(tuán)呈較為整齊的錐形分布,y方向分布在4～13 m,z方向分布在1.5～3 m。t=100 ms時(shí),箔條云已基本成型,x方向上分布在1～3 m,yz方向上的最大散布面積已達(dá)10 m2以上?？梢?jiàn),箔條多體分離用時(shí)很短,能夠在100 ms內(nèi)迅速散開(kāi)成型。300 ms后云團(tuán)分布基本穩(wěn)定,進(jìn)入沉降階段,箔條云團(tuán)整體以一定速度下沉。

表1給出了箔條云團(tuán)軸向及徑向尺寸隨時(shí)間的變化情況,也反映了箔條云的遮蔽面積與時(shí)間的關(guān)系,試驗(yàn)結(jié)果與仿真結(jié)果發(fā)展趨勢(shì)一致。當(dāng)云團(tuán)形態(tài)達(dá)到穩(wěn)態(tài)時(shí),其在軸向上的散布面積可達(dá)13 m2,徑向上的散布面積可達(dá)4 m2。同仿真結(jié)果相同的是,試驗(yàn)中云團(tuán)軸向截面由錐形截面演化為矩形截面。

表1 云團(tuán)尺寸變化試驗(yàn)與仿真結(jié)果對(duì)比Table 1 Comparison of cloud size change in experiment and simulation results

3 云團(tuán)散布特性參數(shù)統(tǒng)計(jì)分析

3.1 數(shù)量分布

數(shù)量分布η是箔條散布過(guò)程中重要的特征量。圖5給出了箔條數(shù)量分別沿x,y,z方向的分布情況,從而可以得出箔條云團(tuán)隨時(shí)間變化的散布規(guī)律。

圖5 不同方向箔條分布個(gè)數(shù)Fig.5 Number of chaff in different directions

由圖5可知,50 ms時(shí)在x方向上箔條集中在1.5～2.5 m之間,隨著云團(tuán)向四周的擴(kuò)散,箔條散布面積變大,分布逐漸均勻;在y方向上箔條的數(shù)量分布在100 ms后基本不再變化,軸向散布的結(jié)構(gòu)已經(jīng)形成,各時(shí)刻數(shù)量分布曲線(xiàn)的整體趨勢(shì)大致相同。z軸上的分布與x軸類(lèi)似,但由于重力作用,云團(tuán)進(jìn)入沉降階段后以一定速度沿重力方向運(yùn)動(dòng),表現(xiàn)為數(shù)量分布曲線(xiàn)整體以相對(duì)固定的趨勢(shì)隨時(shí)間沿z軸負(fù)向移動(dòng)。

3.2 擴(kuò)散分布的數(shù)字特征

通過(guò)討論擴(kuò)散過(guò)程中箔條位置坐標(biāo)D的均值、軸向平面與徑向平面內(nèi)箔條云團(tuán)邊界值和方差3個(gè)數(shù)字特征量,對(duì)云團(tuán)的空間分布進(jìn)行了分析。

3.2.1 箔條群坐標(biāo)均值

箔條集群坐標(biāo)的統(tǒng)計(jì)平均值反映了箔條云整體在不同時(shí)刻的所在位置,其計(jì)算方法為

(23)

式中:D為箔條質(zhì)心的位置坐標(biāo),N為箔條個(gè)數(shù)。通過(guò)計(jì)算可得幾個(gè)特征時(shí)刻箔條群空間坐標(biāo)統(tǒng)計(jì)平均值,如表2所示。

表2 箔條位置坐標(biāo)均值Table 2 Mean value of chaff position

以各方向上箔條位置坐標(biāo)的平均值作為箔條云質(zhì)心的位置坐標(biāo),即可得出云團(tuán)整體速度隨時(shí)間的變化曲線(xiàn),如圖6所示。300 ms后,箔條質(zhì)心速度在各方向上的分量基本不再發(fā)生改變,說(shuō)明云團(tuán)整體運(yùn)動(dòng)進(jìn)入沉降階段。

圖6 箔條云整體速度變化圖Fig.6 Change of overall velocity of chaff cloud

3.2.2 云團(tuán)邊界處箔條坐標(biāo)值

云團(tuán)邊界值是指某一時(shí)刻坐標(biāo)的最小值和最大值,可以對(duì)云團(tuán)外輪廓進(jìn)行描述。由圖7可以看出,y方向上的邊界值曲線(xiàn)在50～300 ms間斜率大幅度降低,這是由于拋撒完成后箔條速度很快衰減進(jìn)而以低速運(yùn)動(dòng);x方向上300 ms之前的云團(tuán)尺寸增長(zhǎng)速率大于300 ms之后的;300 ms后,所有箔條保持勻速運(yùn)動(dòng),使得整個(gè)云團(tuán)也以一定速度沉降,z方向上云團(tuán)上下邊界值同步減小。

圖7 不同方向上云團(tuán)邊界值變化圖Fig.7 Cloud boundary changes in different directions

3.2.3 箔條位置坐標(biāo)偏離均值的方差

將不同時(shí)刻各箔條位置和由式(23)計(jì)算得出的E代入式(24)中,得到箔條偏離平均值的方差。

(24)

圖8為不同時(shí)刻箔條云團(tuán)在各方向上的方差變化曲線(xiàn)。y方向上的方差明顯大于x,z方向,這是由于箔條云在y方向上散布區(qū)域遠(yuǎn)大于其他2個(gè)方向。300 ms前,方差在y方向上的增長(zhǎng)速度較快,這是因?yàn)榇藭r(shí)位于云團(tuán)頭部和尾部箔條的速度相差很大,云團(tuán)內(nèi)部箔條位置分布變化較大。300 ms后,方差在x,y,z方向上都基本不再變化,箔條分布結(jié)構(gòu)已基本形成。

圖8 3個(gè)方向方差變化圖Fig.8 Change in variance in x,y and z directions

4 結(jié)束語(yǔ)

本文在剛體動(dòng)力學(xué)及碰撞概率模型的基礎(chǔ)上,對(duì)箔條在牽連速度為270 m/s下的拋撒散布過(guò)程進(jìn)行了數(shù)值模擬,分析了云團(tuán)散布軌跡,并開(kāi)展試驗(yàn),驗(yàn)證了模型的準(zhǔn)確性,進(jìn)而對(duì)箔條云團(tuán)的整體性能及云團(tuán)內(nèi)部箔條的分布規(guī)律進(jìn)行了分析?？梢缘玫揭韵陆Y(jié)論:

①牽連速度是影響箔條云團(tuán)成型和散布的關(guān)鍵因素。箔條云團(tuán)在0.1 s內(nèi)基本成型,0.3 s后云團(tuán)內(nèi)箔條的數(shù)量分布情況基本不再隨時(shí)間變化。箔條集束能夠在0.3 s內(nèi)形成穩(wěn)態(tài)云團(tuán),云團(tuán)成熟后在軸向上的散布面積可達(dá)13 m2,徑向上的散布面積可達(dá)4 m2,云團(tuán)邊界值的發(fā)展趨勢(shì)在0.3 s時(shí)發(fā)生明顯變化。

②進(jìn)入沉降過(guò)程后,云團(tuán)整體運(yùn)動(dòng)速度穩(wěn)定,其中在y方向上的速度為0.12 m/s左右,z方向上的速度為-0.4 m/s左右。由于y方向上不同位置的箔條速度相差較大,因而在該方向上箔條分布更具不均勻性。

③在牽連速度為270 m/s的情況下拋撒箔條,仿真結(jié)果與實(shí)驗(yàn)結(jié)果都表明箔條云團(tuán)軸向截面經(jīng)歷了由錐形向長(zhǎng)方形變化的過(guò)程。說(shuō)明該模型能很好地描述高牽連速度下拋撒的云團(tuán)散布過(guò)程,為后續(xù)開(kāi)展不同牽連速度的箔條拋撒散布研究提供理論依據(jù)。