飞机气动及飞行性能计算

飞机气动及飞行性能计算

------ 课程设计报告

专业：飞行器设计与工程

班号：01011203

学号：**********

姓名：***

2016.3

 

第一章 预备知识

1.1 翼型的几何特性

参见上图：

中弧线    翼型内切圆中心的轨迹，在最前部内切圆（即决定前缘半径的圆）中

          心之前，则是由该内切圆中心至切点的半径线段

前缘      翼型中弧线的最前点

后缘      翼型中弧线的最后点

弦线      连接前缘与后缘的直线

弦长b(m)  前缘与后缘之间的直线线段长度

厚度c(m)  翼型最大内切圆的直径

相对厚度

最大厚度位置(m)   翼型最大内切圆的中心在翼型弦线上的投影至翼型前缘

                    的距离

最大厚度相对位置  

弯度f(m)  中弧线与弦线之间垂直于弦线的最大线段长度

相对弯度

最大弯度位置(m)   中弧线与弦线之间垂直于弦线的最大线段至翼型前缘的

                    距离

最大弯度相对位置  

前缘半径(m)  翼型最前部内切圆的半径

上弧线     从前缘到后缘，翼型的上部轮廓曲线，以y1=f1(x)表示

下弧线     从前缘到后缘，翼型的下部轮廓曲线，以y2=f2(x)表示

后缘角在后缘处上弧线和下弧线的二切线之间的角度

1.2 机翼的几何特性

参见上图：

飞机基准纵轴  可以取机身纵轴 

机翼基准平面  包含机翼弦线或外露翼根弦线与飞机对称平面垂直的平面

外露机翼      不包括穿越机身部分的机翼

毛机翼        包括穿越机身部分的机翼（穿越机身部分通常由左右机翼的前后

              缘的延长线所构成，如图所示）

机翼面积S(m2)  毛机翼在机翼基准平面上的投影面积

机翼展长(翼展) (m)  左右翼梢之间的距离

外露机翼面积(m2)  外漏机翼在机翼基准平面上的投影面积

毛机翼根弦长b0(m)  毛机翼的根部弦长

翼梢弦长b1(m)  机翼的梢部弦长

机翼局部弦长b(z)  机翼展向翼剖面的弦长，是展向位置z的函数

机翼平均几何弦长  

机翼平均气动弦长   

机翼展弦比     

机翼根梢比      

机翼后掠角  至前缘的距离为弦长一定百分比的点的连线与垂直于

                  弦线的平面之间的夹角。前缘后掠角、1/4弦线后掠角、后

                  缘后掠角分别以、、表示

机翼上反角  机翼1/4弦点线在垂直于弦线的平面上的投影与机翼

                  基准平面之间的角度，当翼梢高于翼根时角度为正

机翼安装角  机翼弦线与飞机基准纵轴之间的角度，若以飞机基准纵

                  轴为参考线，弦线的前缘点高于后缘点时角度为正

机翼剖面扭转角  机翼任意剖面的弦线与机翼基准平面之间的角度。若相

                      对于机翼基准平面，前缘点高于后缘点则角度为正

1.3 机身的几何特性

参见上图：

机身最大横截面积

机身当量最大直径   

机身长度   机身的前端面和后端面之间的距离

机身长径比   

头部（前体）长度  机身的柱形部分以前的长度

头部（前体）长径比   

尾部（后体）长度  机身的柱形部分以后的长度

尾部（后体）长径比   

第二章 飞机的基本情况和本文计算方案

2.1 飞机基本情况简介

F-4“鬼怪II”(Phantom II)是美国麦克唐纳公司(后合并为麦克唐纳·道格拉斯公司，现已并入波音)为美国海军研制的双座舰载战斗轰炸机，后来美国空军也大量采购，成为美国空海军60～70年代的通用主力战斗机，是美国第二代战斗机的典型代表，参加过越南战争和中东战争，也曾经是美国空军的“雷鸟”飞行表演队的表演用机。1956年开始设计，1958年5月原型机试飞，1961年10月开始交付海军使用，1963年11月进入空军服役。F-4不但空战能力好，对地攻击能力也很强，是60年代以来美国生产数量最多的战斗机。

    上图为越战时一架携带的F-4E。F-4其型别众多，其主要型别有：A型，舰队防空型，生产49架；B型，海军全天候型，共生产696架；RF-4B，B型的侦察型，生产46架；C型，空军战术型，共生产583架；RF-C，侦察型，共生产493架；D型，空军的战斗轰炸型，共生产825架；E型，空军制空型，共生产846架；F-4EJ，日本仿制型，共生产109架；RF-4E，西德等国使用的侦察型，共生产102架；F型，西德定购的单座型，共生产175架；J型，海军舰队防空型，共生产518架；K、M型，英国海军和空军的使用型，共生产174架；S型，J型的改进型；G型，空军的反雷达攻击型。至1981年停产为止，美国共生产了5195架各种型号的F-4飞机。使用国家包括美国、英国、日本、西德、伊朗、

希腊、埃及、西班牙、以色列以及韩国等。

F-4B飞机采用大后掠角小展弦比的机翼、其外侧可折叠部分有12度上反角，全动水平尾翼下反角为23度。F-4B装有两台J-79-GE-8型涡喷发动机。飞机上没有机炮，机腹部悬挂有4枚半埋式的麻雀IIIA型中距空空导弹，需要时，可在机翼上增挂2到4枚响尾蛇导弹。此外，他还有多种外挂方案以执行不同的任务。F-4B飞机装有AN/APQ-72雷达火控系统，其搜索距离为36～96公里，跟踪距离为10～40公里。

以下是F-4B飞机的部分原始数据：

表2-1  飞机的重量数据

方    案	载荷情况	重量（公斤）
空    重		12670
全    机		19040
正常载荷	4枚麻雀III	19740
超    载	4枚麻雀III＋2枚响尾蛇	190

表2-2  飞机的载油数据

油箱	最大储油量（公斤）	可用燃油量（公斤）
机身油箱	4260	4183
机翼油箱	1968	1943
飞机基本油量	6228	6126
机翼下副油箱	2 1148	2 1141
机身下副油箱	1857	1851
全机最大总油量	10379	10259

表2-3  飞机的几何数据

全机参数	机长	17.75米
	翼展	11.7米
	机高	4.95米
