zlxn.m

飞机续航性能计算

% No.05
% 飞行性能计算－着陆性能
% 确定飞机的着陆性能，包括滑跑距离、下滑距离、着陆距离、所需跑道长度和着陆时间
% 气动修正考虑了收放起落架、襟翼的影响，但未考虑地面效应，
% 安全高度按15m计算，起飞安全速度取1.2Vld，着陆驾驶误差系数取1.5
% 查某型飞机的技术说明书，进场速度取60m/s，主轮接地速度取48m/s，前轮接地速度取40m/s
% 发动机处于慢车工作状态

function zlxn
clc;
clear;
close;
format short g
% format compact

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


% 着陆质量
m = 7100.0;
G = m*9.81;
S = 23.0;
% 干水泥跑道摩擦系数
ff = 0.03;
% 使用刹车时摩擦系数
fb = 0.25;
% 查起降时气动特性曲线(襟翼25°，起落架放下)
CyAlpha = (0.8-0.2)/14.0;

TCxCy = [
         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]';

% 在0高度上查可用最大状态推力曲线，折算为慢车状态推力
Pky = 30000.0*0.85*0.05;

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

% 进场速度、着陆速度
Vjc = 60.0;
Vjd = 48.0;
Vqj = 40.0;

% 求下滑距离

H1 = 0.0;
Theta = 0.0;
Haq = 15.0;
while(H1 < Haq)
    Theta = Theta + 0.001;
    Cy = 2.0*G*cos(Theta)/(D0*Vjc.^2*S);
    Alpha = (Cy-0.2)/CyAlpha*pi/180.0;
    Cx = interp1(TCxCy(1,:), TCxCy(2,:), Cy, 'cube');
    L1 = quad(@Fun_L1, Vjc, Vjd);
    H1 = L1*sin(Theta);
end

t1 = quad(@Fun_T1, Vjc, Vjd);

disp('下滑距离(m)：');
disp(L1);

% 两轮接地时，取接地迎角8°
Cy = 0.2+CyAlpha*8.0;
Cx = interp1(TCxCy(1,:), TCxCy(2,:), Cy, 'cube');
f = ff;

L21 = quad(@Fun_L2, Vjd, Vqj);
t21 = quad(@Fun_T2, Vjd, Vqj);

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

L22 = quad(@Fun_L2, Vqj, 0.0);
t22 = quad(@Fun_T2, Vqj, 0.0);

L2 = L21+L22;
t2 = t21+t22;
disp('滑跑距离(m)：');
disp(L2);

L = L1+L2;
disp('着陆距离(m)：');
disp(L);

disp('着陆时间(s)：');
disp(t1+t2);

% 所需跑道长度
disp('所需跑道长度(m)：');
disp(L*1.5);

% 被积函数

function y = Fun_L1(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

function y = Fun_T1(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 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 f;

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