火箭彈氣動(dòng)力工程算法計(jì)算模型研究*

許邵杰，陳 雄，胡少青

(南京理工大學(xué)機(jī)械工程學(xué)院，南京 210094)

0 引言

許多火箭彈空氣動(dòng)力計(jì)算相關(guān)文獻(xiàn)中都存在有大量的圖表和試驗(yàn)曲線數(shù)據(jù)，而且在進(jìn)行計(jì)算過程中都大量多次進(jìn)行查圖線取值以及插值計(jì)算[1－2]。圖表數(shù)據(jù)有如下缺點(diǎn)：一是容易受主觀影響，讀值與插值過程中誤差大;二是耗費(fèi)時(shí)間長，且重復(fù)性不好，不利于編程計(jì)算和理論分析。

為克服上述缺點(diǎn)，可采取兩種途徑：一是離散圖表數(shù)據(jù)，用表格形式表示出各個(gè)點(diǎn)的數(shù)據(jù);二是進(jìn)一步對離散數(shù)據(jù)進(jìn)行擬合，得到一系列擬合曲線公式。第一種途徑數(shù)據(jù)量大，不直觀，處理不便;第二種途徑擬合曲線過程中存在擬合誤差，雖然可采用分段擬合、多次方程式擬合等方法減小誤差，但對最后的計(jì)算結(jié)果仍有一定影響[3]。

針對以上問題，文中利用AutoCAD自帶的AutoLISP語言編寫程序?qū)鸺龔椏諝鈩?dòng)力圖表數(shù)據(jù)進(jìn)行選取，并轉(zhuǎn)換為數(shù)據(jù)文本文件，再利用BASIC語言編寫插值程序?qū)?shù)據(jù)文本文件進(jìn)行插值處理，進(jìn)而完成火箭彈空氣動(dòng)力模型的計(jì)算。

1 火箭彈空氣動(dòng)力計(jì)算模型與圖表

許多文獻(xiàn)在計(jì)算火箭彈的空氣動(dòng)力參數(shù)時(shí)均采用了一些經(jīng)驗(yàn)公式，但大部分計(jì)算仍需要大量曲線圖表。例如圖1所示，當(dāng)計(jì)算彈體錐形頭部波阻系數(shù)時(shí)，需要根據(jù)彈體頭部長細(xì)比λn查其隨馬赫數(shù)Ma變化的圖形曲線。

圖1 錐形頭部波阻系數(shù)

確定飛行馬赫數(shù)Ma后，根據(jù)頭部長細(xì)比λn，從圖1中即可查出波阻系數(shù)的大小。

2 利用AutoLISP程序選取圖表數(shù)據(jù)

對于圖1所示波阻系數(shù)曲線模型，通常的手段都是劃分網(wǎng)格平均離散數(shù)據(jù)。文中利用掃描儀、數(shù)碼相機(jī)等截取紙質(zhì)文獻(xiàn)中的曲線圖表并轉(zhuǎn)換為圖形文件，在AutoCAD系統(tǒng)中以光柵圖像形式插入，利用AutoLISP語言中的屏幕坐標(biāo)點(diǎn)輸入函數(shù)(getpoint)來讀取曲線上每一點(diǎn)的實(shí)際坐標(biāo)值，將讀取的數(shù)據(jù)點(diǎn)x和y坐標(biāo)分別乘以x向和y向比例系數(shù)并寫入數(shù)據(jù)文件中，通過編寫的程序循環(huán)讀取鼠標(biāo)點(diǎn)坐標(biāo)。只要用鼠標(biāo)連續(xù)點(diǎn)取曲線上點(diǎn)即可獲得點(diǎn)的坐標(biāo)值，將獲得的點(diǎn)坐標(biāo)值分別用(car)和(cadr)函數(shù)得到點(diǎn)的 x、y值，將x、y值分別乘以該方向的比例系數(shù)即可得到曲線上該點(diǎn)的實(shí)際坐標(biāo)值，并將計(jì)算坐標(biāo)實(shí)時(shí)寫入到數(shù)據(jù)文件中。

AutoLISP程序分為3個(gè)模塊：圖形水平校正(horizontal)、坐標(biāo)系初始化(initialize)和選擇數(shù)據(jù)點(diǎn)(selectpoint)。具體代碼如下：

