一種基于Whitham方程的導(dǎo)彈沖擊波仿真優(yōu)化方法

劉軼煒, 陳蕙心, 蔡 猛, 吳思源, 陸子劼

(上海無線電設(shè)備研究所,上海 201109)

0 引言

物體超音速飛行,前方空氣被急劇壓縮產(chǎn)生音障,當(dāng)物體突破該音障時,由于前方空氣無法傳播而累積成為激波面。這種激波在波面前方形成一個突然的超壓邊界,在后方形成一個欠壓面,壓力隨時間變化的形狀如字母N,稱之為N形波。目前通常采用Whitham方程[1]對略超音速細(xì)長尖頭彈體進(jìn)行沖擊波仿真,對仿真結(jié)果中超壓振幅及周期等參數(shù)進(jìn)行分析,進(jìn)而研究飛行過程中沖擊波[2]對彈體的影響。彈體沖擊波仿真通常選用經(jīng)驗常數(shù)作為Whitham方程的形狀系數(shù)。由于不同類型導(dǎo)彈的形狀系數(shù)具有較大差異,采用上述方法得到的仿真結(jié)果與實際情況存在較大偏差,不具有普遍適用性。本文提出了一種基于Whitham方程的導(dǎo)彈沖擊波仿真優(yōu)化方法,根據(jù)彈體頭部橫截面形狀和彈體長寬比計算形狀系數(shù),提升超音速導(dǎo)彈彈體沖擊波仿真的準(zhǔn)確性。

1 Whitham方程與導(dǎo)彈形狀系數(shù)

以彈體頭部頂點為坐標(biāo)原點,沿彈軸方向自頭向尾為x軸正方向,建立oxyz直角坐標(biāo)系,彈體幾何模型如圖1所示。其中,S(x)為彈體的橫截面面積函數(shù),R為彈體最寬處半徑,H為彈體前端弧度部分長度,L為彈體整體長度。

利用Whitham方程可計算出距彈體一定距離的某點的沖擊波振幅及周期,同時彈體的形狀和尺寸等參數(shù)也被引入到Whitham方程中。基于Whitham方程改寫的彈體橫截面的函數(shù)表達(dá)式為

式中:s為橫截面位置;S″(x)表示對彈體橫截面面積S x（）二次求導(dǎo),即彈體橫截面的函數(shù)表達(dá)式。對式(1)中的變量s進(jìn)行積分,可得到Whitham積分慣性矩[3]

定義彈體形狀系數(shù)Fw,Fw與Whitham積分慣性矩的關(guān)系可表示為

利用形狀系數(shù)Fw可以計算某觀測點處沖擊波的振幅[4-5]

式中:P0為環(huán)境氣壓;γ為比熱容;Ma為觀測點處馬赫數(shù);r為觀測點與彈體的垂直距離。

以小型導(dǎo)彈為研究對象,若形狀系數(shù)

式中:D為導(dǎo)彈直徑。根據(jù)式(4)可得到常用的沖擊波振幅計算公式[6-7]

同時,依據(jù)Whitham方程,N形波的周期可以表示為

式中:c0為觀測點處的音速。

由式(5)和式(7)可知,沖擊波周期與彈長成負(fù)相關(guān)。由于天線罩參與后激波的形成,故上述負(fù)相關(guān)關(guān)系不符合氣動原理,因此給出了另一種關(guān)系式[8]

式(8)未引入形狀系數(shù)Fw,計算過程較簡單,但該公式無法反映導(dǎo)彈形狀對沖擊波的影響。

通過分析可知,當(dāng)導(dǎo)彈彈徑一定時,若彈體的圓柱體部分長度發(fā)生改變,沖擊波振幅基本保持不變,不同彈體長度下振幅的相對誤差僅為±2%,而周期受彈體長度影響改變較大。

對于圖1所示的彈體幾何模型,彈體的形狀可以用函數(shù)f(x)來表示,且基于 Whitham方程的F(s)為關(guān)于彈體橫截面的函數(shù),而ks為關(guān)于彈整體的函數(shù),故計算ks時需對F(s)沿彈軸,即x軸進(jìn)行積分求解。根據(jù)式(1),可以得到ks的表達(dá)式

圖1中S(x)=π(f(x))2,則S″(x)=f(x)。取彈長L=1 m,根據(jù)不同的彈體橫截面函數(shù)y=f(x),可計算得到不同的比例系數(shù)ks,如表1所示。

表1 不同彈體橫截面函數(shù)對應(yīng)的比例系數(shù)

由表1分析可知,不同彈體橫截面函數(shù)對應(yīng)的比例系數(shù)ks取值為1.00±0.15。

2 彈體形狀對沖擊波的影響

