计算目的和内容
计算目的
巩固用简单推力法计算飞机基本飞行性能、续航性能和起飞着陆性能的计算原理、方法和步骤，培养独立分析和解决工程实际问题的能力，熟练掌握 MATLAB 编程在飞行性能计算中的应用，能够对计算结果进行合理分析与验证，形成完整的实验报告。
计算内容
1.2.1 基本飞行性能计算
计算某歼击飞机发动机以最大状态工作时，在 H=0m、1000m、3000m、5000m、8000m、10000m、11000m、12500m、13500m 等 9 个高度上，M=0.20~1.30（间隔 0.05）等 23 个马赫数下的航迹倾角 γ 和上升率 Vv，并绘制各高度上 γ 和 Vv 随马赫数变化的曲线。
计算上述各高度上的最大航迹倾角 γmax 及对应最陡上升马赫数 Mγ、最快上升率 Vv.max 及对应快升马赫数 Mqc，绘制 γmax 和 Vv.max 随高度变化的曲线，确定理论升限 Hmax.a 和实用升限 Hmax.s。
计算各高度上的最大平飞马赫数 Mmax 和最小平飞马赫数 Mmin（含升力限制 Mmin.a 和推力限制 Mmin.T）。
绘制由 Mmin~H、Mmax~H、Mγ~H 和 Mqc~H 组成的飞行包线。
计算飞机从海平面上升到实用升限的最短上升时间 tmin。
1.2.2 续航性能计算
计算该歼击飞机在特定马赫数 M 和发动机转速 n 下巡航飞行的最大航程 Rcr.max 和最久航时 tcr.max，确定对应的远航马赫数 MR.max、远航发动机转速 nR.max 及久航马赫数 Mt.max、久航发动机转速 nt.max。
1.2.3 起飞着陆性能计算
计算飞机起飞地面滑跑段距离 d1 和时间 t1、起飞空中段距离 d2 和时间 t2 以及起飞离地速度 Vlo。
计算飞机着陆空中段距离 d3 和时间 t3、着陆地面滑跑段距离 d4 和时间 t4 以及着陆接地速度 Vtd，进而得到起飞总距离和着陆总距离。
计算原理与方法
发动机可用推力和平飞需用推力
（1）发动机可用推力Ta的计算
当H ≤ 11000 ∼ m时，Ta = T(1 + ΔT̄i)(1 + ΔT̄j)
当H > 11000 ∼ m时，$T_{a} = T_{a.11}\frac{\rho}{\rho_{11}}$
式中，下标11代表11km高度时的相应参数值。
（2）平飞需用推力TR（平飞阻力D）的计算
$$T_{R} = \frac{m_{av}g}{K} = \frac{m_{av}gc_{D}}{c_{L}} = \frac{m_{av}g\left( c_{D0} + Ac_{L}^{2} \right)}{c_{L}} = \frac{\rho M^{2}a^{2}Sc_{D0}}{2} + \frac{2Am_{av}^{2}g^{2}}{\rho M^{2}a^{2}S}$$
（3）剩余推力ΔT的计算
ΔT = Ta − TR
最小平飞速度和最大平飞速度
由图一（a）可知，Ta曲线与TR曲线的左交点对应推力限制的最小平飞速度Mmin , T，右交点对应最大平飞速度Mmax。由图一（b）可知：
若ΔTi < 0，ΔTi + 1 > 0，则第i点和第i + 1点之间的ΔT = 0对应的M数为Mmin .T；
若ΔTi > 0，ΔTi + 1 < 0，则第i点和第i + 1点之间的ΔT = 0对应的M数为Mmax。
（可通过两点线性插值求解）
图1
注：图一（a）中横坐标为马赫数M，纵坐标为推力T，包含可用推力曲线Ta和平飞需用推力曲线TR；图一（b）中横坐标为马赫数M，纵坐标为剩余推力ΔT，标注了ΔTi < 0、ΔTi + 1 > 0、ΔTi > 0、ΔTi + 1 < 0的区间及对应的Mmin .T和Mmax位置。
另外，由升力系数限制的最小平飞速度Mmin .a由气动特性确定：
$$M_{min.a} = \sqrt{\frac{2m_{av}g}{\rho a^{2}Sc_{L.a}}}$$
真正的最小平飞速度Mmin取Mmin .a和Mmin .T（如果存在的话）中大者。
航迹倾角γ和上升率Vv
航迹倾角γ和最大航迹倾角γmax的计算公式为：
$$\gamma = \arcsin\left( \frac{\Delta T}{m_{av}g} \right)$$
$$\gamma_{\max} = \arcsin\left( \frac{\Delta T_{\max}}{m_{av}g} \right)$$
对应γmax的M数为最陡上升M数Mγ。
上升率VV和最大上升率VV, max 的计算公式为：
$$V_{V} = \frac{\Delta TV}{m_{av}g} = \frac{\Delta TMa}{m_{av}g}$$
$$V_{V \cdot \max} = \frac{(\Delta TMa)_{\max}}{m_{av}g}$$
对应VV ⋅ max 的M数为快升M数Mqc。
所以，求γmax和VV ⋅ max 就转化为分别求ΔTmax和Δ(TMa)max。
（初略解可通过求离散点最大值求解；较精确解可通过局部多项式拟合后求极值求解）
最短上升时间
最短上升时间tmin的计算公式为：
$$t_{\min} = \sum_{i = 1}^{n}\left( \frac{\Delta H}{V_{V,\max}} \right)$$
计算tmin时，显然当n越大（即ΔH越小）时计算结果越精确。按前面给的9个高度的VV ⋅ max 和ΔH来决定tmin误差很大，特别是在升限附近误差更大。建议取ΔH ≤ 500m，补充高度上的VV ⋅ max 值可用现有的9个VV.max 值中相应的三点进行插值，也可从Ta和TR开始插值计算出最后结果。
航程和航时
巡航段最大航程Rcr.max
计算公式为：
$$R_{crmax} = \left( \frac{\eta_{11}KMa_{11}}{gc_{f11}} \right)_{\max \cdot \max}\ln\frac{m_{1}}{m_{2}}$$
式中，m1和m2分别为巡航起始和结束时的飞机质量；η11为11km高度上的推力有效系数，近似取为1.0；耗油率cf.11和升阻比K是M和n的函数；K = cL/cD，其中cD = 2Ta.11(M, n)/(ρ11M2a112S)，$c_{L} = \sqrt{\left\lbrack c_{D} - c_{D0}\left( H = 11\sim km,M \right) \right\rbrack/A(M)}$。
求Rcr.max 就转化为求$\left( \frac{\eta_{11}KMa_{11}}{gc_{f.11}} \right)_{\max \cdot \max}$。取一固定的发动机转速值（例如，n = 80%n额定），计算马赫数M = 0.3，0.4，0.5，……1.5时，对应的$\frac{\eta_{11}KMa_{11}}{gc_{f,11}}$值，并通过抛物线插值法从中寻找最大值$\left( \frac{\eta_{11}KMa_{11}}{gc_{f11}} \right)_{\max}$；依次计算发动机转速n = 80%n额定，81%n额定，82%n额定时，对应的$\left( \frac{\eta_{11}KMa_{11}}{gc_{f11}} \right)_{\max}$值，并通过抛物线插值法从中寻找最大值$\left( \frac{\eta_{11}KMa_{11}}{gc_{f.11}} \right)_{\max \cdot \max}$。由此，可得到Rcr.max 以及与之相对应的MR, max 和nRmax。
巡航段最久航时tcr.max
计算公式为：
$$t_{crmax} = \left( \frac{\eta_{11}K}{gc_{f.11}} \right)_{\max \cdot \max}\ln\frac{m_{1}}{m_{2}}$$
求tcr.max 就转化为求$\left( \frac{\eta_{11}K}{gc_{f11}} \right)_{\max \cdot \max}$。计算方法与航程计算方法类似，可得到tcr.max 以及与之相对应的Mtmax和nt.max 。
离地速度和接地速度
离地速度Vlo和接地速度Vtd的计算公式为：
$$V_{lo} = \sqrt{\frac{2mg}{\rho Sc_{L.lo}}}$$
$$V_{td} = k_{1}\sqrt{\frac{2mg}{\rho Sc_{L.td}}}$$
式中，cL.lo为离地升力系数，cL.td为接地升力系数，可分别取为最大允许离地升力系数和接地升力系数；k1为接地速度修正系数，表明Vtd小于升力平衡重力时的速度，可取为0.95。
安全高度处飞行速度
安全高度取为H = 25 ∼ m。
起飞时，飞机上升至安全高度时的速度近似取为VH = 1.3Vlo；
着陆时，飞机下降至安全高度时的速度近似取为VH = 1.2Vtd。
起飞地面滑跑段的距离和时间
起飞地面滑跑段的距离d1和时间t1的计算公式为：
$$d_{1} = \frac{V_{lo}^{2}}{2g\left( \frac{T_{av}}{mg} - f^{'} \right)}$$
$$t_{1} = \frac{V_{lo}}{g\left( \frac{T_{av}}{mg} - f^{'} \right)}$$
式中，Tav为起飞地面滑跑段的发动机平均可用推力，近似取为0.9Ta, 0，其中Ta, 0为M = 0时发动机以最大状态工作时的可用推力；f′ = (f + 1/Klo)/2，其中，f为地面对机轮的摩擦系数，近似取为0.03；Klo为离地瞬间的升阻比，Klo = cL.lo/cD.lo，cL.lo确定后，cD.lo由M = 0.2时的起飞极曲线近似计算得到。
起飞空中段的距离和时间
起飞空中段的距离d2和时间t2的计算公式为：
$$d_{2} = \frac{mg}{\Delta T_{av}}\left( \frac{V_{H}^{2} - V_{lo}^{2}}{2g} + H \right)$$
$$t_{2} = \frac{d_{2}}{V_{av}}$$
式中，Vav = (Vlo + VH)/2；ΔTav = [(Ta.lo − Dlo) + (Ta.H − DH)]/2；Ta.lo和TaH分别取为发动机以最大状态工作时，速度达到Vlo和
着陆空中段的距离和时间
着陆空中段的距离d3和时间t3的计算公式为：
$$d_{3} = K_{av}\left( \frac{V_{H}^{2} - V_{td}^{2}}{2g} + H \right)$$
$$t_{3} = \frac{d_{3}}{V_{av}}$$
式中，Vav = (Vtd + VH)/2；Kav = (Ktd + KH)/2，可由M = 0.2时的着陆极曲线近似计算，cL.H = 0.7cL.td。
着陆地面滑跑段的距离和时间
着陆地面滑跑段的距离d4和时间t4的计算公式为：
$$d_{4} = \frac{V_{td}^{2}}{g\left( f + \frac{1}{K_{td}} \right)}$$
$$t_{4} = \frac{2V_{td}}{g\left( f + \frac{1}{K_{td}} \right)}$$
式中，f为地面对机轮的摩擦系数，包含刹车作用，近似取为0.3；Ktd为接地瞬间的升阻比，Ktd = cL, td/cD, td，由M = 0.2时的着陆极曲线近似确定。
计算过程与结果分析
计算过程
读取标准大气数据、发动机特性数据、气动特性数据和质量特性数据，建立插值函数实现任意高度和马赫数下参数的线性插值求解。基于 MATLAB 编写计算程序，分模块实现基本飞行性能、续航性能和起飞着陆性能的计算，包含数据读取、参数插值、性能计算、结果输出和图形绘制等功能。绘制航迹倾角和上升率随马赫数变化曲线、最大性能参数随高度变化曲线及飞行包线，直观展示计算结果。
最后将计算结果与理论规律对比，验证方法的正确性和结果的合理性。
结果分析
3.2.1 基本飞行性能结果
图 1 各高度上航迹倾角随马赫数变化曲线
图2
各高度上上升率随马赫数变化曲线
图3 最大航迹倾角和最大上升率随高度变化
图4 飞行包线
图1中各高度航迹倾角随马赫数变化曲线中，低高度（0~3000m）航迹倾角整体数值较大，对应马赫数偏低，随着马赫数增大，航迹倾角先增后减。5000~8000m迹倾角整体小于较低高度，对应马赫数较之于低高度有所增大，曲线变化更为平缓， 10000m 以上航迹倾角整体数值偏小，对应马赫数趋于稳定，马赫数接近 1.0 时航迹倾角呈现负向变化。
各高度上升率随马赫数变化呈现同样的分布，低高度上升率整体数值高，曲线上升段斜率大、下降段斜率小，快升对应的马赫数高于同高度最陡上升对应的马赫数，之后上升率整体低于低高度，曲线上升段与下降段斜率相对均衡，过渡区间拓宽，10000m以上上升率整体大幅降低，曲线贴近横轴，马赫数超过 1.0 时上升率呈现负值。
最大航迹倾角和最大上升率随高度变化均呈单调递减趋势，最大航迹倾角随高度增加平稳下降，变化速率逐渐放缓；最大上升率随高度增加呈非线性下降，低高度下降较快，中高度下降速率减缓，高高度下降速率进一步降低，结合最大上升率变化可确定实用升限约 13100m、理论升限约 13500m。
飞行包线中，最大平飞速度随高度变化呈 “先升后降” 的规律，在 8000~10000m 区间达到较高值；最小平飞速度由升力限制和推力限制两条曲线共同构成，5000m 以上推力限制对应的最小平飞速度超过升力限制对应的数值，且随高度增加快速上升；最陡上升速度和快升速度均随高度增加逐渐增大，最终趋于同一马赫数（0.9）。
表 1 基本飞行性能计算结果
| 高度_m | 最大航迹倾角_deg | 最大上升率_m/s | 最陡上升马赫数 | 快升马赫数 | 最大平飞马赫数 | 升力限制最小马赫数 | 推力限制最小马赫数 | 最小平飞马赫数 |
|---|---|---|---|---|---|---|---|---|
| 0 | 26.981662 | 86.19816 | 0.35 | 0.75 | 0.976269 | 0.179101 | - | 0.179101 |
| 1000 | 23.534241 | 79.98368 | 0.35 | 0.75 | 0.98577 | 0.190125 | - | 0.190125 |
| 3000 | 17.697672 | 68.445024 | 0.45 | 0.8 | 1.001542 | 0.215087 | - | 0.215087 |
| 5000 | 13.477478 | 55.744587 | 0.65 | 0.85 | 1.010684 | 0.244848 | 0.215353 | 0.244848 |
| 8000 | 8.433829 | 37.410848 | 0.75 | 0.9 | 1.021062 | 0.301294 | 0.318108 | 0.318108 |
| 10000 | 5.644693 | 25.45332 | 0.85 | 0.9 | 1.021598 | 0.349398 | 0.433359 | 0.433359 |
| 11000 | 4.094234 | 18.63295 | 0.85 | 0.9 | 1.015721 | 0.37748 | 0.515656 | 0.515656 |
| 12500 | 1.826968 | 8.466439 | 0.9 | 0.9 | 0.982962 | 0.424518 | 0.668233 | 0.668233 |
| 13500 | 0.449156 | 2.081789 | 0.9 | 0.9 | 0.925257 | 0.458961 | 0.816161 | 0.816161 |
如表1所示随着高度增加，最大航迹倾角和最大上升率均单调递减，符合大气密度随高度增加而减小的物理规律，发动机推力下降导致上升能力减弱。最陡上升马赫数和快升马赫数随高度增加逐渐增大，在 12500m 以上均稳定在 0.9，表明高空需更高马赫数才能实现最优上升性能。最小平飞马赫数随高度增加而增大，且在 5000m 以上由推力限制主导，低于该高度由升力限制主导。最大平飞马赫数在 8000~10000m 达到峰值（约 1.02），之后随高度增加略有下降。
3.2.2 续航性能结果
表 2 续航性能计算结果
| 指标 | 数值 |
|---|---|
| 最大航程_km | 2279.456332 |
| 远航马赫数 | 0.9 |
| 远航转速_percent | 90 |
| 最久航时_min | 161.95941 |
| 久航马赫数 | 0.6 |
| 久航转速_percent | 80 |
远航马赫数（0.9）高于久航马赫数（0.6），符合续航性能一章中所学内容：远航需兼顾速度和燃油效率，久航则以最低燃油消耗为目标，采用更低马赫数。
远航转速（90%）高于久航转速（80%），因更高转速可提供更大推力，满足远航时的速度需求，而久航需降低转速减少燃油消耗。
最大航程约 2280km，最久航时约 162min，符合歼击机的典型续航性能特征。
3.2.3 起飞着陆性能结果
表 3 起飞着陆性能计算结果
| 指标 | 数值 |
|---|---|
| 离地速度_m/s | 91.056778 |
| 起飞滑跑距离_m | 1015.278545 |
| 起飞滑跑时间_s | 22.299901 |
| 起飞空中段距离_m | 888.935245 |
| 起飞空中段时间_s | 8.489068 |
| 起飞总距离_m | 1904.21379 |
| 接地速度_m/s | 73.73464 |
| 着陆空中段距离_m | 544.036149 |
| 着陆空中段时间_s | 6.707544 |
| 着陆滑跑距离_m | 973.944219 |
| 着陆滑跑时间_s | 26.417549 |
| 着陆总距离_m | 1517.980369 |
离地速度（91.06m/s）高于接地速度（73.73m/s），因着陆时采用更高升力系数和速度修正系数，降低接地速度以缩短着陆距离。
起飞总距离（1904.21m）大于着陆总距离（1517.98m），符合飞机起降性能中：起飞需克服更大阻力加速至离地速度，着陆时可通过刹车和阻力装置缩短滑跑距离。
起飞滑跑时间（22.30s）短于着陆滑跑时间（26.42s），因着陆时需控制减速速率。
总结及感想
本次飞机飞行性能计算实验，通过运用简单推力法，系统完成了飞机基本飞行性能、续航性能和起飞着陆性能的计算，不仅巩固了发动机推力、气动阻力、航迹倾角、上升率等核心参数的计算原理，更熟练掌握了 MATLAB 编程在数据读取、插值运算、图形绘制和结果分析中的应用技巧。对我而言不仅是《无人机飞行动力学》课程理论知识的实践，更是一次将课堂所学与实践经历互相打通的关系。此前参与大创项目的时候有自己设计过一架固定翼无人机，在性能计算和优化的时候遇到了许多问题：如何在有限的机身重量和发动机功率下，平衡续航里程与爬升能力？为何按经验设定的巡航速度，实际测试时续航总是达不到设计预期？这些设计中的困惑，当时只是简单地找网课和现有的算例与模型来解决，并没有很深刻地理解这其中的联系，但通过本次实验的计算与分析，不仅提高了我运用MATLAB的能力，也让我对飞行动力学的认识更清晰了。
通过本次实验，我不仅熟练掌握了飞行性能计算的原理和 MATLAB 编程技巧，更实现了从 “经验主义” 到 “科学计算” 的思维转变。以前参与无人机设计时，更多是依赖前辈的经验参数和反复试错，不仅效率低下，还容易留下性能隐患。如今能够从推力、阻力、升力等核心力学参数出发，通过计算模型精准预判无人机的飞行性能，再结合实际需求进行参数优化，设计过程更加高效、可靠。这次实验也让我深刻认识到，无人驾驶航空器的设计是理论与实践的高度结合，每一个参数的设定都需要坚实的力学理论支撑。未来，我将把本次实验所学的计算方法和思维方式运用到更多无人机设计项目中。
参考文献
[1] 方振平，陈万春，张曙光。航空飞行器飞行动力学 [M]. 北京：北京航空航天大学出版社，2015.
[2] 杨一栋。飞机飞行动力学 [M]. 南京：南京航空航天大学出版社，2012.
[3] 王正平，李为吉。飞行器性能与计算 [M]. 西安：西北工业大学出版社，2010.
[4] 谢础。航空概论 [M]. 北京：北京航空航天大学出版社，2018.
[5] 张志涌. MATLAB 教程：基于 R2018b [M]. 北京：北京航空航天大学出版社，2019.
[6] 中国民用航空局。民用航空器飞行性能要求 [M]. 北京：中国民航出版社，2016.
[7] 李军。无人机飞行控制与导航技术 [M]. 北京：国防工业出版社，2017.
[8] 刘刚，王敏。小型固定翼无人机设计与性能优化 [M]. 西安：西北工业大学出版社，2019..
附录
MATLAB代码
%% 飞机飞行性能计算
clear; clc; close all;
%% 数据
dataFile = '2飞行性能计算原始数据.xlsx';
atm_data_raw = readmatrix(dataFile, 'Sheet', '1-1_标准大气数据');
H_atm = atm_data_raw(:, 1);
rho_atm = atm_data_raw(:, 2);
a_atm = atm_data_raw(:, 3);
T_data_raw = readmatrix(dataFile, 'Sheet', '2-1_最大状态时的推力特性');
M_T = T_data_raw(2:end, 1);
H_T = [0, 1000, 3000, 5000, 8000, 10000, 11000, 12500, 13500];
T_table = T_data_raw(2:end, 2:end);
DTi_data_raw = readmatrix(dataFile, 'Sheet', '2-2_最大状态时由进气道引起的推力损失系数');
M_DTi = DTi_data_raw(2:end, 1);
H_DTi = [0, 1000, 3000, 5000, 8000, 10000, 11000];
DTi_table = DTi_data_raw(2:end, 2:end);
if length(M_DTi) >= 22 && M_DTi(22) == 1.15 && M_DTi(21) == 1.20
M_DTi(22) = 1.25;
end
[M_DTi, unique_idx] = unique(M_DTi, 'stable');
DTi_table = DTi_table(unique_idx, :);
DTj_data_raw = readmatrix(dataFile, 'Sheet', '2-3_最大状态时尾喷口对推力的影响系数');
M_DTj = DTj_data_raw(:, 1);
DTj = DTj_data_raw(:, 2);
T11_data_raw = readmatrix(dataFile, 'Sheet', '2-4_11km 高度上的推力特性');
M_T11 = T11_data_raw(:, 1);
n_levels = [0.80, 0.85, 0.90, 1.00];
T11_table = T11_data_raw(:, 2:end);
cf11_data_raw = readmatrix(dataFile, 'Sheet', '2-5_11km 高度上的耗油率特性');
M_cf11 = cf11_data_raw(:, 1);
cf11_table = cf11_data_raw(:, 2:end);
CD0_data_raw = readmatrix(dataFile, 'Sheet', '3-1_零升阻力系数');
M_CD0 = CD0_data_raw(2:end, 1);
H_CD0 = [0, 1000, 3000, 5000, 8000, 10000, 11000, 12500, 13500];
CD0_table = CD0_data_raw(2:end, 2:end);
A_data_raw = readmatrix(dataFile, 'Sheet', '3-2_升致阻力因子');
M_A = A_data_raw(:, 1);
A_values = A_data_raw(:, 2);
CL_data_raw = readmatrix(dataFile, 'Sheet', '3-3_升力系数');
M_CL = CL_data_raw(:, 1);
CL_a = CL_data_raw(:, 2);
CL_lo = 0.56;
CL_td = 0.52;
polar_data_raw = readmatrix(dataFile, 'Sheet', '3-4_起飞着陆极曲线');
CD_polar = polar_data_raw(:, 1);
CL_takeoff = polar_data_raw(:, 2);
CL_landing = polar_data_raw(:, 3);
mass_data = readmatrix(dataFile, 'Sheet', '4-1_质量特性数据');
m_av = mass_data(2, 2:10);
m1 = mass_data(5, 2);
m2 = mass_data(6, 2);
m_takeoff_data = mass_data(9, 2);
m_landing_data = mass_data(10, 2);
S = mass_data(13, 2);
g = 9.80665;
H_calc = [0, 1000, 3000, 5000, 8000, 10000, 11000, 12500, 13500];
M_calc = 0.20:0.05:1.30;
[H_grid_T, M_grid_T] = meshgrid(H_T, M_T);
F_T = griddedInterpolant(M_grid_T, H_grid_T, T_table, 'linear', 'linear');
[H_grid_DTi, M_grid_DTi] = meshgrid(H_DTi, M_DTi);
F_DTi = griddedInterpolant(M_grid_DTi, H_grid_DTi, DTi_table, 'linear', 'linear');
[H_grid_CD0, M_grid_CD0] = meshgrid(H_CD0, M_CD0);
F_CD0 = griddedInterpolant(M_grid_CD0, H_grid_CD0, CD0_table, 'linear', 'linear');
%% 飞行性能计算
nH = length(H_calc);
nM = length(M_calc);
gamma_matrix = zeros(nH, nM);
Vv_matrix = zeros(nH, nM);
Ta_matrix = zeros(nH, nM);
TR_matrix = zeros(nH, nM);
deltaT_matrix = zeros(nH, nM);
for i = 1:nH
H = H_calc(i);
m_avg = m_av(i);
rho = interp1(H_atm, rho_atm, H, 'linear', 'extrap');
a = interp1(H_atm, a_atm, H, 'linear', 'extrap');
for j = 1:nM
M = M_calc(j);
V = M * a;
if H <= 11000
T_base = F_T(M, H);
DTi = F_DTi(M, min(H, 11000));
DTj_val = interp1(M_DTj, DTj, M, 'linear', 'extrap');
Ta = T_base * (1 + DTi) * (1 + DTj_val);
else
T_base_11 = F_T(M, 11000);
DTi_11 = F_DTi(M, 11000);
DTj_val = interp1(M_DTj, DTj, M, 'linear', 'extrap');
Ta_11 = T_base_11 * (1 + DTi_11) * (1 + DTj_val);
rho_11 = interp1(H_atm, rho_atm, 11000, 'linear');
Ta = Ta_11 * (rho / rho_11);
end
CD0 = F_CD0(M, H);
A = interp1(M_A, A_values, M, 'linear', 'extrap');
q = 0.5 * rho * V^2;
CL = (m_avg * g) / (q * S);
CD = CD0 + A * CL^2;
TR = q * S * CD / g;
deltaT = Ta - TR;
sin_gamma = max(min(deltaT / m_avg, 1), -1);
gamma = asind(sin_gamma);
Vv = deltaT * g * V / (m_avg * g);
Ta_matrix(i,j) = Ta;
TR_matrix(i,j) = TR;
deltaT_matrix(i,j) = deltaT;
gamma_matrix(i,j) = gamma;
Vv_matrix(i,j) = Vv;
end
end
gamma_max = zeros(nH, 1);
M_gamma = zeros(nH, 1);
Vv_max = zeros(nH, 1);
M_qc = zeros(nH, 1);
M_max = zeros(nH, 1);
M_min_T = zeros(nH, 1);
M_min_a = zeros(nH, 1);
M_min = zeros(nH, 1);
for i = 1:nH
H = H_calc(i);
m_avg = m_av(i);
rho = interp1(H_atm, rho_atm, H, 'linear', 'extrap');
a = interp1(H_atm, a_atm, H, 'linear', 'extrap');
[gamma_max(i), idx_gamma] = max(gamma_matrix(i,:));
M_gamma(i) = M_calc(idx_gamma);
[Vv_max(i), idx_Vv] = max(Vv_matrix(i,:));
M_qc(i) = M_calc(idx_Vv);
deltaT_row = deltaT_matrix(i,:);
M_max_found = false;
for j = 1:nM-1
if deltaT_row(j) > 0 && deltaT_row(j+1) < 0
M_max(i) = M_calc(j) + (0 - deltaT_row(j)) * (M_calc(j+1) - M_calc(j)) / (deltaT_row(j+1) - deltaT_row(j));
M_max_found = true;
break;
end
end
if ~M_max_found
M_max(i) = M_calc(end);
end
M_min_T_found = false;
for j = 1:nM-1
if deltaT_row(j) < 0 && deltaT_row(j+1) > 0
M_min_T(i) = M_calc(j) + (0 - deltaT_row(j)) * (M_calc(j+1) - M_calc(j)) / (deltaT_row(j+1) - deltaT_row(j));
M_min_T_found = true;
break;
end
end
if ~M_min_T_found
M_min_T(i) = NaN;
end
CL_allow = interp1(M_CL, CL_a, M_calc, 'linear', 'extrap');
M_min_a(i) = sqrt(2 * m_avg * g / (rho * a^2 * S * max(CL_allow)));
if isnan(M_min_T(i))
M_min(i) = M_min_a(i);
else
M_min(i) = max(M_min_a(i), M_min_T(i));
end
end
H_fine = 0:100:15000;
Vv_max_interp = interp1(H_calc, Vv_max, H_fine, 'pchip');
idx_theory = find(Vv_max_interp <= 0, 1, 'first');
H_max_theory = H_fine(min(idx_theory, length(H_fine)));
idx_practical = find(Vv_max_interp <= 5, 1, 'first');
H_max_practical = H_fine(min(idx_practical, length(H_fine)));
H_climb = 0:500:H_max_practical;
Vv_climb = interp1(H_calc, Vv_max, H_climb, 'pchip');
Vv_climb = max(Vv_climb, 0.1);
t_min = 0;
for k = 1:length(H_climb)-1
dH = H_climb(k+1) - H_climb(k);
Vv_avg = (Vv_climb(k) + Vv_climb(k+1)) / 2;
t_min = t_min + dH / Vv_avg;
end
%% 续航性能计算
H_cruise = 11000;
rho_11 = interp1(H_atm, rho_atm, H_cruise, 'linear');
a_11 = interp1(H_atm, a_atm, H_cruise, 'linear');
eta_11 = 1.0;
[n_grid_T11, M_grid_T11] = meshgrid(n_levels, M_T11);
F_T11 = griddedInterpolant(M_grid_T11, n_grid_T11, T11_table, 'linear', 'linear');
F_cf11 = griddedInterpolant(M_grid_T11, n_grid_T11, cf11_table, 'linear', 'linear');
n_range = 0.80:0.01:1.00;
M_range = 0.3:0.1:1.5;
R_param_max_max = 0; M_R_max = 0; n_R_max = 0;
t_param_max_max = 0; M_t_max = 0; n_t_max = 0;
for n_idx = 1:length(n_range)
n = n_range(n_idx);
R_param_max = 0; t_param_max = 0;
M_R_temp = 0; M_t_temp = 0;
for M_idx = 1:length(M_range)
M = M_range(M_idx);
V = M * a_11;
if M < min(M_T11) || M > max(M_T11), continue; end
T_11 = F_T11(M, n);
cf_11 = F_cf11(M, n);
if T_11 <= 0 || cf_11 <= 0, continue; end
CD0_11 = F_CD0(M, H_cruise);
A_11 = interp1(M_A, A_values, M, 'linear', 'extrap');
CD = 2 * T_11 * g / (rho_11 * V^2 * S);
CL_sq = (CD - CD0_11) / A_11;
if CL_sq < 0, continue; end
CL = sqrt(CL_sq);
K = CL / CD;
R_param = eta_11 * K * M * a_11 / (cf_11 / 3600);
t_param = eta_11 * K / (cf_11 / 3600);
if R_param > R_param_max, R_param_max = R_param; M_R_temp = M; end
if t_param > t_param_max, t_param_max = t_param; M_t_temp = M; end
end
if R_param_max > R_param_max_max
R_param_max_max = R_param_max; M_R_max = M_R_temp; n_R_max = n;
end
if t_param_max > t_param_max_max
t_param_max_max = t_param_max; M_t_max = M_t_temp; n_t_max = n;
end
end
R_cr_max = R_param_max_max * log(m1 / m2) / 1000;
t_cr_max = t_param_max_max * log(m1 / m2) / 3600;
%% 起飞着陆性能计算
rho_0 = rho_atm(1);
a_0 = a_atm(1);
m_takeoff = m_takeoff_data;
m_landing = m_landing_data;
f_takeoff = 0.03;
f_landing = 0.30;
H_safe = 25;
V_lo = sqrt(2 * m_takeoff * g / (rho_0 * S * CL_lo));
M_lo = V_lo / a_0;
V_H_takeoff = 1.3 * V_lo;
valid_takeoff = ~isnan(CL_takeoff) & ~isnan(CD_polar);
CD_lo = interp1(CL_takeoff(valid_takeoff), CD_polar(valid_takeoff), CL_lo, 'linear', 'extrap');
K_lo = CL_lo / CD_lo;
T_a0 = F_T(0.2, 0);
DTi_0 = F_DTi(0.2, 0);
DTj_0 = interp1(M_DTj, DTj, 0.2, 'linear', 'extrap');
T_av_takeoff = 0.9 * T_a0 * (1 + DTi_0) * (1 + DTj_0);
f_prime = (f_takeoff + 1/K_lo) / 2;
thrust_weight_ratio = T_av_takeoff / m_takeoff;
d1 = V_lo^2 / (2 * g * (thrust_weight_ratio - f_prime));
t1 = V_lo / (g * (thrust_weight_ratio - f_prime));
T_a_lo = F_T(M_lo, 0) * (1 + F_DTi(M_lo, 0)) * (1 + interp1(M_DTj, DTj, M_lo, 'linear', 'extrap'));
D_lo = CD_lo * 0.5 * rho_0 * V_lo^2 * S / g;
M_H = V_H_takeoff / a_0;
T_a_H = F_T(min(M_H, max(M_T)), 0) * (1 + F_DTi(min(M_H, max(M_DTi)), 0)) * (1 + interp1(M_DTj, DTj, min(M_H, max(M_DTj)), 'linear', 'extrap'));
CL_H = m_takeoff * g / (0.5 * rho_0 * V_H_takeoff^2 * S);
CD_H = interp1(CL_takeoff(valid_takeoff), CD_polar(valid_takeoff), CL_H, 'linear', 'extrap');
D_H = CD_H * 0.5 * rho_0 * V_H_takeoff^2 * S / g;
deltaT_av = ((T_a_lo - D_lo) + (T_a_H - D_H)) / 2;
V_av_takeoff = (V_lo + V_H_takeoff) / 2;
d2 = m_takeoff / deltaT_av * ((V_H_takeoff^2 - V_lo^2) / (2*g) + H_safe);
t2 = d2 / V_av_takeoff;
V_td = 0.95 * sqrt(2 * m_landing * g / (rho_0 * S * CL_td));
M_td = V_td / a_0;
V_H_landing = 1.2 * V_td;
valid_landing = ~isnan(CL_landing) & ~isnan(CD_polar);
CD_td = interp1(CL_landing(valid_landing), CD_polar(valid_landing), CL_td, 'linear', 'extrap');
K_td = CL_td / CD_td;
CL_H_land = 0.7 * CL_td;
CD_H_land = interp1(CL_landing(valid_landing), CD_polar(valid_landing), CL_H_land, 'linear', 'extrap');
K_H_land = CL_H_land / CD_H_land;
K_av = (K_td + K_H_land) / 2;
V_av_landing = (V_td + V_H_landing) / 2;
d3 = K_av * ((V_H_landing^2 - V_td^2) / (2*g) + H_safe);
t3 = d3 / V_av_landing;
d4 = V_td^2 / (g * (f_landing + 1/K_td));
t4 = 2 * V_td / (g * (f_landing + 1/K_td));
%% 图
set(0, 'DefaultAxesFontName', 'SimHei');
set(0, 'DefaultTextFontName', 'SimHei');
colors = lines(nH);
figure('Position', [100, 100, 900, 600]);
hold on;
for i = 1:nH
plot(M_calc, gamma_matrix(i,:), 'LineWidth', 1.5, 'Color', colors(i,:), 'DisplayName', sprintf('H=%dm', H_calc(i)));
end
hold off;
xlabel('马赫数 M'); ylabel('航迹倾角 \gamma (°)');
title('各高度上航迹倾角随马赫数变化曲线');
legend('Location', 'northeast'); grid on; xlim([0.2, 1.3]);
saveas(gcf, '图1_航迹倾角曲线.png');
figure('Position', [100, 100, 900, 600]);
hold on;
for i = 1:nH
plot(M_calc, Vv_matrix(i,:), 'LineWidth', 1.5, 'Color', colors(i,:), 'DisplayName', sprintf('H=%dm', H_calc(i)));
end
hold off;
xlabel('马赫数 M'); ylabel('上升率 V_v (m/s)');
title('各高度上升率随马赫数变化曲线');
legend('Location', 'northeast'); grid on; xlim([0.2, 1.3]);
saveas(gcf, '图2_上升率曲线.png');
figure('Position', [100, 100, 900, 500]);
subplot(1,2,1);
plot(gamma_max, H_calc/1000, 'b-o', 'LineWidth', 2, 'MarkerSize', 6);
xlabel('最大航迹倾角 \gamma_{max} (°)'); ylabel('高度 H (km)');
title('最大航迹倾角随高度变化'); grid on;
subplot(1,2,2);
plot(Vv_max, H_calc/1000, 'r-o', 'LineWidth', 2, 'MarkerSize', 6);
hold on;
yline(H_max_practical/1000, 'k--', sprintf('实用升限 %.0fm', H_max_practical), 'LineWidth', 1.5);
yline(H_max_theory/1000, 'k:', sprintf('理论升限 %.0fm', H_max_theory), 'LineWidth', 1.5);
hold off;
xlabel('最大上升率 V_{v,max} (m/s)'); ylabel('高度 H (km)');
title('最大上升率随高度变化'); grid on;
sgtitle('最大航迹倾角和最大上升率随高度变化');
saveas(gcf, '图3_最大性能参数曲线.png');
figure('Position', [100, 100, 1000, 700]);
hold on;
H_plot = H_calc / 1000;
plot(M_max, H_plot, 'b-', 'LineWidth', 2.5, 'DisplayName', 'M_{max} (最大平飞速度)');
plot(M_min_a, H_plot, 'r--', 'LineWidth', 2, 'DisplayName', 'M_{min,a} (升力限制)');
valid_idx = ~isnan(M_min_T);
if sum(valid_idx) > 1
plot(M_min_T(valid_idx), H_plot(valid_idx), 'r:', 'LineWidth', 2, 'DisplayName', 'M_{min,T} (推力限制)');
end
plot(M_gamma, H_plot, 'g-', 'LineWidth', 2, 'DisplayName', 'M_\gamma (最陡上升速度)');
plot(M_qc, H_plot, 'm-', 'LineWidth', 2, 'DisplayName', 'M_{qc} (快升速度)');
yline(H_max_practical/1000, 'k--', 'LineWidth', 1.5, 'HandleVisibility', 'off');
text(0.9, H_max_practical/1000 + 0.3, sprintf('实用升限 H_{max,s}=%.0fm', H_max_practical), 'FontSize', 10);
hold off;
xlabel('马赫数 M'); ylabel('高度 H (km)');
title('飞行包线', 'FontSize', 16, 'FontWeight', 'bold');
legend('Location', 'southeast'); grid on;
xlim([0, 1.5]); ylim([0, 15]);
saveas(gcf, '图4_飞行包线.png');
%% 结果
results_basic = table(H_calc', gamma_max, Vv_max, M_gamma, M_qc, M_max, M_min_a, M_min_T, M_min, ...
'VariableNames', {'高度_m', '最大航迹倾角_deg', '最大上升率_m_s', '最陡上升马赫数', ...
'快升马赫数', '最大平飞马赫数', '升力限制最小马赫数', '推力限制最小马赫数', '最小平飞马赫数'});
writetable(results_basic, '飞行性能计算结果.xlsx', 'Sheet', '基本性能');
results_cruise = table([R_cr_max; M_R_max; n_R_max*100; t_cr_max*60; M_t_max; n_t_max*100], ...
'VariableNames', {'数值'}, 'RowNames', {'最大航程_km', '远航马赫数', '远航转速_percent', ...
'最久航时_min', '久航马赫数', '久航转速_percent'});
writetable(results_cruise, '飞行性能计算结果.xlsx', 'Sheet', '续航性能', 'WriteRowNames', true);
results_takeoff_landing = table([V_lo; d1; t1; d2; t2; d1+d2; V_td; d3; t3; d4; t4; d3+d4], ...
'VariableNames', {'数值'}, 'RowNames', {'离地速度_m_s', '起飞滑跑距离_m', '起飞滑跑时间_s', ...
'起飞空中段距离_m', '起飞空中段时间_s', '起飞总距离_m', '接地速度_m_s', '着陆空中段距离_m', ...
'着陆空中段时间_s', '着陆滑跑距离_m', '着陆滑跑时间_s', '着陆总距离_m'});
writetable(results_takeoff_landing, '飞行性能计算结果.xlsx', 'Sheet', '起降性能', 'WriteRowNames', true);