(defun c：horizontal(/pt1 pt2)

(prompt" 水平校正")

(setq pt1(getpoint" 選擇圖像水平線上第一點(diǎn)："))

(setq pt2(getpoint" 選擇圖像水平線上第二點(diǎn)："))

(command"rotate""all"""pt1"r"pt1 pt2 0)

)

(defun c：initialize(/po xs ys px xe py ye)

(prompt" 圖像原點(diǎn)設(shè)定")

(setq po(getpoint" 選擇曲線圖中坐標(biāo)原點(diǎn)："))

(setq xs(getreal" 輸入坐標(biāo)原點(diǎn)x軸起點(diǎn)數(shù)值："))

(setq ys(getreal" 輸入坐標(biāo)原點(diǎn)y軸起點(diǎn)數(shù)值："))

(prompt" 比例設(shè)定")

(setq px(getpoint" 選擇曲線圖中 x軸終點(diǎn)："))

(setq xe(getreal" 輸入x軸終點(diǎn)數(shù)值："))

(setq py(getpoint" 選擇曲線圖中 y軸終點(diǎn)："))

(setq ye(getreal" 輸入y軸終點(diǎn)數(shù)值："))

(setvar"userr1"(/(－xe xs)(－(car px)(car po))));計(jì)算x向比例系數(shù)

(setvar"userr2"(/(－ye ys)(－(cadr py)(cadr po))));計(jì)算y向比例系數(shù)

(setvar"userr3"xs)

(setvar"userr4"ys)

(command"move""all"""po(list 0 0))

(command"zoom""e")

(princ)

)

(defun c：selectpoint(/f1 k fname p0 pt)

(setq f1 nil k't)

(setq fname(getstring" 輸入數(shù)據(jù)文件名稱："))

(while k

(if(findfile fname)

(progn

(princ(strcat" "fname"文件已存在，重新輸入："))

(setq fname(getstring))

)

(setq k nil)

)

)

(setq f1(open fname"w"))

(setq p0(getpoint" 選取曲線起始點(diǎn)："))

(princ(+(getvar"userr3")(*(－(car p0)(getvar"userr3"))(getvar"userr1")))f1)

(princ""f1)

(princ(+(getvar"userr4")(*(－(cadr p0)(getvar"userr4"))(getvar"userr2")))f1)

(print(+(getvar"userr3")(*(－(car p0)(getvar"userr3"))(getvar"userr1"))))

(princ"")

(princ(+(getvar"userr4")(*(－(cadr p0)(getvar"userr4"))(getvar"userr2"))))

(while(setq pt(getpoint p0" 選取下一點(diǎn)或按Enter結(jié)束："))

(setq p0 pt)

(print(+(getvar"userr3")(*(－(car pt)(getvar"userr3"))(getvar"userr1")))f1)

(princ"")

(princ(+(getvar"userr4")(*(－(cadr pt)(getvar"userr4"))(getvar"userr2")))f1)

(print(+(getvar"userr3")(*(－(car pt)(getvar"userr3"))(getvar"userr1"))))

(princ"")

(princ(+(getvar"userr4")(*(－(cadr pt)(getvar"userr4"))(getvar"userr2"))))

)

(princ" 數(shù)據(jù)保存在")(princ fname)(princ"文件中")

(princ)

)

在圖形中選取數(shù)據(jù)點(diǎn)過程中，采用分段非平均取點(diǎn)，對于線性度較好的曲線段，例如圖1中各條曲線的上升段和下降段，選擇較少數(shù)據(jù)點(diǎn)，滿足圖形精度即可。對于圖形的轉(zhuǎn)折點(diǎn)及變化較大處，例如圖1中的馬赫數(shù)1.0到1.5附近曲線段，取點(diǎn)間距縮短，保證數(shù)據(jù)的可靠性。這樣處理在保證數(shù)據(jù)精度的前提下大大減小了數(shù)據(jù)處理量。實(shí)際使用時(shí)，為了提高數(shù)據(jù)的采集精度，在以交互方式選取曲線圖上坐標(biāo)點(diǎn)時(shí)，可利用 AutoCAD系統(tǒng)提供的視窗縮放功能(ZOOM)來提高顯示及采集精度(放大后沿曲線縱向中部取點(diǎn)，保證取點(diǎn)連線不超出曲線徑向范圍)。實(shí)際采集精度只與原曲線圖的繪制精度有關(guān)，幾乎不存在主觀采集誤差。表1為λn=2時(shí)的錐形波阻系數(shù)隨馬赫數(shù)Ma變化的數(shù)據(jù)。

表1 λn=2的錐形波阻系數(shù)隨馬赫數(shù)Ma變化的數(shù)據(jù)

3 數(shù)據(jù)計(jì)算

1)用input語句讀取經(jīng)過AutoLISP程序轉(zhuǎn)換的數(shù)據(jù)文本文件(txt格式)，并將x向坐標(biāo)數(shù)據(jù)和y向坐標(biāo)數(shù)據(jù)分別寫入兩個(gè)單獨(dú)數(shù)組(定義為數(shù)組A、數(shù)組B)。