机翼参数	全翼面积	49.24米2
	外露翼面积	35.21米2
	翼展	11.7米
	前缘后掠角	52
	1/4弦线处后掠角	45
	上反角（外翼部分）	12
	上反角（内翼部分）	0
	安装角	0
	展弦比	2.79
	梯形比（根梢比）	5.48
	平均相对厚度	5.1
	平均气动弦长	5.02米
	副翼面积	2 1.32米2
	副翼偏度	向上0  向下30
	最大厚度线处后掠角	41 29’
水平尾翼	全面积	8.31米2
	外露面积	6.60米2
	翼展	5.00米
	前缘后掠角	42.5
	下反角	23
	最大厚度线处后掠角	33 55’
	安装角	0
	展弦比	2.93
	外露翼梯形比	4.46
	相对厚度	3.62
	平均气动弦长	1.711米
	平尾偏角范围	向上14.5  向下20
	尾臂（1/2bA处至飞机重心）	7.19～7.63米
	尾容量	0.192

垂直尾翼	面积	5.52米2
	翼展	1.75米
	前缘后掠角	65.5
	展弦比	0.51
	梯形比	4.1
	平均气动弦长	3.44
	相对厚度	3.61
	尾臂（1/2bA处至飞机重心）	6.74～7.18米
	尾容量	0.151
	方向舵面积	1.28米2
	方向舵偏角范围	左右28
	最大厚度线处后掠角	52 47’
机身	全长	17.75米
	最大当量直径	2.13米
	最大截面积	3.57米2
	进气口面积	2 0.405米2
	头	4.9米
	头细比	2.93
	柱段长	8.50米
	柱段长细比	3.94
	尾段长	4.35米
	尾段长细比	2.04
	长细比	8.34
	减速板面积	2 1.04米2
	尾段收缩比	0

表2-4  J-79-GE-8最大状态可用推力表（单台）

	0.2	0.4	0.6	0.8	1.0	1.2	1.4
0km	3700	3590	3600	3680	3820	3980
5km	2310	2310	2410	2560	2700	2910	3160
8km		10	1760	1900	2100	2300	2520
10km				1520	1700	10	2080
11km				1360	1550	1710	1900

表2-5  J-79-GE-8全加力状态可用推力表（单台）

	0.4	0.6	0.8	1.0	1.2	1.4	1.6	1.8	2.0	2.2	2.4
0km	6500	7000	7700	8650	9800	10750
5km	4150	4500	5200	5950	6750	7550	8200
8km	2900	3300	3750	4450	5250	6100	6800	7200
10km			3000	3600	4250	5050	5950	6550	7000
11km			2650	3200	3800	4550	5500	6250	6750	6900	6950

注：1).11km以上的推力数据可按公式进行计算，

    2).表4和表5中推力P的单位为公斤(kg)。

2.2 本文计算方案

为高效实现该飞机气动及飞行性能计算，本文采用Matlab软件进行编程，对数据批量处理。首先根据第二章的飞机原始数据建立飞机基本参数数据库，进而通过程序实现所需参数的计算。对于涉及查图得到参数的内容时，采用合适的插值方法（如径向基函数），完成图中曲线拟合并求取未知点坐标。

第三章 飞机气动特性估算

飞机升阻特性是飞机最为重要的原始数据之一，在性能计算、飞行仿真等方面必不可少。

在飞机设计过程中，特别是方案论证或方案设计初期，气动布局等总体参数通常是变化的，翼型等参数尚未完全确定，因此计算精确的气动数据较为困难。通常采用工程方法进行气动估算，以获得进一步计算分析的原始参数。另外对于国外设计的飞机，由于无法得到精确的翼型等外形参数，也只能够对其进行气动估算以获得其气动参数。

3.1 升力特性估算

作用在飞机上的升力可以表示为：

其中

            升力系数

            机翼参考面积

            动压

            

对于没有增升装置的对称翼型，升力系数可以表示为：

        升力线斜率

            迎角

对于非对称翼型，升力系数可以表示为：

            零升迎角，取决于机翼的弯度等特性

从上式可以看出，描述飞机升力特性的参数主要是和。

飞机的机翼升力特性如图3-1所示：

图3-1 机翼升力特性

通过下面的工程估算方法可以计算得到飞机的升力线斜率，按照机翼、平尾、机身三部分分别求解并叠加。

3.1.1 单独机翼升力估算

对于单独的机翼，其升力线斜率可以表示为以下参数的函数：

其中

            展弦比

        1/2弦线的后掠角

            机翼相对厚度

            尖削比， 

            根梢比或称梯形比

其函数关系较为复杂，可以由指导资料中图3曲线查出。

有时机翼的几何参数数据给出机翼的前缘后掠角，则1/n弦线的后掠角可以由下式求出：

其中

            前缘的后掠角

对于大展弦比的后掠翼来说，其升力线斜率可以表示为：

其中

            翼型效率，可取0.95

估算过程：

因为机翼展弦比=2.79，属小展弦比，故采用查图3的方式求升力线斜率。

得到：

对指导资料中图3(c)进行处理，令y轴左侧为负值，以便于插值。采用分段插值方法，将图3(c)分为两段，分界处为。第一段选择8个点进行RBF插值拟合曲线，径向基函数选为逆多二次，分散度为各个插值点的平均距离2.66；第二段选择7个点进行RBF插值拟合曲线，径向基函数选为逆多二次，分散度为2.90。插值后函数图像如图3-2所示：

图3-2 机翼升力线斜率计算

得到的各个马赫数下单独机翼升力线斜率为：

Ma
0.4	0.0178	0.0497
0.6	0.0185	0.0515
0.8	0.0200	0.0557
1.2	0.0204	0.0568
1.4	0.0179	0.0500
1.6	0.0162	0.0451
1.8	0.0167	0.0467
2.0	0.0146	0.0407
2.2	0.0129	0.0360

3.1.2 机身升力估算

机身升力主要由头部及尾部两部分构成，对于圆柱形状的机身，有：

其中

        机身的升力线斜率

        头部产生的升力线斜率

            尾部收缩比

            底部面积，若尾部形状为锥形，则底部面积为零

            机身面积，即尾部的最大面积

            修正系数，取决于雷诺数、马赫数、尾部形状、尾翼布局等参数

            可取0.15～0.20。例如可取0.17或0.18

可按照下式查资料中图4曲线得出：

