qfxn.m

飞机起飞性能计算

% No.04
% 飞行性能计算－起飞性能
% 确定飞机的起飞性能，包括滑跑距离、起飞距离、起飞时间、离地速度
% 气动修正考虑了收放起落架、襟翼的影响，但未考虑地面效应
% 安全高度按15m计算，起飞安全速度取1.2Vld

function qfxn
clc;
clear;
close;

global m;
global Pky;
global Cx;
global Cy;
global D0;
global S;
global f;
global Theta;
global Alpha;


% 起飞质量
m = 8150.0;
G = m*9.81;
S = 23.0;
% 干水泥跑道摩擦系数
f = 0.03;
% 查起降时气动特性曲线(襟翼25°，起落架放下)
TCxCy = [
    0.0       0.0551
    0.1014    0.0516
    0.2029    0.0553
    0.3007    0.0691
    0.3986    0.0859
    0.5000    0.1084
    0.6014    0.1390
    0.6993    0.1848
    0.7971    0.2400
    0.8986    0.3097]';
CyAlpha = (0.8-0.2)/14.0;

% 在0高度上查可用推力曲线，由于推力变化在小M下比较平缓，可不考虑速度的影响，取平均值
Pky_full = 50000.0; % 全加力状态
Pky_min = 40000.0;  % 小加力状态
Pky_max = 30000.0;  % 最大状态
TPky = [Pky_max, Pky_min, Pky_full];

% 0高度大气密度、音速
D0 = 1.225;
A0 = 340.294;

% 取前轮离地速度为70m/s
Vql = 70.0;

key = menu('请选择起飞时发动机工作状态', '最大', '小加力' , '全加力');
Pky = TPky(key);

% 求滑跑时间、滑跑距离

% 前轮离地前取滑跑迎角0.33°
Cy = 0.2+CyAlpha*0.33;
Cx = interp1(TCxCy(1,:), TCxCy(2,:), Cy, 'cube');

t11 = quad8(@Fun_T1, 0.0, Vql);
L11 = quad8(@Fun_L1, 0.0, Vql);

% 前轮离地后取滑跑迎角10°
Cy = 0.2+CyAlpha*10.0;
Cx = interp1(TCxCy(1,:), TCxCy(2,:), Cy, 'cube');
% 求离地速度、起飞安全速度
Vld = sqrt(2.0*G/(Cy*D0*S))
Vaq = 1.2*Vld;

t12 = quad8(@Fun_T1, Vql, Vld);
L12 = quad8(@Fun_L1, Vql, Vld);

t1 = t11+t12;
L1 = L11+L12;
disp('滑跑时间(s)：');
disp(t1);
disp('滑跑距离(m)：');
disp(L1);

% 求上升时间、滑跑距离

H2 = 0.0;
Theta = 0.0;
Haq = 15.0;
while(H2 < Haq)
    Theta = Theta + 0.001;
    Cy = 2.0*G*cos(Theta)/(D0*Vld.^2*S);
    Alpha = (Cy-0.2)/CyAlpha*pi/180.0;
    Cx = interp1(TCxCy(1,:), TCxCy(2,:), Cy);
    L2 = quad(@Fun_L2, Vld, Vaq);
    H2 = L2*sin(Theta);
end

t2 = quad(@Fun_T2, Vld, Vaq);

disp('上升时间(s)：');
disp(t2);
disp('上升距离(m)：');
disp(L2);

disp('起飞时间(s)：');
disp(t1+t2);
disp('起飞距离(m)：');
disp(L1+L2);


% 被积函数

function y = Fun_T1(V)
	global m;
	global Pky;
	global Cx;
	global Cy;
	global D0;
	global S;
	global f;
	
	y = m./(Pky-Cx*0.5*D0*V.^2*S-(m*9.81-Cy*0.5*D0*V.^2*S)*f);
return

function y = Fun_L1(V)
	global m;
	global Pky;
	global Cx;
	global Cy;
	global D0;
	global S;
	global f;
	
	y = m*V./(Pky-Cx*0.5*D0*V.^2*S-(m*9.81-Cy*0.5*D0*V.^2*S)*f);
return

function y = Fun_T2(V)
	global m;
	global Pky;
	global Cx;
	global Cy;
	global D0;
	global S;
    global Theta;
    global Alpha;
	
	y = m./(Pky*cos(Alpha)-Cx*0.5*D0*V.^2*S-m*9.81*sin(Theta));
return

function y = Fun_L2(V)
	global m;
	global Pky;
	global Cx;
	global Cy;
	global D0;
	global S;
    global Theta;
    global Alpha;
	
	y = m*V./(Pky*cos(Alpha)-Cx*0.5*D0*V.^2*S-m*9.81*sin(Theta));
return