aerodymic.m

本段代码是计算

H=20000;%高度
x=10;%攻角
deta=4.5;%舵偏角
M=5;%马赫数

g0=9.80665;%地球表面重力加速度
R1=287;%气体常数
r=6371210;%地球半径
%以下为计算不同高度处温度和大气压
if H<=11000
    T=288.15-0.0065*H;
    P=exp(log(1.0132e5)+1/(R1*0.0065)*log((288.15-0.0065*H)/288.15));
else 
    if 11000<H<=20000
        T=216.15;
        P=exp(log(22610)-1/(R1*216.65)*(H-11000));
    else
        T=216.15+0.001*(H-20000);
        P=exp(log(5460)-1/(R1*0.001)*log((216.65+0.001*(H-20000))/216.65));
    end
end

g=g0*(r/(r+H))^2;%H高度处重力加速度
p=P/(g*R1*T);%H高度处大气密度
u=1.8247e-6*(T/273.15)^(3/2)*(273.15+110.4)/(T+110.4);%H高度处大气粘性系数
v=u/p;%运动粘性系数
Va=20.04*T^0.5;%H高度处音速

%以下为导弹外形参数测量值
DB=0.533;%弹体柱段直径
D=0.533;%弹体尾部直径
d=0.412;%弹头圆弧处直径
s=1.063;%弹头圆弧处直径位置与头柱段结合处距离
LW_1=1.72;%毛弹翼展长
LW_2=1.39;%毛舵面展长
Lb=4.83;%导弹长度
br1=0.754;%翼根长度
br2=0.330;%舵根长度
bt1=0.306;%翼梢长度
bt2=0.189;%舵梢长度
br_2=DB/2*bt2+br2;%毛舵根长度
xbr1=2.681;%翼根前缘到导弹头部顶点距离
xbr2=4.299;%舵根前缘到导弹头部顶点距离
Ln=1.208;%头部长度
Lc=2.657;%柱段长度
Lt=0.965;%尾部长度
n=4;%翼片数目
kq=0.85;%速度阻滞系数

%以下为导弹参数计算
R=DB/2;%圆柱段半径
SB=pi*DB^2/4;%圆柱段截面积
fn=Ln/D;%头部长细比
fc=Lc/D;%柱段长细比
ft=Lt/D;%尾部长细比
fb=Lb/D;%导弹长细比
LW1=LW_1-D;%翼展
LW2=LW_2-D;%舵展
SW1=(bt1+br1)*LW1/2;%翼面面积
SW2=(bt2+br2)*LW2/2;%舵面面积
SW_2=(bt2+br_2)*LW_2/2;%毛舵面面积
gsb1=br1/bt1;%翼面根梢比
gsb2=br2/bt2;%舵面根梢比
gsb2_=br_2/bt2;%毛舵面根梢比
sgb1=1/gsb1;%翼面梢根比
sgb2=1/gsb2;%舵面梢根比
sgb2_=1/gsb2_;%毛舵面梢根比
zxb1=LW1^2/SW1;%翼面展弦比
zxb2=LW2^2/SW2;%舵面展弦比
zxb2_=2*LW_2/(br_2+bt2);%毛舵面展弦比
bA1=4*SW1*(1-sgb1/(1+sgb1)^2)/(3*LW1);%外露翼平均气动力弦
bA2=4*SW2*(1-sgb2/(1+sgb2)^2)/(3*LW2);%外露舵平均气动力弦
b_1=SW1/LW1;%外露翼平均几何弦
b_2=SW2/LW2;%外露舵平均几何弦
X01=(br1-bt1)/(LW1/2);%翼面前缘后掠角正切
X02=(br2-bt2)/(LW2/2);%舵面前缘后掠角正切
X1_4_1=X01-(1-sgb1)/(zxb1*(1+sgb1));%翼面1/4弦线处后掠角正切
X1_4_2=X02-(1-sgb2)/(zxb2*(1+sgb2));%舵面1/4弦线处后掠角正切
Xt1=X01-2*(1-sgb1)/(zxb1*(1+sgb1));%翼面1/2弦线处后掠角正切
Xt2=X02-2*(1-sgb2)/(zxb2*(1+sgb2));%舵面1/2弦线处后掠角正切
X11=X01-4*(1-sgb1)/(zxb1*(1+sgb1));%翼面后缘后掠角正切
X12=X02-4*(1-sgb2)/(zxb2*(1+sgb2));%舵面后缘后掠角正切
xA1=LW1*(1+2*sgb1)*X01/(6*(1+sgb1));%bA1到外露翼根弦的距离
xA2=LW2*(1+2*sgb2)*X02/(6*(1+sgb2));%bA2到外露舵根弦的距离
cosX01=1/(1+X01^2)^0.5;%翼面前缘后掠角余弦
cosX02=1/(1+X02^2)^0.5;%舵面前缘后掠角余弦
cosX11=1/(1+X11^2)^0.5;%翼面1/2弦线后掠角余弦
cosX12=1/(1+X12^2)^0.5;%舵面面1/2弦线后掠角余弦
cosXt1=1/(1+Xt1^2)^0.5;%翼面最大厚度处后掠角余弦
cosXt2=1/(1+Xt2^2)^0.5;%舵面最大厚度处后掠角余弦