其中

            头细比

            机身圆柱部分长细比

如果机身截面形状为椭圆形，则其升力线斜率按照圆柱形进行修正：

其中

            机身最大截面的宽度

估算过程：

同样，对资料中图4进行处理，令y轴左侧为负值，以便于插值。选择10个点进行RBF插值拟合曲线，径向基函数选为逆多二次，分散度为各个插值点的平均距离1.1。插值后函数图像如图3-3所示：

图3-3 具有锥形头部旋成体的头部升力计算曲线

其他相关量如图：

Ma
0.4	0.0360	0.0301	0.0300
0.6	0.0363	0.0303	0.0303
0.8	0.0369	0.0310	0.0309
1.2	0.0463	0.0404	0.0403
1.4	0.0490	0.0430	0.0429
1.6	0.0507	0.0448	0.0447
1.8	0.0519	0.0459	0.0458
2.0	0.0526	0.0467	0.0466
2.2	0.0532	0.0472	0.0471

3.1.3 翼身组合体的升力估算

现代超音速战斗机的直径对翼展的比值可以达到0.3～0.5，在这种情况下用单独的机翼代替翼身组合体就会带来很大误差。

通常计算翼身组合体的升力如下：

其中

        机翼外露部分（外露翼）升力，再考虑机身对机翼升力影响的修正

        单独机身的升力

这里，我们忽略了机翼对机身升力的影响。

外露翼升力系数在考虑机身的影响后，可以写为：

其中

        外露翼部分升力线斜率，按单独机翼的方法（见第一节）进行计算，所有参数均按照外露翼取值

    机翼升力系数，考虑了机身的影响，参考面积按照外露翼面积

f            修正系数

d            机身直径

l            翼展

估算过程：

将飞机看做翼身组合体，则修正系数：

f = 1.4951

得到的外露翼升力系数：

Ma
0.4	0.0744
0.6	0.0770
0.8	0.0833
1.2	0.0850
1.4	0.0748
1.6	0.0675
1.8	0.0698
2.0	0.0609
2.2	0.0538

3.1.4 尾翼升力估算

尾翼分为水平尾翼和垂直尾翼，只有水平尾翼产生升力。尾翼升力线斜率首先按照单独机翼的升力线斜率估算方法，计算出单独尾翼的升力线斜率，再进行修正，主要修正下洗和阻滞。

其中

        按单独尾翼计算的升力线斜率

                尾翼处的气流下洗角，近似认为等于机翼处的气流下洗角

                气流阻滞系数，可根据尾翼布局按照资料中表6确定

对于三角形机翼后气流下洗角的计算可以通过资料中图5由和计算，对于根梢比为无穷大的、后缘具有不大后掠角的机翼，可以采用同样方法确定。对于梯形机翼（）产生的下洗角可以对三角形机翼的下洗进行修正：

            不考虑机翼根梢比的下洗系数

A                尖削比对下洗的影响系数，通过资料中图6确定

            由单独机翼计算的参数

估算过程：

尾翼中仅水平尾翼产生升力。

因为尾翼展弦比=2.93，属小展弦比，故采用查图3的方式求升力线斜率。

得到：

同样，对资料中图3(b)进行处理，令y轴左侧为负值，以便于插值。采用分段插值方法，将图3(b)分为两段，分界处为。

第一段选择8个点进行RBF插值拟合曲线，径向基函数选为逆多二次，分散度为4.5；第二段选择7个点进行RBF插值拟合曲线，径向基函数选为逆多二次，分散度为各个插值点的平均距离1.。插值后函数图像如图3-4所示：

图3-4 尾翼升力线斜率计算

得到的各马赫数下单独尾翼升力线斜率为：

Ma
0.4	0.0028	0.0082
0.6	0.0031	0.0092
0.8	0.0039	0.0114
1.2	0.0210	0.0615
1.4	0.0218	0.0639
1.6	0.0176	0.0516
1.8	0.0151	0.0443
2.0	0.0132	0.0387
2.2	0.0117	0.0343

F-4B的机翼不是三角形，因此需要通过对三角形机翼下洗进行修正。首先对图5进行处理，令y轴左侧为负值，以便于插值。通过资料中图5计算下洗系数。选择10个点进行RBF插值拟合曲线，径向基函数选为逆多二次，分散度为各个插值点的平均距离3.67。插值后函数图像如图3-5所示：

图3-5 确定三角机翼后面气流下洗角的曲线

对资料中图6进行处理，令y轴左侧为负值，以便于插值。选择9个点进行RBF插值拟合曲线，径向基函数选为逆多二次，分散度为6.5。插值后函数图像如图3-6所示：

图3-6 确定参数A所用的曲线

得到的下洗系数、影响系数A及气流下洗角，进而得到的尾翼升力线斜率为：

Ma		A
0.4	34.8728	0.9624	0.5983	0.0030
0.6	34.8044	0.9604	0.6170	0.0032
0.8	34.5616	0.9521	0.6568	0.0035
1.2	32.1883	0.8446	0.5537	0.0247
1.4	29.32	0.8274	0.4359	0.0324
1.6	26.7059	0.8211	0.3548	0.0300
1.8	25.3313	0.8168	0.34	0.0261
2.0	24.8205	0.8131	0.2945	0.0246
2.2	24.97	0.8128	0.2611	0.0228

3.1.5 合升力线斜率计算

以上计算的各个部件的升力系数其参考面积均为各自的参考面积，例如机身的参考面积一般采用机身截面的面积，机翼的参考面积为外露翼部分面积，尾翼的参考面积为尾翼外露面积，这样为求得合升力系数，必须对其参考面积进行转化后再叠加，其计算公式如下：

其中

        外露翼面积

            机身截面积

            平尾面积

            全翼面积

估算过程：

首先根据资料中图2的升力线斜率与马赫数关系，插值得到待求的升力线斜率：

根据图2，选择7点插值拟合曲线，用RBF插值方法，径向基函数选为逆多二次，选择分散度为0.45。输出得到下表：

Ma
0.4	0.0726
0.6	0.0749
0.8	0.0817
1.2	0.0956
1.4	0.0854
1.6	0.0741
1.8	0.0629
2.0	0.0524
2.2	0.0452

按照公式叠加得到的升力线斜率为：

