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

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

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

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

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

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

1 計算模型

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

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

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

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

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

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

ua=u0+(I·ω)×ra

(1)

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

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

(2)

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

Ap=πdlasinβ

(3)

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

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

CD=0.71+0.07RN

(4)

(5)

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

(6)

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

(7)

(8)

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

(9)

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

(10)

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

(11)

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

(12)

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

1.2 箔條間碰撞模型

為考慮箔條間碰撞,對計算空間進(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ù)目;Dp為箔條分段的等效直徑;Gab為箔條分段a和b的相對速度;go為徑向分布函數(shù);Δt為時間步長。

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

(14)

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

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

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

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

式中:下標(biāo)P代表平行,V代表垂直;ui,uj為箔條碰撞前碰點速度;vi0,vj0為箔條碰撞前質(zhì)心速度;Ii,Ij為坐標(biāo)轉(zhuǎn)換余弦矩陣;ωi0,ωj0為碰撞前箔條在體坐標(biāo)系中的角速度;ri,rj為碰點到箔條質(zhì)心的位置矢量;v為碰點最大壓縮形變時速度;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)動慣量;Ji,Jj為箔條在定坐標(biāo)系中的轉(zhuǎn)動慣量;e為碰撞過程中的恢復(fù)系數(shù),當(dāng)彈性碰撞時e=1,當(dāng)完全非彈性碰撞時e=0。

1.3 計算過程

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

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

2 數(shù)值計算和試驗結(jié)果對比

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

為了驗證模型的合理性,進(jìn)行了近地箔條拋撒試驗。針對本文試驗背景需求,設(shè)計了固體火箭發(fā)動機(jī)攜帶干擾彈高速拋撒的箭載試驗系統(tǒng)。點火后火箭加速至干擾彈發(fā)射所需的模擬速度時,采用延時點火炬點燃干擾彈發(fā)射藥,箔條集束被拋出。試驗時高速攝像機(jī)布置在彈道側(cè)方,用于記錄箔條集束分離、云團(tuán)擴(kuò)散的全過程。試驗中通過六通道同步觸發(fā)器控制火箭與高速攝像機(jī)的同步點火與觸發(fā)。圖3為箔條拋撒后不同時刻yz截面云團(tuá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中不同時刻云團(tuán)散布對比可知,當(dāng)所有箔條拋撒完畢,即t=50 ms時,箔條云團(tuán)呈較為整齊的錐形分布,y方向分布在4～13 m,z方向分布在1.5～3 m。t=100 ms時,箔條云已基本成型,x方向上分布在1～3 m,yz方向上的最大散布面積已達(dá)10 m2以上?？梢?箔條多體分離用時很短,能夠在100 ms內(nèi)迅速散開成型。300 ms后云團(tuán)分布基本穩(wěn)定,進(jìn)入沉降階段,箔條云團(tuán)整體以一定速度下沉。

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

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

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

3.1 數(shù)量分布

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

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

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

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

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

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

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

(23)

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

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

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

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

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

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

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

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

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

(24)

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

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

4 結(jié)束語

本文在剛體動力學(xué)及碰撞概率模型的基礎(chǔ)上,對箔條在牽連速度為270 m/s下的拋撒散布過程進(jìn)行了數(shù)值模擬,分析了云團(tuán)散布軌跡,并開展試驗,驗證了模型的準(zhǔn)確性,進(jìn)而對箔條云團(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ù)量分布情況基本不再隨時間變化。箔條集束能夠在0.3 s內(nèi)形成穩(wěn)態(tài)云團(tuán),云團(tuán)成熟后在軸向上的散布面積可達(dá)13 m2,徑向上的散布面積可達(dá)4 m2,云團(tuán)邊界值的發(fā)展趨勢在0.3 s時發(fā)生明顯變化。

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

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