火箭弹气动力工程算法计算模型研究*

2012-12-10许邵杰胡少青

弹箭与制导学报 2012年3期

许邵杰，陈雄，胡少青

(南京理工大学机械工程学院，南京 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。可以看出，两者误差小于百分之一，计算速度快，精度有保证，完全满足工程设计精度要求。