Ma		（图2）
0.4	0.0558	0.0726
0.6	0.0578	0.0749
0.8	0.0624	0.0817
1.2	0.0678	0.0956
1.4	0.0620	0.0854
1.6	0.0566	0.0741
1.8	0.0577	0.0629
2.0	0.0510	0.0524
2.2	0.0457	0.0452

将两种结果进行比较。如图3-7所示：

图3-7 计算的不同马赫数下升力系数（This Paper）和资料中（Ref.fig）的对比

从图3-7的对比可知，计算的升力系数基本与资料中的变化趋势一致，但在Ma>1.8时存在一个较大的上升，如果不考虑计算误差，目前推断可能是资料中图3的曲线不连贯所导致的。此外，计算的升力线斜率普遍偏小。

3.2 升阻极曲线的估算

作用在飞机上的气动阻力可以表示为

其中阻力系数可以表示为

或

其中

            零升阻力系数

A                诱导阻力因子

阻力系数与升力系数的关系可以用极曲线表示，图7给出了极曲线的两种形式。

图3-8 极曲线的两种形式，左图为无弯度机翼，右图为有弯度机翼

3.2.1 亚音速零升阻力估算

亚音速范围内，飞机的零升阻力主要由表面摩擦阻力和气流分离引起的压差阻力组成。

其中

            摩擦阻力系数

            压差阻力系数

3.2.1.1 全机摩擦阻力估算

其中

、、、分别为机翼、机身、平尾、    垂尾（立尾）的厚度修正系数

            机身浸润面积

                垂尾（立尾）面积

、、、分别为机翼、机身、平尾、    垂尾（立尾）的摩擦系数，它们与表面附面层状态、沿表面压力分布梯度及表面粗糙情况有关，同时也与基于各部件特征长度的雷诺数有关。当飞机在大气中飞行时，基于各部件特征长度的飞行雷诺数通常是相当大的，加上由于工艺水平等原因，飞机表面不可能做到理想的光滑，因此可以把飞机附面层近似看成是全湍流附面层。

对于光滑平板，具有全湍流附面层的表面摩擦系数可以用下面的半经验公式表示：

其中

Re            基于各部件特征长度计算的雷诺数

受工艺水平所限，飞机不可能做到理想的光滑，诸如铆钉头、螺钉头、缝隙、蒙皮台阶以及表面喷漆、划伤等因素，使得飞机相当粗糙。另外，飞机上还常有如天线、空速管、通风口鼓包等附加物。对此，在方案设计阶段通常用一个系数来考虑这些由表面粗糙和附加物产生的对阻力的影响，这就是§3.2.1的公式中1.1的来历。对于轻型战斗机，也可以用1.15。

厚度修正系数、、的计算公式如下，考虑了马赫数对摩擦影响的修正

其中

            翼型最大厚度线的弦向位置，无量纲

        最大厚度线的后掠角

对于机身，的计算公式如下

其中

            机身长度

            机身直径

机身的浸润面积计算公式如下

其中

、、            头部、尾部、柱段长度

估算过程：

按照高度为0km选择空气密度与声速。各参数计算如下：

机翼、机身、平尾、垂尾的摩擦系数为： 

Ma
0.4	0.2759e-03	0.2309e-03	0.3245 e-03	0.2917 e-03
0.6	0.2602e-03	0.2186e-03	0.3049 e-03	0.2748 e-03
0.8	0.2498e-03	0.2104e-03	0.2920 e-03	0.2636 e-03
1.2	0.2361e-03	0.1996e-03	0.2750 e-03	0.24 e-03
1.4	0.2312e-03	0.1957e-03	0.26 e-03	0.2436 e-03
1.6	0.2270e-03	0.1924e-03	0.2638 e-03	0.2391 e-03
1.8	0.2235e-03	0.15e-03	0.2594 e-03	0.2352 e-03
2.0	0.2203e-03	0.1870e-03	0.2555 e-03	0.2319 e-03
2.2	0.2175e-03	0.1848e-03	0.2521 e-03	0.22 e-03

修正系数：

Ma
0.4	0.2375	1.1245	0.2444	0.2237
0.6	0.4928	1.1245	0.5071	0.41
0.8	0.8271	1.1245	0.8511	0.7790
1.2	1.7160	1.1245	1.7658	1.6163
1.4	2.28	1.1245	2.3305	2.1331
1.6	2.8801	1.1245	2.9607	2.7127
1.8	3.5603	1.1245	3.6637	3.3533
2.0	4.3038	1.1245	4.4287	4.0536
2.2	5.1092	1.1245	5.2576	4.8122

由此得到的全机亚音速摩擦阻力：

Ma
0.4	0.00066
0.6	0.00078
0.8	0.00093

3.2.1.2 亚音速压差阻力估算

在计算压差阻力时，由于机翼及尾翼的压差阻力非常小，所以只考虑机身的压差阻力。飞机在超音速飞行时，压差阻力实际上就是波阻，所以不单独计算压差阻力。压差阻力可以按照下式分为头部阻力、尾部阻力、底部阻力、附加阻力四部分。

其中

        头部阻力系数，取决于头细比、马赫数，见资料中图8。

        尾部阻力系数，可以通过资料中图9由尾细比、收缩比、马赫数确定。(由于纵坐标没有刻度，故此项可暂时忽略)

        底部阻力系数，通常超音速战斗机发动机安装在尾部，所以此项为0。

        附加阻力系数，通常取0.007～0.01。

估算过程：

首先计算头部阻力系数，拟合资料中图8的亚音速区：

选择5个点进行RBF插值拟合曲线，径向基函数选为高斯基函数，分散度为4。插值后函数图像如图3-9所示：

图3-9 抛物线母线头部的阻力系数与马赫数的关系

尾部阻力系数，因为没有刻度，故暂时忽略。

底部阻力系数为0；

附加阻力系数取0.009。

最后的亚音速压差阻力：

Ma
0.4	-0.00148
0.6	-0.00153
0.8	-0.00119

综上，亚音速的零升阻力为：

Ma
0.4	-0.000
0.6	-0.00084
0.8	-0.00028

3.2.2 超音速零升波阻估算

飞机某一部件在局部马赫数超过1.0时，就会有波阻的存在，这个飞行状态的马赫数称之为临界马赫数，计算飞机的波阻时，必须首先确定临界马赫数。

3.2.2.1 临界马赫数的确定