%以下部分为计算翼面和舵面零升阻力系数
%首先计算摩擦阻力系数
t1=0.10;%翼面相对剖面厚度
t2=0.04;%舵面剖面相对厚度
xt1=0.15;%翼面上层流附面层转化为紊流附面层的转捩点的相对位置（估计值）
xt2=0.08;%舵面上层流附面层转化为紊流附面层的转捩点的相对位置（估计值）
Re=M*Lb*Va/v;%雷诺数

%查附图2-2-2 M=0时，平板摩擦阻力系数，进行插值
xt=[0,0.1,0.2,0.3,0.5,0.7,0.9,1.0];%转捩点相对位置变量
RE=[1e6,4e6,1e7,4e7,1e8,2e8,4e8,1e9];%雷诺数变量
Cf=[4.50e-3,4.28e-3,4.00e-3,3.73e-3,3.12e-3,2.45e-3,1.71e-3,1.33e-3;
    3.48e-3,3.30e-3,3.06e-3,2.80e-3,2.23e-3,1.64e-3,1.02e-3,0.68e-3;
    3.00e-3,2.80e-3,2.60e-3,2.36e-3,1.87e-3,1.33e-3,0.74e-3,0.42e-3;   
    2.41e-3,2.25e-3,2.08e-3,1.83e-3,1.43e-3,1.00e-3,0.48e-3,0.21e-3;    
    2.11e-3,1.96e-3,1.80e-3,1.59e-3,1.21e-3,0.81e-3,0.37e-3,0.13e-3;
    1.93e-3,1.78e-3,1.60e-3,1.43e-3,1.08e-3,0.72e-3,0.30e-3,0.10e-3;
    1.76e-3,1.61e-3,1.49e-3,1.30e-3,0.99e-3,0.65e-3,0.26e-3,0.07e-3
    1.60e-3,1.47e-3,1.33e-3,1.20e-3,0.89e-3,0.59e-3,0.23e-3,0.02e-3];%平板摩擦阻力系数
%查附图2-2-4 M数对平板摩擦阻力系数的影响，进行插值
xtt=[0,0.1,0.3,0.5,0.7,0.8,1];%转捩点相对位置变量
M1=[0,1,2,3,4,5,6];%马赫数变量
nM=[0,0,0,0,0,0,0;
    0.918,0.922,0.925,0.930,0.934,0.940,0.975;
    0.755,0.768,0.775,0.790,0.821,0.835,0.970;
    0.620,0.630,0.640,0.665,0.705,0.726,0.940;
    0.510,0.526,0.548,0.565,0.620,0.645,0.905;
    0.421,0.440,0.471,0.486,0.535,0.575,0.860;
    0.350,0.382,0.410,0.428,0.470,0.520,0.800];%计及压缩性影响的修正系数
