許邵杰，陳 雄，胡少青

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

0 引言

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

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

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

1 火箭彈空氣動力計算模型與圖表

許多文獻在計算火箭彈的空氣動力參數(shù)時均采用了一些經(jīng)驗公式，但大部分計算仍需要大量曲線圖表。例如圖1所示，當(dāng)計算彈體錐形頭部波阻系數(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ù)碼相機等截取紙質(zhì)文獻中的曲線圖表并轉(zhuǎn)換為圖形文件，在AutoCAD系統(tǒng)中以光柵圖像形式插入，利用AutoLISP語言中的屏幕坐標(biāo)點輸入函數(shù)(getpoint)來讀取曲線上每一點的實際坐標(biāo)值，將讀取的數(shù)據(jù)點x和y坐標(biāo)分別乘以x向和y向比例系數(shù)并寫入數(shù)據(jù)文件中，通過編寫的程序循環(huán)讀取鼠標(biāo)點坐標(biāo)。只要用鼠標(biāo)連續(xù)點取曲線上點即可獲得點的坐標(biāo)值，將獲得的點坐標(biāo)值分別用(car)和(cadr)函數(shù)得到點的 x、y值，將x、y值分別乘以該方向的比例系數(shù)即可得到曲線上該點的實際坐標(biāo)值，并將計算坐標(biāo)實時寫入到數(shù)據(jù)文件中。

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

(defun c：horizontal(/pt1 pt2)

(prompt" 水平校正")

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

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

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

)

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

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

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

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

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

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

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

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

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

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

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

(setvar"userr2"(/(－ye ys)(－(cadr py)(cadr po))));計算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" 選取曲線起始點："))

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

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

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

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

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

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

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

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

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

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

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

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

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

4 結(jié)論

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

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

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

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

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

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

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

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