机翼临界马赫数主要取决于机翼剖面形状、展弦比、后掠角等因素：

其中

        临界马赫数

        机翼剖面的临界马赫数，通过资料中图10，由机翼升力系数、相对厚度和翼型最大厚度线的弦向位置所决定。

        展弦比对临界马赫数的影响，由资料中图11根据零升临界迎角查得。

    后掠角对临界马赫数的影响，由资料中图11根据零升临界迎角查得。

估算过程：

用升力与马赫数关系推导出马赫数与升力的关系。

第一部分，机翼剖面临界马赫数：

通过计算最大厚度线的弦向位置，对资料中图10(b)进行拟合，注意此时自变量为升力系数。选择8个点进行RBF插值拟合曲线，径向基函数选为高斯基函数，分散度为0.23。拟合后函数图像如图3-10所示。

图3-10 剖面临界马赫数与升力系数的关系

第二部分，展弦比的影响：

对资料中图11左图拟合，选择7个点进行RBF插值拟合曲线，径向基函数选为高斯基函数，分散度为0.2。拟合后函数图像如图3-11所示。

图3-11 展弦比与后掠角对临界马赫数影响曲线

第三部分，后掠角的影响：

对资料中图11右图曲线拟合，选择7个点进行RBF插值拟合曲线，径向基函数选为高斯基函数，分散度为0.23。拟合后函数图像如图3-12所示。

图3-12 展弦比与后掠角对临界马赫数影响曲线

将三个结果叠加，得到的临界马赫数为：

0.1	0.7819
0.2	0.6374
0.3	0.5315
0.4	0.4977
0.5	0.4204
0.6	0.3811

确定了临界马赫数之后，就可以将飞行状态按马赫数分为三个阶段：、、，其中前者可以按亚音速方法处理，后者则完全按超音速方法处理，对于的情况，通常难以进行估算，为获得其数据可以利用图解法，由另外两种情况计算结果曲线进行光滑过渡而得出。

3.2.2.2 M>1时零升阻力系数

实践证明，超音速摩擦阻力的计算可以使用前面的亚音速摩擦阻力计算方法。在超音速情况下，摩擦阻力几乎与剖面形状无关，不需要进行剖面形状修正，因此在厚度修正系数表达式中可以认为相对厚度值为零。

超音速零升阻力的另一部分是零升波阻，零升波阻可以表示为各部件波阻之和：

其中

        零升波阻

、、、分别为机翼、机身、平尾、垂尾的波阻系数

单独机翼的波阻与飞行马赫数、机翼剖面形状和平面形状有关。资料中图12以组合参数形式给出了计算机翼波阻的工作曲线。每一张曲线对应菱形剖面和给定的尖削比。图中点划线是利用超音速线性理论计算的结果，而实线是根据实验数据整理的结果。由图可见，当和时两组曲线有较大差别，参数对波阻系数有显著影响。在做机翼波阻时，宜取实线值。平尾与垂尾的波阻系数也可以按照此理论进行计算。

对于非菱形机翼，其波阻计算式为

其中

        菱形剖面机翼的波阻系数，由资料中图12查得

            非菱形剖面的修正因子，由资料中表7确定

            由机翼最大厚度线的后掠角所确定的修正因子，由资料中图13确定

飞机机身的波阻系数分别由头部波阻、尾部波阻和头部对尾部的干扰阻力所组成：

其中

        头部波阻，可以查资料中图14得出

        尾部波阻，可以查资料中图15得出

        头部对尾部的干扰阻力，当柱段长度大于2倍直径时，认为头部对尾部的干扰很小，其阻力可以忽略不计。

估算过程：

超音速摩擦阻力，与亚音速摩擦阻力计算方法相同。通过之前公式可以得到：

Ma
1.2	0.00128
1.4	0.00152
1.6	0.00178
1.8	0.00206
2.0	0.00237
2.2	0.00270

接下来计算零升波阻。选择翼型形状为四角形，从而得到机翼、平尾、垂尾的非菱形剖面修正因子为K = 1.0423，1.2352和1.66。

计算菱形剖面机翼的波阻系数，则需要拟合资料中图12的曲线。具体参数：

机翼：选择11个点进行RBF插值拟合曲线，径向基函数选为高斯基函数，分散度为样本点之间的平均距离2.9670。拟合后函数图像如图3-13(a)所示。

平尾：选择11个点进行RBF插值拟合曲线，径向基函数选为高斯基函数，分散度为样本点之间的平均距离3.0693。拟合后函数图像如图3-13(b)所示。

垂尾：选择11个点进行RBF插值拟合曲线，径向基函数选为高斯基函数，分散度为样本点之间的平均距离0.5406。拟合后函数图像如图3-13(c)所示。

图3-13(a) 菱形机翼波阻计算图

图3-13(b) 菱形平尾波阻计算图

图3-13(c) 菱形垂尾波阻计算图

机翼最大厚度线后掠角修正因子，同样对资料中图13曲线进行拟合。选择11个点进行RBF插值拟合曲线，径向基函数选为高斯基函数，分散度为样本点之间的平均距离3.2020。拟合后函数图像如图3-14所示。

图3-14 菱形垂尾波阻计算图

飞机机身的波阻系数求解：

头部波阻，查资料中图14，选择10个点进行RBF插值拟合曲线，径向基函数选为线性基函数，分散度为样本点之间的平均距离0.3481。拟合后函数图像如图3-15所示。

图3-15 尖拱形头部跨、超音速波阻系数

尾部波阻，查资料中图15，选择10个点进行RBF插值拟合曲线，径向基函数选为逆多二次基函数，分散度为样本点之间的平均距离0.79。拟合后函数图像如图3-16所示。

图3-16 拱形尾部跨、超音速波阻系数

最终得到的零升波阻：

Ma
1.2	0.0321
1.4	0.0299
1.6	0.0275
1.8	0.0251
2.0	0.0230
2.2	0.0212

3.2.3 亚音速升致阻力估算

飞机在正常飞行状态下，升力主要由机翼产生，因此，在对飞机进行气动估算时，可以近似采用机翼的升致阻力代替全机的升致阻力。

飞机升致阻力可以由升致阻力因子所描述，对于升力沿展向椭圆分布的机翼，。实际机翼升力沿展向分布受机翼平面形状影响：

其中

            奥斯瓦德因子，是机翼展弦比和后掠角的函数

对于直机翼