%查附图2-2-1 翼型相对厚度对Cxm的修正系数，进行插值
xttt=[0,0.20,0.40];%转捩点相对位置变量
t=[0,0.02,0.04,0.06,0.08,0.10,0.12];%相对厚度变量
nt=[1.000,1.000,1.000;1.051,1.038,1.025;1.103,1.079,1.051;1.160,1.122,1.081;1.220,1.172,1.115;1.282,1.226,1.149;1.360,1.285,1.187];%计及剖面相对厚度影响对Cxm的修正系数

%以下为分别插值后翼面和舵面各系数的值
Cf1=interp2(xt,RE,Cf,xt1,Re);
Cf2=interp2(xt,RE,Cf,xt2,Re);
nM1=interp2(xtt,M1,nM,xt1,M);
nM2=interp2(xtt,M1,nM,xt2,M);
nt1=interp2(xttt,t,nt,xt1,t1);
nt2=interp2(xttt,t,nt,xt2,t2);

Cxmw1=2*Cf1*nM1*nt1;%单独翼面的摩擦阻力系数
Cxmw2=2*Cf2*nM2*nt2;%单独舵面的摩擦阻力系数

%以下部分为计算翼面波阻系数
K=1;%任意翼剖面波阻系数与菱形翼剖面波阻系数之比
fai=1;%辅助系数

if M>=1
zuhexishu1=zxb1*(M^2-1)^0.5-zxb1*Xt1;%组合变量
zuhexishu2=zxb2*(M^2-1)^0.5-zxb2*Xt2;%组合变量
zuhexishu3=zxb1*(M^2-1)^0.5;%翼面横轴实际取值
zuhexishu4=zxb2*(M^2-1)^0.5;%舵面横轴实际取值
zuhexishu5=zxb1*Xt1;
zuhexishu6=zxb2*Xt2;
zuhexishu7=zxb1*t1^(1/3);
zuhexishu8=zxb2*t2^(1/3);%作为选择哪条曲线进行一维插值的依据（由于原图非常复杂，图线断续严重，故其他变量作估计，进行一维插值）

%由于翼面和舵面根梢比都接近1，附图2-2-9 菱形翼型三维翼波阻系数 根梢比=1，插值分别得：
zuhexishu3_4=[0,1,2,3,4,5,6,7,20];
zuhexishu9_10_1=[2.72,2.27,1.60,1.12,0.82,0.68,0.59,0.51,0.35];
zuhexishu9=interp1(zuhexishu3_4,zuhexishu9_10_1,zuhexishu3);
zuhexishu10=interp1(zuhexishu3_4,zuhexishu9_10_1,zuhexishu4);

CxBw1=zuhexishu9*zxb1*t1^2*(1+fai*(K-1));%翼面波阻系数
CxBw2=zuhexishu10*zxb2*t2^2*(1+fai*(K-1));%舵面波阻系数
else
    CxBw1=0;
    CxBw2=0;
end

%以下部分为计算钝前缘波阻系数
l1=t1;%翼面钝前缘宽度近似
l2=t2;%舵面钝前缘宽度近似
Mn1=M*cosX01;%翼面垂直来流马赫数
Mn2=M*cosX02;%舵面垂直来流马赫数
M_l=[0,0.5,0.9,1,2,3,4,6,10];%马赫数变量
Cxu=[0,0,0,0.47,0.98,1.05,1.07,1.09,1.10];%超音速钝头平板波阻系数
Cxu1=interp1(M_l,Cxu,Mn1);%翼面超音速钝头平板波阻系数
Cxu2=interp1(M_l,Cxu,Mn2);%舵面超音速钝头平板波阻系数
Suw1=0;%翼面前缘面积近似
Suw2=0;%舵面前缘面积近似

Cxuw1=Cxu1*cosX01^3*Suw1/SW1;%弹翼前缘钝度引起的波阻系数
Cxuw2=Cxu2*cosX02^3*Suw2/SW2;%舵面前缘钝度引起的波阻系数

%以下部分为计算钝后缘波阻系数
%根据附图2-2-14 tb/xita=10^2时紊流附面层二维翼钝后缘底部压强系数，插值
M2=[1,1.2,1.4,1.6,1.8,2,2.4,2.8,3.2,3.6,4,10];%马赫数变量
Cp=[0.800,0.520,0.380,0.319,0.270,0.240,0.180,0.150,0.120,0.090,0.070,0.060];%底部压强系数
Cp11=interp1(M2,Cp,M);
Cp12=interp1(M2,Cp,M);