2)根據(jù)所求的頭部外形(頭部長徑比)，選擇某一條或兩條曲線進(jìn)行插值計(jì)算。

3)根據(jù)飛行條件(輸入馬赫數(shù)Ma)，依據(jù)其在數(shù)組A中的位置對應(yīng)在數(shù)組B中進(jìn)行線性插值。

以圖1為例，求錐形頭部波阻系數(shù)的具體計(jì)算流程圖如圖2所示。

圖2 錐形頭部波阻系數(shù)計(jì)算流程圖

圖3 線性插值(一次插值)模型

數(shù)據(jù)處理過程中采用線性插值模型，如圖3所示，插值公式為：

已知兩點(diǎn)坐標(biāo)值，通過點(diǎn)斜式線性插值公式確定連點(diǎn)連線段上任意點(diǎn)數(shù)值。

以文獻(xiàn)[1]中算例為例，通過圖表計(jì)算得出，在攻角為5°、馬赫數(shù)為3.5時(shí)，該算例的升力系數(shù)為1.06，阻力系數(shù)為0.473，壓力中心為0.750。采用擬合公式的方法最小二乘多項(xiàng)式擬合，在四次方程形式下最大擬合誤差達(dá)到了2.3%。采用本方法進(jìn)行編程計(jì)算，結(jié)果如圖4所示，在相同條件下，得出升力系數(shù)為1.0539，阻力系數(shù)為0.47525，壓力中心為0.7480?？梢钥闯觯瑑烧哒`差小于百分之一，計(jì)算速度快，精度有保證，完全滿足工程設(shè)計(jì)精度要求。

圖4 某模型計(jì)算結(jié)果

4 結(jié)論

應(yīng)用本方法編制氣動(dòng)力計(jì)算程序，通過對若干文獻(xiàn)中的空氣動(dòng)力計(jì)算的大量圖表進(jìn)行處理，并與其他各種方法實(shí)際比較計(jì)算，可以得出如下結(jié)論：

1)相比于文獻(xiàn)中手工查圖線計(jì)算火箭彈氣動(dòng)力的傳統(tǒng)方法，速度更快，使用更方便，結(jié)果更精確;

2)相比于擬合公式計(jì)算，計(jì)算精度更高，誤差僅限于圖表的繪制精度和掃描的清晰程度;

3)采用分段非平均選取數(shù)據(jù)點(diǎn)，數(shù)據(jù)量大大減小，計(jì)算方便，且與原圖吻合程度高;

4)編程中采用線性插值方法，算法簡單，針對某些圖表具有多個(gè)限制條件，采用多次線性插值，不受插值順序影響，使用方便。

[1]臧國才，李樹常.彈箭空氣動(dòng)力學(xué)[M].北京：兵器工業(yè)出版社，1989.

[2]周長省，鞠玉濤，朱福亞，等.火箭彈設(shè)計(jì)理論[M].北京：北京理工大學(xué)出版社，2005.

[3]陳軍.火箭彈快速空氣動(dòng)力計(jì)算模型研究[J].彈箭與制導(dǎo)學(xué)報(bào)，2001，21(2)：45－47.