对于后掠翼

通常情况，升致阻力系数可能无法表示为升致阻力因子的形式，则其升致阻力系数可以表示为：

估算过程：

飞机在正常飞行状态下，升力主要由机翼产生，因此，在对飞机进行气动估算时，可以近似采用机翼的升致阻力代替全机的升致阻力。

对升致阻力因子描述的升致阻力，需要计算升致阻力因子：

对于后掠翼，计算的升致阻力因子为：

升致阻力系数与不同的升力系数有关，可以表示为：

0.1	0.0108
0.2	0.0458
0.3	0.1093
0.4	0.2069
0.5	0.3456
0.6	0.5346

3.2.4 超音速升致阻力估算

超音速情况下，机翼的升致阻力系数可以表示为

其中

            前缘吸力对升致阻力的影响，按资料中图16

            修正系数，按资料中图17确定

上式只适应于小迎角范围，后一项表示前缘吸力对升致阻力的影响，可以用于有限根梢比的的机翼。经验表明，吸力的实际数值比理论值要低得多，故引入修正系数k。

升力线斜率是马赫数的函数，所以升致阻力系数不仅随变化，同时也随马赫数变化，计算时必须加以考虑。

估算过程：

超音速情况下，机翼的升致阻力系数与升力系数及升力线斜率都有关系，因此需要计算不同Ma和升力系数下对应的升致阻力系数。其中，需要从资料中图16中求处前缘吸力对升致阻力的影响，以及从资料中图17中求出修正系数。

图16，选择10个点进行RBF插值拟合曲线，径向基函数选为逆多二次基函数，分散度为样本点之间的平均距离0.3047。值得注意的是，其中会出现横轴大于1的情况，对此将最终结果修正为0。拟合后函数图像如图3-17所示。

图3-17 计算的曲线

图17，选择9个点进行RBF插值拟合曲线，径向基函数选为逆多二次基函数，分散度为样本点之间的平均距离0.2595。拟合后函数图像如图3-18所示。

图3-18修正系数k

得到的结果如下（横向为升力系数从小到大，纵向为马赫数从小到大）：

	0.1	0.2	0.3	0.4	0.5	0.6
1.2	0.0023	0.0094	0.0218	0.0395	0.0625	0.0908
1.4	0.0026	0.0106	0.0244	0.0438	0.0690	0.1000
1.6	0.0030	0.0122	0.0276	0.0491	0.0769	0.1108
1.8	0.0030	0.0121	0.0272	0.0484	0.0757	0.10
2.0	0.0034	0.0137	0.0308	0.0547	0.0854	0.1230
2.2	0.0038	0.0153	0.0343	0.0610	0.0954	0.1373

3.2.5 不同马赫数下的升阻极曲线

从上述结论中，可以得到不同马赫数的升阻极曲线，如图3-19所示：

(a) Ma=0.4

(b) Ma=0.6

(c) Ma=0.8

(d) Ma=1.2

(e) Ma=1.4

(f) Ma=1.6

(g) Ma=1.8

(h) Ma=2.0

(i) Ma=2.2

图3-19 不同Ma下的升阻极曲线

3.3 结果汇总

各马赫数下的升力线斜率及升力特性曲线见1.1.5。

各升力系数下的临界马赫数见1.2.2.1。

各马赫数，各升力系数下的阻力系数见1.2.4，不同马赫数下的升阻极曲线见1.2.5。

第四章 飞机基本飞行性能计算

4.1 速度-高度范围

飞机速度－高度范围包括最大飞行速度（马赫数）、最小飞行速度（马赫数）、静升限等速度、高度性能，是飞机基本飞行性能的重要组成部分，也是飞机的主要战术技术指标。这些性能指标的计算通常由飞机质点运动方程出发，采用简单推力法进行计算。

飞机定直平飞时的运动方程：

忽略 、 p的影响（），可简化为：

首先根据运动方程计算平飞需用推力Ppx，在给定飞行状态（H，M）下，由法向力方程Y=G，可得飞机在此飞行状态下的升力系数：

由极曲线可根据求得，则可以计算平飞需用推力：

根据飞机在同一高度不同速度下的平飞需用推力计算结果可以绘制出平飞需用推力曲线，在同一张图上，同时绘制飞机在最大推力状态下或全加力状态下的推力（称之为可用推力），则构成推力曲线图（资料中图18），飞机在此飞行状态下的飞行速度范围就可以通过平飞需用推力与可用推力曲线的交点确定。在某种推力状态下（最大或加力），需用推力曲线与可用推力曲线左侧的交点决定了最小飞行马赫数，右侧的交点决定了最大飞行马赫数。

飞机的最小平飞马赫数还取决于失速迎角等因素的，以上采用简单推力法所确定的只是由推力所的最小平飞马赫数，实际上略大于真实值。

根据不同高度下的飞行马赫数范围，绘制H-Ma曲线，则构成了飞行包线（资料中图19）。注意，前面的计算只考虑了推力，实际上飞行包线的边界还受到失速迎角（气动边界）、最大飞行马赫数（气动加热边界）、最大动压（结构强度边界）等因素。

估算过程：

在给定飞行状态（H，M）下，由法向力方程得到飞机的升力系数：

	0.4	0.6	0.8	1.2	1.4	1.6	1.8	2.0	2.2
0	0.3462	0.1539	0.0865	0.0385	0.0283	0.0216	0.0171	0.0138	0.0114
5	0.94	0.2886	0.1623	0.0722	0.0530	0.0406	0.0321	0.0260	0.0215
10	1.3269	0.58	0.3317	0.1474	0.1083	0.0829	0.0655	0.0531	0.0439
13	2.1245	0.9442	0.5311	0.2361	0.1734	0.1328	0.1049	0.0850	0.0702
16	3.4098	1.5155	0.8524	0.37	0.2783	0.2131	0.1684	0.13	0.1127
20	6.4074	2.8477	1.6019	0.7119	0.5231	0.4005	0.31	0.2563	0.2118

进一步求解各个状态下的阻力系数，通过3.2.1中公式进行计算。

	0.4	0.6	0.8	1.2	1.4	1.6	1.8	2.0	2.2