為了驗證彈體形狀對N形波參數(shù)的影響,設(shè)導(dǎo)彈飛行高度為海拔8.8 km,飛行馬赫數(shù)為1.7,使用CFD軟件[9]對不同尺寸及半徑長度比的彈體進(jìn)行沖擊波仿真。用els A工具計算彈體附近氣動尾跡的速度場及壓力場。在標(biāo)準(zhǔn)大氣中,平行于彈體的壓力分布會傳播到海平面,此條件下選用非線性傳播TRAPS代碼[10]計算N形波。仿真結(jié)果如表2所示。

表2 不同尺寸彈體氣動尾跡的速度場及壓力場CFD仿真結(jié)果

經(jīng)上述推導(dǎo)可知,綜合使用els A工具和TRAPS代碼計算所得的N形波是不對稱的。

設(shè)彈體的半徑長度比δ=R/L,則式(5)的形狀系數(shù)可以表示為

當(dāng)δ一定時,ΔP,ΔT與L0.75成正比。對照表2可以發(fā)現(xiàn):對于同一δ,超壓部分增長為L0.72,周期部分增長為L0.80;若給定彈體半徑R,其沖擊波振幅相對誤差僅為±2%。

3 非均勻大氣下沖擊波

根據(jù)式(3)和式(8)對彈體進(jìn)行沖擊波影響的仿真分析時,假設(shè)前提是大氣為均勻的。而大型導(dǎo)彈沖擊波發(fā)生在一定的海拔z處,通常采樣點在海平面處,故如采用上述理論對大型導(dǎo)彈的沖擊波影響進(jìn)行分析,存在一定的困難。

對于距彈體同樣距離r的兩個點M0和Mz,其中Mz位于與導(dǎo)彈相同海拔的位置z處,M0位于其他海拔位置,設(shè)這兩個點的大氣壓力分別為P0和Pz,音速分別為c0和cz。已知導(dǎo)彈飛行馬赫數(shù),由式(4)和式(7)可以得到Mz處的沖擊波振幅和周期為ΔPz和ΔTz。設(shè)M0處沖擊波的振幅和周期為ΔP0和ΔT0,若M0和Mz兩處的N形波聲能相等,則有

上述公式已在標(biāo)準(zhǔn)大氣環(huán)境下通過TRAPS傳播代碼進(jìn)行測試驗證。

結(jié)合某地空導(dǎo)彈[11]實際情況,對在地面提取的導(dǎo)彈飛行試驗沖擊波數(shù)據(jù)進(jìn)行仿真。設(shè)導(dǎo)彈飛行海拔8 km,飛行馬赫數(shù)1.7,推導(dǎo)可得其水平飛行時M0處的沖擊波振幅ΔP0以及周期ΔT0,表達(dá)式為

式中:CL為升力系數(shù);FwΔP和FwΔT分別為N形波振幅和周期對應(yīng)的的彈體形狀系數(shù)。

FwΔP和FwΔT可表示為

考慮到導(dǎo)彈飛行過程中的自體旋轉(zhuǎn),計算時將升力系數(shù)設(shè)為1。使用三種仿真方法對地空導(dǎo) 彈沖擊波數(shù)據(jù)進(jìn)行仿真,結(jié)果如表3所示。

表3 不同仿真方法下地空導(dǎo)彈的沖擊波數(shù)據(jù)

由表3可知,充分考慮彈體形狀和大氣的不均勻性,采用相應(yīng)的形狀系數(shù)計算沖擊波振幅和周期,可以明顯提高沖擊波仿真的準(zhǔn)確性,得到的沖擊波仿真結(jié)果相比Whitham方程更貼近CFD軟件仿真結(jié)果。根據(jù)研究的導(dǎo)彈目標(biāo)采用相應(yīng)的形狀系數(shù),可以明顯提高沖擊波仿真的準(zhǔn)確性;相應(yīng)的仿真變化式也能夠更準(zhǔn)確地反映不同大氣情況下的導(dǎo)彈沖擊波,提高了仿真的普遍適用性。

4 結(jié)論

本文針對采用Whitham方程進(jìn)行導(dǎo)彈沖擊波仿真時,由于經(jīng)驗形狀系數(shù)影響致使仿真結(jié)果與實際情況不符的問題,提出了一種基于Whitham方程的導(dǎo)彈沖擊波仿真優(yōu)化方法。該方法考慮不均勻大氣環(huán)境,結(jié)合彈徑和彈體長度比計算形狀參數(shù),并利用其進(jìn)行沖擊波參數(shù)仿真計算。與CFD軟件的仿真結(jié)果對照表明,該優(yōu)化方法大幅度減小了使用Whitham方程仿真導(dǎo)彈沖擊波時由于彈體形狀導(dǎo)致的誤差,在沒有精確的CFD模型的前提下,能夠針對不同類型的導(dǎo)彈給出更加準(zhǔn)確的沖擊波振幅以及周期仿真結(jié)果。