bcp1=SW1/LW1;%翼面平均几何弦长
bcp2=SW2/LW2;%舵面平均几何弦长

if M>=1
xita1=0.037*v^0.2*(M*Va)^(-0.2)*bcp1^0.8;%超音速时在翼面底部截面处的附面层动量损失厚度
xita2=0.037*v^0.2*(M*Va)^(-0.2)*bcp2^0.8;%超音速时在舵面底部截面处的附面层动量损失厚度
else
    xita1=0.0664*v^0.5*(M*Va)^(-0.5)*bcp1^0.5;%亚音速时在翼面底部截面处的附面层动量损失厚度
    xita2=0.0664*v^0.5*(M*Va)^(-0.5)*bcp2^0.5;%亚音速时在舵面底部截面处的附面层动量损失厚度
end

%根据附图2-2-15 仅是tb/xita函数的二维钝后缘底部压强系数，插值
t_xita=[1,10,20,30,40,50,60,70,80,90,100,200,300,400,500,600,700];%t/xita变量
Cp2=[0.100,0.040,0.023,0.016,0.011,0.008,0.005,0.004,0.003,0.001,0.000,-0.005,-0.008,-0.009,-0.010,-0.011,-0.012];%仅是t_xita函数的二维翼钝后缘底部压强系数
Cp21=interp1(t_xita,Cp2,t1/xita1);
Cp22=interp1(t_xita,Cp2,t2/xita2);
Cpb1=Cp11+Cp21;%翼面总的底部压强系数
Cpb2=Cp12+Cp22;%舵面总的底部压强系数
Sb1=0;%翼面钝后缘面积
Sb2=0;%舵面钝后缘面积
if M>=1
Cxbw1=-Cpb1*cosX11*Sb1/SW1;%翼面钝后缘底部阻力系数
Cxbw2=-Cpb2*cosX12*Sb2/SW2;%舵面钝后缘底部阻力系数
else
    Cxbw1=0;
    Cxbw2=0;
end

Cx0w1=Cxmw1+CxBw1+Cxuw1+Cxbw1;%弹翼零升阻力系数
Cx0w2=Cxmw2+CxBw2+Cxuw2+Cxbw2;%舵面零升阻力系数

%以下部分为计算弹体法向力系数
%该弹体头部为钝锥体，d/D约为0.8，查附图2-1-36 钝锥体法向力系数斜率d/D=0.8,插值
deta_s=[2,10,20,30,40,50];%半锥角变量
Ms=[0,0.6,1,1.2,1.4,2,3,4,7];%马赫数变量
Cy1n1=[1.60,1.42,1.00,0.70,0.80,0.75;
       1.80,1.55,1.10,0.80,0.61,0.46;
       2.60,2.23,1.40,0.80,0.51,0.41;
       2.45,2.19,1.80,0.92,0.47,0.39;
       2.20,1.92,1.59,1.10,0.52,0.31;
       1.75,1.50,1.17,0.91,0.61,0.30;
       1.39,1.20,0.99,0.82,0.70,0.41;
       1.22,1.10,0.95,0.85,0.72,0.51;
       1.15,1.08,0.95,0.85,0.75,0.60];
   
Cy1n=1/57.3*interp2(deta_s,Ms,Cy1n1,atan((D-d)/2/s)*180/pi,M);%弹体头部法向力系数

Cy1t=0;%弹体尾部法向力系数

if M>=1
    CX=0.35;%超音速时绕圆柱体的阻力系数
else
    CX=1.2;%亚音速时绕圆柱体的阻力系数
end

Cy1f=CX*4*(sin(x*pi/180))^2*(Lc+Lt)/R/pi;%弹体的粘性法向力系数

Cy1B=Cy1n+Cy1t+Cy1f;%单独弹体法向力系数

%以下部分为计算弹体零升阻力系数，从而也就得出导弹零升阻力系数以及单独弹体升力系数