0	0.0148	0.0023	0.0007	0.0350	0.0318	0.0295	0.0273	0.0255	0.0239
5	0.0543	0.0101	0.0032	0.0392	0.0327	0.0301	0.0277	0.0257	0.0241
10	0.2298	0.0447	0.0141	0.0588	0.0370	0.0329	0.0292	0.0268	0.0249
13	0.5904	0.1160	0.0367	0.1007	0.0461	0.0385	0.0325	0.0291	0.0266
16	1.5222	0.3000	0.0949	0.2122	0.0701	0.0530	0.0409	0.0351	0.0310
20	5.3774	1.0616	0.3359	0.5943	0.1716	0.1133	0.0757	0.0596	0.04

由此得到不同飞行高度和不同马赫数下的平飞需用推力。在实际计算中，由于之前的误差，在平飞需用推力的计算时，在Ma=1.2出现一个较大的偏差，使曲线的趋势不符合实际情况，因此，通过拟合曲线对Ma=1.2的平飞需用推力进行修正，结果如下：

	0.4	0.6	0.8	1.2	1.4	1.6	1.8	2.0	2.2
0	8.27	2.85	1.56	176.12	217.45	2.03	308.96	355.93	404.78
5	16.19	6.75	3.78	105.17	119.38	143.67	166.87	191.63	217.46
10	26.50	14.67	8.24	57.16	66.13	76.63	86.33	97.84	110.02
13	43.76	23.76	13.36	42.55	51.45	56.02	59.92	66.330	73.34
16	56.36	38.30	21.54	40.84	48.73	48.08	46.99	49.73	53.15
20	66.236	42.113	30.562	56.147	63.460	54.733	46.276	44.978	44.69

计算最大推力和全加力推力，得到的不同高度下推力曲线图如图4-1所示，（PQJL为全加力状态，PZD为最大加力状态，PPX为平飞需用推力）：

(a) H=0Km

(b) H=5m

(c) H=10Km

(d) H=13Km

(e) H=16Km

(f) H=20Km

图4-1 不同高度下的推力曲线

由推力曲线，可以得到飞机由推力特征的飞行包线如图4-2：

图4-2 飞机的飞行包线

4.2 定常上升性能

飞机在定常直线爬升飞行时（忽略迎角及发动机安装角）

其中

            上升率

根据以上方程可以得出飞机在不同飞行状态（H，M）下的上升率

在不同高度下绘制曲线，则得到了上升率曲线图（资料图20）。

飞机在同一高度下的最大上升率为

上式中，由于和均为马赫数的函数，所以求解比较麻烦，通常利用曲线图直接读出某一飞行高度下的最大上升率，对应的速度则为此飞行高度下的快升速度。

按照同样的方法可以确定每个飞行高度下的最大爬升角：

最大爬升角对应的速度为最陡上升速度。

根据上面的计算结果，可以绘制曲线（资料图21），图中每条曲线与H轴的交点对应于的情况，这一点的高度，刚好是飞机能完成定直平的最大高度，这就是飞机静升限，对应于最大推力状态和最大加力状态下存在两个不同的静升限。对应于上升率为5m/s的高度则为实用升限。

估算过程：

定常直线爬升飞行中，不同飞行状态（H，M）下的上升率，需要通过全加力状态和最大加力状态两种情况下的计算得到不同高度下上升率关于Ma的曲线，即不同高度的上升率曲线见图4-3。

(a) H=0Km

(b) H=5Km

(c) H=10Km

(d) H=13Km

(e) H=16

图4-3不同飞行高度下的图

由上图可知飞机在不同高度下，由不同发动机状态决定的最大上升率为：

H	最大加力（m/s）	全加力（m/s）
0	99.7	210.6
5	62.6	131.7
10	28.5	160.4
13	9.36	92.1
16	/	31.8

此时对应的快升速度为：

H	最大加力（m/s）	全加力（m/s）
0	265.4262	268.9992
5	247.1286	268.2836
10	250.7079	357.0761
13	250.8095	592.9432
16	/	592.7956

每个高度下最大爬升角为：

H	最大加力（rad）	全加力（rad）
0	0.3790	0.60
5	0.2536	0.5325
10	0.1228	0.3159
13	0.0448	0.1701
16	/	0.0542

此时对应的最陡上升速度为：

H	最大加力（rad）	全加力（rad）
0	265.4262	259.3010
5	238.1538	253.2187
10	253.4929	267.4178
13	255.2355	268.9563
16	/	575.3865

根据以上数据开始绘制曲线：

图4-4曲线

从画出的曲线，可以找到该飞机由全加力状态和最大加力状态约束的的静升限分别为：14.8Km和18Km。寻找两种状态的实用升限分别为：13.7Km和17.5Km。

4.3 爬升方式

根据上一节的结果，如果在每一飞行高度下都达到最大爬升率，则飞机能达到最佳的爬升效果（最小时间）。但是，上述结论是在定常爬升条件下得出的，由于每个高度下的快升速度不同，显然，在爬升的过程中要不断按照快升速度调整油门，这对飞行来说较为困难。因此，通常在飞机爬升过程中要保持某一参数不变，例如等表速爬升、等真速爬升、等马赫数爬升等。一般在中低空（H<11km），按照亚音速等表速爬升，而在11km以上，则按照超音速等马赫数爬升，这种爬升方式既便于操作，也能够获得较短的爬升时间，因此常被采用。

爬升方式的计算，需要涉及到求解曲线与坐标轴围成图形的面积，本文采取的方法是面积积分法，通过将计算的点拟合成曲线，在计算曲线与X轴围成的面积，可以得到三个状态下的飞行时间。

4.3.1 亚音速等表速爬升

飞机的表速可以按照下式计算：

其中

            表速

            真速

            海平面大气密度

            本地大气密度

飞机做定常上升时，上升率为：

加速爬升时上升率小于定常上升率，可表示为定常上升率减去上升率偏差的形式：

其中上升率偏差为：

亚音速等表速爬升时：

其中：

综上：

通过以上方法确定的关系后，可应通过数值积分方法得出爬升时间，在没有计算机的情况下，通常也采用图解积分法。首先绘制曲线，在爬升高度的范围内曲线与横轴所围成的面积就是爬升时间。

估算过程：

从上述公式可知，首先需要计算发动机推力与平飞需用推力之差，进而求解得到上升率；接下来计算表速，可以得到上升率偏差。最后绘制曲线，并通过面积积分法求解面积。

针对最大加力状态，得到的曲线为图4-5。

