[转]南航高人编写飞机设计迭代MATLAB程序深度剖析

由于飞机设计是一个需要多次反复迭代的重复过程，需要进行大量的计算，因此我们自己编写了一些MATLAB程序。下面是程序源代码

[1]全机总体参数计算程序
%本程序计算全机各部分参数
%符号说明
%A_ht：平尾尾容量；Angle_to：起飞迎角（deg）；A_vt：垂尾尾容量
%AR_ht：平尾展弦比；AR_vt：垂尾展弦比；AR_w：机翼展弦比
%Cl_max：全机最大升力系数（失速时）；Cl_min：全机最小升力系数（最大速度时）
%Cl_cru：全机巡航升力系数（即设计升力系数）
%H_fus：机身中段高度（m）
%K_to：起飞升阻比；K_V_max：最大平飞速度时的升阻比
%L_fus：机身长度（m）；L_ht：平尾尾力臂（m）；L_vt：垂尾尾力臂（m）
%P_to：起飞时螺旋桨输出功率（w）；P_mot_to：起飞时电机输出功率（w）
%P_max：螺旋桨最大输出功率（w）；P_mot_max：电机最大输出功率（w）
%PWR_to：起飞功重比（W/kg）；PWR_max：最大功率重量比（W/kg）
%Re_max：全机最大雷诺数；Re_cru：全机巡航雷诺数；Re_min：全机最小雷诺数
%S_w：机翼面积(m^2)；S_ail：副翼面积（m^2）；S_ht：平尾面积(m^2)
%S_ele：升降舵面积(m^2)；S_vt：垂尾面积(m^2)；S_rud：方向舵面积(m^2)
%TR_w：机翼根梢比；TR_ht：平尾根梢比；TR_vt：垂尾根梢比
%T_to：起飞时螺旋桨拉力（kg）；T_max：最大速度时螺旋桨拉力（kg）
%V_max：最大飞行速度（m/s）；V_min：最小飞行速度（m/s）；V_cru：巡航飞行速度（m/s）
%WL：翼载荷(kg/m^2)；Wto：起飞重量(g)；W_fus：机身宽度（m）
%W_w：机翼结构重量（g）；W_fus：机身结构重量（g）；W_ht：平尾结构重量（g）
%W_vt：垂尾结构重量（g）；W_lg：起落架结构重量（g）；Wstr：结构重量（g）
%Wfix：固定设备重量（g）；Wele：电子设备（舵机、接收机、连杆等）重量（g）
%Weff：有效载荷（g）；Wthr：推进系统重量（g）；Wcel：电池重量（g）
%Wmot：电机、减速器、调速器、安装架重量和（g）；Wbld：螺旋桨和桨罩重量（g）
%X0_vt：垂尾前缘后掠角（deg）；X4_vt：垂尾1/4弦线后掠角（deg）
%c_w：机翼平均气动弦长(m)；c_ail：副翼弦长（m）；c_ht：平尾平均气动弦长(m)
%c_ele：升降舵弦长(m)；ca_vt：垂尾平均气动弦长(m)；c_rud：方向舵弦长(m)
%den:大气密度（kg/m^3）;dih_i:机翼中央翼段上反角（deg）;dih_o:机翼外翼段上反角(deg)
%eta_bld：螺旋桨效率；eta_gear：减速器效率；eta_mot：电机效率
%g：重力加速度（m/s^2）
%l_w_i_eff：机翼中央翼段有效展长（m）；l_w_i：机翼中央翼段实际长度（m）
%l_w_o_eff：机翼外翼段有效展长（m）；；l_w_o：机翼外翼段实际长度（m）
%l_w：机翼展长(m)；l_ail：副翼展长（m）；l_ht：平尾展长(m)
%l_ele：升降舵展长(m)；l_vt：垂尾展长(m)；l_rud：方向舵展长(m)
%mu：空气粘性系数（kg/(m*s)）
%r1：副翼面积与机翼面积之比；r2：副翼弦长与机翼平均气动弦长之比
%r3：升降舵面积与平尾面积之比 
%tc：机翼相对厚度(%)
%全机总体参数及机翼
clc;
clear;
disp('*************全机主要参数及机翼参数计算**************');
disp(' ');
disp('************本程序只支持无扭转的矩形机翼*************');
disp(' ');
Wto=input('请输入初设全机重量值（g） Wto=');
WL=input('请输入初设翼载荷值（kg/m2） WL=');
AR_w=input('请输入机翼初设展弦比 AR_w=');
l_w_i_eff=input('请输入机翼中央翼段初设有效展长（m） l_w_i_eff=');
dih_i=input('请输入机翼中央翼段初设上反角（deg） dih_i=');
dih_o=input('请输入机翼外翼段初设上反角（deg） dih_o=');
V_min=input('请输入最小平飞速度（m/s） V_min=');
V_cru=input('请输入巡航速度（m/s） V_cru=');
V_max=input('请输入最大平飞速度（m/s） V_max=');
r1=input('请输入副翼面积与机翼面积之比 r1=');
r2=input('请输入副翼弦长与机翼平均气动弦长之比 r2=');
disp(' ');
disp('****************************************************');
%Wto=2400;WL=4.17;AR_w=8;dih_i=0;dih_o=3;
den=1.225;mu=1.7894*10^(-5);g=9.81;
%V_min=8;V_cru=15;V_max=20;
S_w=Wto/1000/WL;
l_w=(S_w*AR_w)^0.5;
%l_w_i=1.00;
l_w_i=l_w_i_eff/cos(dih_i*pi/180);
l_w_o_eff=(l_w-l_w_i_eff)/2;
l_w_o=l_w_o_eff/cos(dih_o*pi/180);
c_w=l_w/AR_w;
Cl_max=2*Wto*g/(den*V_min^2*S_w)/1000;
Cl_cru=2*Wto*g/(den*V_cru^2*S_w)/1000;
Cl_min=2*Wto*g/(den*V_max^2*S_w)/1000;
Re_max=den*V_max*c_w/mu;
Re_cru=den*V_cru*c_w/mu;
Re_min=den*V_min*c_w/mu;
fprintf('\n翼载荷为 %.3f kg/m^2\n\n',WL);
fprintf('最大升力系数为 %.3f \n巡航升力系数为 %.3f\n最小升力系数为 %.3f\n\n',Cl_max,Cl_cru,Cl_min);
fprintf('最大雷诺数为 %d \n巡航雷诺数为 %d\n最小雷诺数为 %d\n\n',Re_max,Re_cru,Re_min);
fprintf('机翼面积 %.4f m^2\n机翼展长 %.3f m\n机翼平均气动弦长 %.3f m\n\n',S_w,l_w,c_w);
fprintf('中央翼段有效展长 %.3f m\n中央翼段上反角 %.3f deg\n中央翼段实际长度 %.3f m\n\n',l_w_i_eff,dih_i,l_w_i);
fprintf('外翼段有效展长 %.3f m\n外翼段上反角 %.3f deg\n外翼段实际长度 %.3f m\n\n',l_w_o_eff,dih_o,l_w_o);
%r1=0.08;r2=0.19;
S_ail=r1*S_w/2;
c_ail=r2*c_w;
l_ail=S_ail/c_ail;
fprintf('副翼面积 %.4f m^2\n副翼展长 %.3f m\n副翼弦长 %.3f m\n\n',S_ail,l_ail,c_ail);
%动力系统
disp('********************动力装置计算********************');
disp(' ');
eta_bld=input('请输入螺旋桨效率 eta_bld＝');
eta_gear=input('请输入减速器效率 eta_gear＝');
eta_mot=input('请输入电机效率 eta_mot＝');
K_to=input('请输入起飞离地时升阻比 K_to＝');
Angle_to=input('请输入起飞离地时迎角（deg） Angle_to＝');
K_V_max=input('请输入最大速度时升阻比 K_V_max＝');
disp(' ');
disp('****************************************************');
%eta_bld=0.6;eta_gear=0.98;eta_mot=0.7;K_to=9;Angle_to=10;K_V_max=3.8;
Cl_ld=0.85*Cl_max;       %Cl_ld与Cl_max的比值因飞机而异，约在0.8～0.9之间
V_ld=(2*Wto*g/(1000*den*S_w*Cl_ld))^0.5;
T_to=Wto*(cos(Angle_to*pi/180)+K_to*sin(Angle_to*pi/180))/(1000*K_to);
P_to=T_to*g*V_ld;
eta=eta_bld*eta_gear*eta_mot;
P_mot_to=P_to/eta;
PWR_to=P_mot_to/(Wto/1000);
fprintf('\n起飞离地时螺旋桨拉力为 %.3f kg\n起飞离地时螺旋桨输出功率为 %.3f W\n\n',T_to,P_to);
fprintf('起飞离地时电机输出功率为 %.3f W\n起飞离地时功率重量比为 %.3f W/kg\n\n',P_mot_to,PWR_to);
T_max=Wto/1000/K_V_max;
P_max=T_max*g*V_max;
P_mot_max=P_max/eta;
PWR_max=P_mot_max/(Wto/1000);
fprintf('最大速度时螺旋桨拉力为 %.3f kg\n螺旋桨最大输出功率为 %.3f W\n\n',T_max,P_max);
fprintf('电机最大输出功率为 %.3f W\n最大功率重量比为 %.3f W/kg\n\n',P_mot_max,PWR_max);
%平尾
disp('********************平尾参数计算**********************');
disp(' ');
disp('*****************此程序只支持矩形平尾*****************');
disp(' ');
A_ht=input('请输入平尾尾容量 A_ht＝');
L_ht=input('请输入平尾尾力臂（m） L_ht＝');
AR_ht=input('请输入平尾展弦比 AR_ht＝');
r3=input('请输入升降舵面积与平尾面积之比 r3=');
disp(' ');
disp('*****************************************************');
%A_ht=0.39;L_ht=0.741;AR_ht=3.4;r3=0.2;
S_ht=S_w*c_w*A_ht/L_ht;
l_ht=(S_ht*AR_ht)^0.5;
c_ht=l_ht/AR_ht;
S_ele=r3*S_ht;
c_ele=r3*c_ht;
l_ele=S_ele/c_ele;
fprintf('\n平尾面积 %.4f m^2\n平尾展长 %.3f m\n平尾平均气动弦长 %.3f m\n\n',S_ht,l_ht,c_ht);
fprintf('升降舵面积 %.4f m^2\n升降舵展长 %.3f m\n升降舵弦长 %.3f m\n\n',S_ele,l_ele,c_ele);
%垂尾
disp('********************垂尾参数计算*********************');
disp(' ');
A_vt=input('请输入垂尾尾容量 A_vt＝');
L_vt=input('请输入垂尾尾力臂（m） L_vt＝');
AR_vt=input('请输入垂尾展弦比 AR_vt＝');
TR_vt=input('请输入垂尾根梢比 TR_vt＝');
X4_vt=input('请输入垂尾1/4弦线后掠角（deg） X4_vt＝');
r4=input('请输入方向舵面积与垂尾面积之比 r4=');
disp(' ');
disp('*****************************************************');
%A_vt=0.014;L_vt=0.735;AR_vt=0.92;TR_vt=2.55;X4_vt=30;r4=0.2;
S_vt=S_w*l_w*A_vt/L_vt;
l_vt=(S_vt*AR_vt)^0.5;
c0_vt=2*S_vt/(l_vt*(1+1/TR_vt));
X0_vt=atan(tan(X4_vt*pi/180)+(1-1/TR_vt)/(AR_vt*(1+1/TR_vt)))*180/pi;
ca_vt=2*c0_vt*(1+TR_vt+TR_vt^2)/(3*(1+TR_vt));
S_rud=r4*S_vt; 
l_rud=l_vt;
c_rud=S_rud/l_rud;
fprintf('\n垂尾面积 %.4f m^2\n垂尾展长 %.3f m\n',S_vt,l_vt);
fprintf('垂尾根弦长 %.3f m\n垂尾梢弦长 %.3f m\n垂尾平均气动弦长 %.3f m\n',c0_vt,c0_vt/TR_vt,ca_vt);
fprintf('垂尾前缘后掠角 %.3f deg\n1/4弦线后掠角 %.3f deg\n\n',X0_vt,X4_vt);
fprintf('方向舵面积 %.4f m^2\n方向舵展长 %.3f m\n方向舵弦长 %.3f m\n\n',S_rud,l_rud,c_rud);
%全机及各部分重量
disp('****************全机及各部分重量计算******************');
disp(' ');
W_fus=input('请输入中段机身宽度（m） W_fus＝');
H_fus=input('请输入中段机身高度（m） H_fus＝');
L_fus=input('请输入机身长度（m） L_fus＝');
tc=input('请输入机翼翼型最大相对厚度(%) tc＝');
Wele=input('请输入电子设备（舵机、接收机等）重量（g） Wele=');
Weff=input('请输入有效载荷重量（g） Weff=');
Wmot=input('请输入电机、调速器、减速器重量和（g） Wmot=');
Wcel=input('请输入电池组重量（g） Wcel=');
Wbld=input('请输入螺旋桨重量（g） Wbld=');
disp(' ');
disp('*****************************************************');
%H_fus=0.075;L_fus=1.124;W_fus=0.071;tc=9.154;
%机翼
W_w=453.6*0.14676*(S_w*3.2808^2)^0.4852*AR_w^0.7082*tc^(-0.2210);
%机身
W_fus=453.6*0.07092*(W_fus/H_fus)^0.04832*(L_fus*3.2808)^1.6566;
%平尾
W_ht=453.6*0.1570*(S_ht*3.2808^2)^0.1939;
%垂尾
W_vt=453.6*0.1393*(S_vt*3.2808^2)^0.6729;
%起落架
W_lg=0.07*Wto;
%全机结构
Wstr=W_w+W_fus+W_ht+W_vt+W_lg;
%固定设备
%Wele=220;Weff=500;
Wfix=Wele+Weff;
%推进系统
%Wmot=311;Wcel=332;Wbld=64;
Wthr=Wmot+Wcel+Wbld;
%全机重量
Wto_1=Wstr+Wfix+Wthr;
error=(Wto-Wto_1)/Wto;
if abs(error)<=0.005;
    fprintf('\n机翼重量为 %d g\n',W_w);
    fprintf('机身重量为 %d g\n',W_fus);
    fprintf('平尾重量为 %d g\n',W_ht);
    fprintf('垂尾重量为 %d g\n',W_vt);
    fprintf('起落架重量为 %d g\n',W_lg);
    fprintf('结构重量为 %d g\n',Wstr);
    fprintf('电子设备（操纵系统）重量为 %d g\n',Wele);
    fprintf('有效载荷重量为 %d g\n',Weff);
    fprintf('固定设备重量为 %d g\n',Wfix);
    fprintf('电机、调速器、减速器重量和为 %d g\n',Wmot);
    fprintf('电池组重量为 %d g\n',Wcel);
    fprintf('螺旋桨重量为 %d g\n',Wbld);
    fprintf('推进系统重量为 %d g\n',Wthr);
    fprintf('全机重量为 %d g\n\n',Wto);
else if error<-0.005;
    fprintf('初设重量太轻，请重新确定\n');
end
if error>0.005;
    fprintf('初设重量太重，请重新确定\n');
end
end