%按照求翼面摩擦阻力系数时进行插值
xt3=0;%转捩点相对位置
Cf3=interp2(xt,RE,Cf,xt3,Re);
nM3=interp2(xtt,M1,nM,xt3,M);
Sn=R*(Ln^2+R^2)^0.5;%弹体头部侧面积
Sc=pi*D*Lc;%弹体柱段侧面积
St_=pi*D*Lt;%弹体尾部侧面积
Sce=Sn+Sc+St_;%弹体侧表面面积
ncxb=1.05;%形状修正系数估计值
Cxf=Cf3*nM3*ncxb*Sce/SB;%弹体摩擦阻力系数

%由于钝锥形头部d/D约为0.8，故查附图2-2-30 钝锥体前体阻力系数d/D=0.8,进行插值求头部压差阻力系数
fn1=[2,2.25,2.5,3,3.5,4,5];%头部长细比变量
M3=[0,0.5,1,1.25,1.5,2,3,4,5,7.5,10];%马赫数变量
Cxn1=[-0.005,-0.005,-0.005,-0.005,-0.005,-0.005,-0.005;
      0,0,0,0,0,0,0;
      0.100,0.100,0.100,0.060,0.058,0.050,0.040;
      0.192,0.170,0.142,0.120,0.093,0.079,0.055;
      0.200,0.172,0.147,0.115,0.090,0.073,0.051;
      0.182,0.151,0.129,0.099,0.077,0.061,0.042;
      0.156,0.128,0.108,0.080,0.062,0.050,0.036;
      0.145,0.119,0.100,0.073,0.058,0.046,0.031;
      0.140,0.112,0.095,0.069,0.053,0.042,0.029;
      0.130,0.104,0.087,0.061,0.048,0.038,0.024;
      0.128,0.102,0.085,0.060,0.048,0.037,0.024;];

Cxn=interp2(fn1,M3,Cxn1,fn,M); %弹体头部压差阻力系数

%查附图2-2-39 尾部不收缩时旋成体底部压强系数，进行插值
M4=[0,0.8,0.9,1,1.2,1.4,1.6,1.8,2,3,4,5,6,7,8,9,10];%马赫数变量
Cxb_=[0,0,0.145,0.290,0.273,0.269,0.266,0.231,0.198,0.103,0.067,0.049,0.037,0.030,0.024,0.021,0.020];
Cxb=interp1(M4,Cxb_,M);%尾部不收缩时底部阻力系数

Cx0B=Cxf+Cxn+Cxb;%弹体零升阻力系数


Cx0=Cx0w1*n*SW1/SB+Cx0B+Cx0w2*kq*n*SW2/SB;%全弹零升阻力系数

CyB=Cy1B-Cx0B*x*pi/180;%单独弹体的升力系数


tao1=R/(LW_1/2);%比值1
tao2=R/(LW_2/2);%比值2

%查附图图2-1-43 翼体干扰因子，进行插值
tao=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1];%比值变量
KWB=[1,1.08,1.16,1.25,1.34,1.44,1.54,1.66,1.77,1.89,2];%xx情况有弹体影响下的外露弹翼升力与单独外露弹翼升力之比
KBW=[0,0.12,0.28,0.43,0.60,0.80,1.00,1.21,1.46,1.72,2];%xx情况弹翼对弹体的干扰升力与单独外露弹翼升力之比
kWB=[1,0.97,0.95,0.94,0.935,0.94,0.945,0.95,0.97,0.99,1];%deta0情况有弹体影响下的外露弹翼升力与单独外露弹翼升力之比
kBW=[0,0.10,0.20,0.305,0.41,0.515,0.60,0.70,0.80,0.90,1];%deta0情况弹翼对弹体的干扰升力与单独外露弹翼升力之比
KWB1=interp1(tao,KWB,tao1);
KBW1=interp1(tao,KBW,tao1);
kWB1=interp1(tao,kWB,tao1);
kBW1=interp1(tao,kBW,tao1);

KWB2=interp1(tao,KWB,tao2);
KBW2=interp1(tao,KBW,tao2);
kWB2=interp1(tao,kWB,tao2);
kBW2=interp1(tao,kBW,tao2);

%由于zxb1*Xt1约为1，故查附图2-1-11 单独翼升力系数斜率，插值