图4-5 最大加力状态下等表速爬升的曲线

针对全加力状态，得到的曲线为图4-6。

图4-6 全加力状态下等表速爬升的曲线

通过对上图进行面积积分，得到最大加力状态下用时为t = 224.16s，全加力状态下用时为t=134.22s。

4.3.2 超音速等马赫数爬升

当H>11km时，a=const，又有Ma=const，

爬升时间的计算方法同上节。

估算过程：

在最大推力状态，根据飞行包线可知，飞机无法爬升到15Km，并且速度也无法达到1.2Ma，因此仅仅能够实现等表速爬升。

对于全加力状态，超音速等马赫数爬升同样采用面积积分法。得到的曲线为图4-7。

图4-7 全加力状态下等马赫数爬升的曲线

通过对上图进行面积积分，得到全加力状态下用时为t=179.94s。

4.3.3 平飞加速段的求解方法

飞机在11km以下采用亚音速等表速爬升，在11km以上采用等马赫数爬升，两者之间存在速度差距，因此在11km高度需要进行平飞加速，从亚音速加速到超音速。这个加速过程可以按照以下方法求解出加速时间。

根据飞机能量高度的定义，有如下关系式：

因此可以得到：

化简得：

由不同飞行马赫数下的，可以绘制曲线，采用图解积分法就可以计算出平飞加速时间了。

估算过程：

平飞加速段，绘制曲线如图4-8。

图4-8 全加力状态下等马赫数爬升的曲线

通过对上图进行面积积分，得到全加力状态下，平飞加速段用时为t=192.87s。

4.3.4 总用时

最大推力状态，仅能完成等表速爬升，时间为：t = 224.16s。

全加力状态，将三种状态用时相加，有：t = t1 + t2 + t3 = 134.22 + 179.94 + 192.87 = 507.03s，为整个飞行所需时间。

第五章 自主编写的Matlab代码

5.1 RBF径向基函数插值方法实现

RBF径向基函数插值，通过若干RBF径向基函数对自变量进行变换，再将其结果进行线性叠加得到最终的输出。插值过程如下：

1.在一条曲线上选择若干插值点，视情况而定；

2.将各个插值点作为径向基函数的基点，并将各点对应的横坐标作为这个基点的中心；

3.计算每个插值点相对基点的欧氏距离，并进行径向基函数变换；

4.通过各点纵坐标值组成的向量与各个基点变换后的值相除，即可得到对应的权值向量。此时完成整个插值。

5.对于曲线上任意点，给定横坐标后，计算其相对与各个基点的欧式距离并进行基函数变换，在将各个基点处变换后的输出乘以插值中确定的权值，就可以得到该任意点对应的函数值。

Matlab代码：

clear all                            %清理空间变量

load stat.dat                        %数据文件，保存了插值点的横纵坐标

load input.dat                    %待求的未知点输出

                                 

A = stat(:,1);                        %保存横纵坐标

B = stat(:,2);

B = B';

 

e = dist(A,A');                    %求取中心点之间的欧式距离

[m,~] = size(e);

sp = (sum(sum(e))/(m*(m-1)));    %选取径向基函数的作用半径

F = invmulq(e,sp);

W = B/F;

 

k = 0;                            %做出函数图

for i = 0.7:0.005:1

    ex = sqrt((i-A).^2);

    F1 = invmulq(ex,sp);

    k = k+1;

    Y(1,k) = i;

    Y(2,k) = W*F1;

end

plot(Y(1,:),Y(2,:),'b'); 

legend('曲线拟合结果')

x = input;                        %对未知输入量求解函数值

o = length(x);

for i = 1:o

    ex = x(i)-A;

    F1 = invmulq(ex,sp);

    Z(i) = W*F1;

end

output(:,1) = x;

output(:,2) = Z;

disp('输入的未知量')

disp(x')

disp('未知点预测结果')

disp(Z)

文中采用的三种基函数：

1.高斯基函数：

function F = gauss(e,sp)

  F = radbas(-e.^2/(2*sp^2));

end

2.逆多二次基函数：

function F = invmulq(e,sp)

    [m,n]=size(e);

    A = ones(m,n)+(e/sp).^2;

    F = 1./A;

end

3.线性基函数：

function F = lin(e,sp)

  F = abs(e)/sp;

end

5.2 气动计算及性能计算

由于时间仓促，计算程序难免冗余，涉及诸多插值结果及变量定义，在此不一并写出。

第六章 心得体会

从3月12号决定开始完成课程设计到今天3月19号，整整一周的时间，最终终于将报告整理完成。这个过程中自己有所收获，也有一定的思考。

开始拿到课程设计，觉得里面内容繁杂，有无从下手之感。然而真正开始做的时候，才发现资料中已经提供了足够的数据，有的只剩下自己进行迭代求解或者曲线拟合的过程，这一定程度降低了难度，但也意味着更多的计算过程。

看到那么多图表，我的第一反应就是通过插值方法做一定可以提高效率，否则太麻烦。由于自己研究过RBF神经网络，对RBF径向基函数有一定的认识，加之课题组也有RBF径向基函数插值的研究基础，因而我觉得自己通过Matlab编程以实现插值程序。一个晚上的时间，我将程序完成，并且和课题组以前通过Fortran语言编写的RBF插值结果进行了对比，取得了不错精度，因此将其用在了自己的课程设计中。

在计算过程中，为了防止手动计算的计算量，我将所有计算过程都交给了电脑。这个过程是整个课程设计最耗时的一段时间，因为在计算中需要涉及多种变量的定义和批量计算。而对于需要查表的量，只能通过曲线拟合首先得到输出值，再将这个值以文本的形式导入Matlab中，以保证程序的简洁性。因此，不断的重复插值、带回、编程的过程，也让我感觉有点很忙乱。不过等实现了之后，还是会有一种成就感。

课程设计的结束，本科生活仅剩下毕业设计一个部分了。能够在大四上学期进入课题组开始自己的研究工作，并且将之前所用到课程设计中，使我感到自己做的事情都是有意义的。而通过课程设计，进一步熟悉了Matlab，Tecplot，Word的使用，也为自己未来的硕士生活打下坚实的基础。

本课程设计的参考资料主要是《飞机气动及飞行性能计算》的课程设计指导资料和学校的《飞行器性能计算》讲义，同时也参考了一些课题组相关的学术论文以完成程序。