火箭彈氣動力工程算法計算模型研究*

許邵杰，陳 雄，胡少青

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

0 引言

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

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

針對以上問題，文中利用AutoCAD自帶的AutoLISP語言編寫程序對火箭彈空氣動力圖表數據進行選取，并轉換為數據文本文件，再利用BASIC語言編寫插值程序對數據文本文件進行插值處理，進而完成火箭彈空氣動力模型的計算。

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

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

圖1 錐形頭部波阻系數

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

2 利用AutoLISP程序選取圖表數據

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

AutoLISP程序分為3個模塊：圖形水平校正(horizontal)、坐標系初始化(initialize)和選擇數據點(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" 圖像原點設定")

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

(setq xs(getreal" 輸入坐標原點x軸起點數值："))

(setq ys(getreal" 輸入坐標原點y軸起點數值："))

(prompt" 比例設定")

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

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

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

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

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

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

(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" 輸入數據文件名稱："))

(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結束："))

(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" 數據保存在")(princ fname)(princ"文件中")

(princ)

)

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

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

3 數據計算

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

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

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

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

圖2 錐形頭部波阻系數計算流程圖

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

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

已知兩點坐標值，通過點斜式線性插值公式確定連點連線段上任意點數值。

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

圖4 某模型計算結果

4 結論

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

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

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

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

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

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

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

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