aerodymic.m

本段代码是计算

zhxs1=[0,0.25,0.5,1.0,1.5];%zxb*t^(1/3)
zhxs2=[-3,-1.5,-1,0,0.5,1,1.5,2,3,4,5,7,20];%zxb*(M^2-1)^0.5
Cyw_x11=[0.0200,0.0200,0.0200,0.0200,0.0200;
         0.0220,0.0228,0.0232,0.0232,0.0232;
         0.0240,0.0240,0.0255,0.0253,0.0242;
         0.0306,0.0299,0.0296,0.0269,0.0231;
         0.0330,0.0312,0.0300,0.0260,0.0221;
         0.0325,0.0310,0.0289,0.0245,0.0211;
         0.0290,0.0280,0.0262,0.0229,0.0200;
         0.0270,0.0254,0.0240,0.0212,0.0190;
         0.0230,0.0220,0.0190,0.0181,0.0163;
         0.0160,0.0160,0.0160,0.0160,0.0160;
         0.0130,0.0130,0.0130,0.0130,0.0130;
         0.0094,0.0094,0.0094,0.0094,0.0094;
         0.0040,0.0040,0.0040,0.0040,0.0040];
if M>=1
Cyw_x1=interp2(zhxs1,zhxs2,Cyw_x11,zxb1*t1^(1/3),zxb1*(M^2-1)^0.5);
else
    Cyw_x1=interp2(zhxs1,zhxs2,Cyw_x11,zxb1*t1^(1/3),-zxb1*(1-M^2)^0.5);
end

Cyw1=Cyw_x1*zxb1;%翼面升力系数斜率

%舵面为大展弦比，其升力系数斜率计算过程如下
Cyw_x2_1=pi*zxb2/57.3/(1+(1+(1-(M*cosXt2)^2)*(zxb2/(2*cosXt2))^2)^0.5);%亚音速升力系数斜率公式
Cyw_x2_1_=pi*zxb2_/57.3/(1+(1+(1-(M*cosXt2)^2)*(zxb2_/(2*cosXt2))^2)^0.5);%亚音速升力系数斜率公式


M0x=(10+0.91*zxb2^3)/(10+zxb2^3);
M0x_=(10+0.91*zxb2_^3)/(10+zxb2_^3);
Mx=M0x+(1-M0x)*(1-cosXt2)^2;%舵面临界马赫数
Mx_=M0x_+(1-M0x_)*(1-cosXt2)^2;%毛舵面临界马赫数
Cyw_x2_x=pi*zxb2/57.3/(1+(1+(1-(Mx*cosXt2)^2)*(zxb2/(2*cosXt2))^2)^0.5);%舵面临界升力系数斜率
Cyw_x2_x_=pi*zxb2_/57.3/(1+(1+(1-(Mx*cosXt2)^2)*(zxb2_/(2*cosXt2))^2)^0.5);%毛舵面临界升力系数斜率
y=(1+pi*zxb2)/(3+pi*zxb2)*(2-2/3*gsb2^(-0.5)-gsb2^(-2));%翼面中间参数1
y_=(1+pi*zxb2_)/(3+pi*zxb2_)*(2-2/3*gsb2_^(-0.5)-gsb2_^(-2));%毛翼面中间参数1
z=Mx*Cyw_x2_x+zxb2^2/((3*pi*zxb2)*(pi*zxb2/Cyw_x2_x-1)*cosXt2^(2/3)/Cyw_x2_x);%翼面中间参数2
z_=Mx_*Cyw_x2_x_+zxb2_^2/((3*pi*zxb2_)*(pi*zxb2_/Cyw_x2_x_-1)*cosXt2^(2/3)/Cyw_x2_x_);%毛翼面中间参数2
beta_=(M-Mx)*(1+(Mx/M)^y)^2;
beta_1=(M-Mx_)*(1+(Mx_/M)^y_)^2;
Cyw_x2_2=1/57.3/((1/Cyw_x2_x)*(Mx/M)^z+beta_/4);%舵面跨音速升力系数斜率公式
Cyw_x2_2_=1/57.3/((1/Cyw_x2_x_)*(Mx_/M)^z+beta_1/4);%毛舵面跨音速升力系数斜率公式

%超音速时，查附图2-1-3 二维翼型剖面法向力系数斜率，进行插值
t_2=[0.25,0.15,0.03];%舵面相对厚度变量
M_2=[1.2,1.5,2,2.5,5,7.5,10];%马赫数变量
Cyw_x2_3_=[0.250,0.200,0.120;
          0.120,0.074,0.062;
          0.050,0.044,0.040;
          0.036,0.033,0.030;
          0.022,0.018,0.016;
          0.020,0.014,0.011;
          0.020,0.013,0.010];
Cyw_x2_3=interp2(t_2,M_2,Cyw_x2_3_,t2,M);%舵面超音速时升力系数斜率
Cyw_x2_3_=interp2(t_2,M_2,Cyw_x2_3_,t2,M);%毛舵面超音速时升力系数斜率

if M<=Mx
    Cyw_x2=Cyw_x2_1;
    Cyw_x2_=Cyw_x2_1_;
else
    if Mx<M<=1.5
        Cyw_x2=Cyw_x2_2;
        Cyw_x2_=Cyw_x2_2_;
    else
        Cyw_x2=Cyw_x2_3;
        Cyw_x2_=Cyw_x2_3_;
    end
end

Cyw2=Cyw_x2*zxb2;%舵面升力系数斜率（大展弦比）
Cyw2_=Cyw_x2_*zxb2_;%毛舵面升力系数斜率（大展弦比）
CyT_=Cyw2_*x;%舵面升力系数（大展弦比）


e=2;%舵面处下洗角
deta_e=0;%翼面安装角
CyBW=Cyw1*(KBW1*x+kBW1*deta_e)*SW1/SB;%弹翼对弹体的干扰升力系数
CyWB=Cyw1*(KWB1*x+kWB1*deta_e)*SW1/SB;%有弹体时作用在外露弹翼上的升力系数
CyBT=Cyw2*(KBW2*(x-e)+kBW2*deta)*kq*SW2/SB;%尾翼对弹体的干扰升力系数
CyTB=Cyw2*(KWB2*(x-e)+kWB2*deta)*kq*SW2/SB;%有弹体时作用在外露尾翼上的升力系数

CyWB=CyB+CyBW+CyWB;%翼体组合体升力系数

Cy=CyWB+CyTB+CyBT%全弹升力系数
Y=0.5*p*(M*Va)^2*SB*Cy%全弹升力

%以下部分求导弹诱导阻力系数，从而得到导弹阻力系数
CxiBW=CyBW*tan(x*pi/180);%由于弹翼干扰在弹体上产生升力而引起的诱导阻力系数
CxiBT=CyBT*tan(x*pi/180);%由于尾翼干扰在弹体上产生升力而引起的诱导阻力系数

if M>=1
E=1.5/(1+fn);
else
    E=-0.2;
end

Cx1B_fj=E*(Cy1n+Cy1t)*sin(x*pi/180);%附加切向力系数
CxiB=Cy1B*sin(x*pi/180)+Cx1B_fj*cos(x*pi/180);%单独弹体诱导阻力系数

%查附图2-2-59 弹翼前缘形状对吸力的修正系数，进行插值
CY=[0,0.05,0.10,0.15,0.20,0.25,0.30,0.35,0.40,0.45,0.50,0.55,0.60,0.80,1.00,1.50,2.00,2.50,Inf];%升力系数变量
ybsl0=[0.61,0.48,0.38,0.31,0.25,0.21,0.17,0.14,0.12,0.10,0.08,0.07,0.06,0.05,0.04,0.03,0.02,0.01,0.00];%弹翼前缘形状对吸力的修正系数变量
ybsl=interp1(CY,ybsl0,Cy);
zuhecanshu=0;
CF=ybsl*zuhecanshu*(CyWB)^2;


CyW=Cyw1*x;%翼面升力系数
CyT=Cyw2*x;%舵面升力系数
CxiW=0.38*CyW^2/(zxb1-0.8*CyW*(zxb1-1))*(zxb1/cosXt1+4)/(zxb1+4)*SW1/SB;%单独弹翼诱导阻力系数
CxiT=0.38*CyT^2/(zxb2-0.8*CyT*(zxb2-1))*(zxb1/cosXt2+4)/(zxb2+4)*SW2/SB;%单独尾翼诱导阻力系数
CyWB_=CyWB-CyW*SW1/SB;%由于弹体干扰在弹翼上引起的升力增量
CyTB_=CyTB-CyT*SW2/SB;%由于弹体干扰在尾翼上引起的升力增量

if M>=1
CxiWB=CyWB*tan((x+deta_e)*pi/180)-CF;%超音速时有弹体干扰时外露翼的诱导阻力系数
else
CxiWB=CxiW+CyWB_*tan((x+deta_e)*pi/180);%亚音速时有弹体干扰时外露翼的诱导阻力系数
end

%对于大展弦比的舵面，由于sgb2_约为0.5，故查附图2-2-55 弹翼诱导阻力因子1/gsb=0.5
deta_1=1.01;%1+deta
if M>=1
    CxiT_=CyT_*tan(x*pi/180)*SW_2/SB;
else
CxiT_=deta_1/(pi*zxb2_)*CyT_^2*SW_2/SB;
end
Cxi=CxiWB+CxiBW+CxiB+CxiT_;%全弹诱导阻力系数

Cx=Cx0+Cxi%全弹阻力系数
X=0.5*p*(M*Va)^2*SB*Cx%全弹阻力

%以下过程为计算导弹压力中心系数
%查附图2-3-26 钝锥体压力中心系数d/D=0.8，进行插值
deta_ss=[2,10,20,30,40];%半锥角变量
Mss=[0,0.6,1,1.2,2,3,5,8];%马赫数变量
xdn_1=[0.40,0.60,1.10,1.45,2.20;
       0.30,0.55,0.95,1.32,2.00;
       0.37,0.59,0.92,1.23,1.90;
       0.50,0.71,1.01,1.31,1.92;
       0.49,0.68,0.97,1.27,1.70;
       0.48,0.68,0.97,1.27,1.70;
       0.48,0.68,0.97,1.27,1.70;
       0.48,0.68,0.97,1.27,1.70];%钝锥体压力中心系数
xdn_=interp2(deta_ss,Mss,xdn_1,atan((D-d)/2/s)*180/pi,M);%弹体头部压力中心到弹体头部顶点的距离
xdn=xdn_*Ln;%弹体头部压力中心系数

%以下为求导弹弹体粘性法向力压力中心过程
%首先求单独弹体压力中心系数
%由于是锥形头部，查附图2-3-31 锥形头部-圆柱体因子和附图2-3-33 横流干扰因子，分别插图求G1和G2
fn_=[2,4,6,8];%弹体头部长细比变量
fb_=[3,5,10,15,20];%弹体长细比变量
G1_=[1.1,0.5,0.1,0;
     2.3,1.8,1.5,1.0;
     5.2,5.0,4.6,4.1;
     7.7,7.5,7.2,6.9;
     10.5,10.2,10.0,9.8];%锥形头部圆柱体因子
G1=interp2(fn_,fb_,G1_,fn,fb);%锥形头部圆柱体因子实际取值

x_=[0,5,10,15,20,25];%攻角变量
M_x=[0,1,2,3,4,5,7,12];%马赫数变量
G2_=[0.19,0.42,0.70,1.29,2.00,2.05;
     0.20,0.45,0.73,1.31,1.91,1.93;
     0.27,0.55,0.90,1.39,1.81,1.81;
     0.42,0.82,1.19,1.45,1.62,1.62;
     0.70,1.03,1.27,1.38,1.42,1.42;
     0.90,1.10,1.18,1.22,1.24,1.24;
     0.87,1.00,1.01,1.02,1.02,1.02;
     0.90,1.00,1.00,1.00,1.00,1.00];%横流干扰因子
G2=interp2(x_,M_x,G2_,x,M);%横流干扰因子实际取值

Cy1VB=Cy1f;%弹体粘性法向力系数
mzV=-G1*G2*(sin(x*pi/180))^2;
xdV_=-mzV/Cy1VB;%粘性法向力压力中心系数
xdV=xdV_*Lb;%粘性法向力压力中心到导弹头部顶点距离

xdB=(Cy1n*xdn+Cy1VB*xdV)/Cy1B;%单独弹体压力中心到弹体头部顶点距离
xdB_=xdB/Lb;%单独弹体压力中心系数

%然后求单独翼面、舵面压力中心系数
%由于zxb1*X01,zxb2*X02分别大概为2和1，故查附图2-3-9和2-3-8 单独翼压力中心系数，分别插值
gsb=[Inf,5,3,1];%根梢比变量
xdA_11=[0.30,0.28,0.26,0.21;
        0.32,0.30,0.26,0.19;
        0.35,0.32,0.26,0.17;
        0.41,0.38,0.32,0.22;
        0.43,0.41,0.37,0.29;
        0.44,0.43,0.40,0.35;
        0.45,0.44,0.42,0.38;
        0.45,0.44,0.43,0.41;
        0.46,0.46,0.45,0.44;
        0.47,0.47,0.47,0.47;
        0.47,0.47,0.47,0.47;
        0.48,0.48,0.48,0.48;
        0.50,0.50,0.50,0.50];%当攻角为0-5度时xdA_值

    xdA_11_=[0,0,0,0;
             0,0.03,0.06,0.08;
             0,0.04,0.07,0.09;
             0,0.03,0.06,0.08;
             0,0.03,0.05,0.07;
             0,0.02,0.03,0.05;
             0,0.02,0.02,0.03;
             0,0.01,0.01,0.02;
             0,0,0,0.01;
             0,0,0,0;
             0,0,0,0;
             0,0,0,0;
             0,0,0,0];%当攻角为20度时，弹翼压力中心移动量
    xdA_22=[0.35,0.32,0.28,0.22;
            0.38,0.35,0.29,0.20;
            0.41,0.37,0.32,0.20;
            0.46,0.43,0.38,0.27;
            0.48,0.46,0.42,0.33;
            0.49,0.48,0.45,0.38;
            0.49,0.49,0.47,0.42;
            0.49,0.49,0.48,0.45;
            0.49,0.49,0.49,0.48;
            0.49,0.49,0.49,0.49;
            0.48,0.48,0.48,0.48;
            0.48,0.48,0.48,0.48;
            0.47,0.47,0.47,0.47];
    xdA_22_=[0,0,0,0;
         0,0,0.02,0.05;
         0,0.01,0.03,0.06;
         0,0.01,0.03,0.06;
         0,0,0.02,0.05;
         0,0,0.02,0.04;
         0,0,0.01,0.02;
         0,0,0,0;
         0,0,0,0;
         0,0,0,0;
         0,0,0,0;
         0,0,0,0;
         0,0,0,0];
     if M>=1
         f1=zxb1*(M^2-1)^0.5;
         f2=zxb2*(M^2-1)^0.5;
     else
         f1=zxb1*(1-M^2)^0.5;
         f2=zxb2*(1-M^2)^0.5;
     end
 xdA_111=interp2(gsb,zhxs2,xdA_11,gsb1,f1);
 xdA_111_=interp2(gsb,zhxs2,xdA_11_,gsb1,f1);
 xdA_222=interp2(gsb,zhxs2,xdA_22,1000,f2);
 xdA_222_=interp2(gsb,zhxs2,xdA_22_,1000,f2);
     if x>=5
         xdA_1=xdA_111+(x-5)/15*xdA_111_;
         xdA_2=xdA_222+(x-5)/15*xdA_222_;%攻角大于或等于5时公式
     else
        xdA_1=xdA_111;
        xdA_2=xdA_222;%攻角小于5时公式
    end
    
xdW=xbr1+xA1+bA1*xdA_1;%单独翼面压力中心
xdT=xbr2+xA2+bA2*xdA_2;%单独舵面压力中心
xdW_=xdW/Lb;%单独翼面压力中心系数
xdT_=xdT/Lb;%单独舵面压力中心系数

%查附图2-3-34 计算f1之值，和2-3-35 计算f2之值，分别插值，求f1和f2
zhxs6=[0,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1];%D/LW_
f1_lw=[0,0.0325,0.0380,0.0402,0.0401,0.0382,0.0360,0.0325,0.0285,0.0240,0.0200,0.0160,0.0120,0.0085,0.0060,0.0040,0.0022,0.0017,0.0010,0.0005,0];%2*f1/LW
zhxs7=[0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1];%D/LW_
f2_DB=[1.04,1.00,0.96,0.92,0.875,0.83,0.79,0.75,0.71,0.67,0.63];%2*f2/D
f1_lw1=interp1(zhxs6,f1_lw,DB/LW_1);
f1_lw2=interp1(zhxs6,f1_lw,DB/LW_2);
f2_DB1=interp1(zhxs7,f2_DB,DB/LW_1);
f2_DB2=interp1(zhxs7,f2_DB,DB/LW_2);

f11=f1_lw1*LW1/2;%翼面f1
f12=f1_lw2*LW2/2;%舵面f1
f21=f2_DB1*DB/2;%翼面f2
f22=f2_DB2*DB/2;%舵面f2


if M>=1
    xdWB=(xdW+(KWB1-1)*(xdW-f11*X1_4_1))/KWB1;%超音速时有弹体干扰时在外露弹翼的压力中心
    xdTB=(xdW+(KWB2-1)*(xdW-f12*X1_4_2))/KWB2;%超音速时有弹体干扰时在外露尾翼的压力中心
else
    xdWB=(xdW+(KWB1-1)*(xdW-f11*Xt1))/KWB1;%亚音速时有弹体干扰时在外露弹翼的压力中心
    xdTB=(xdW+(KWB2-1)*(xdW-f12*Xt2))/KWB2;%亚音速时有弹体干扰时在外露尾翼的压力中心
end

%以下求由弹翼和尾翼干扰在弹体上引起的升力的压力中心
A1=0.2;
%查附图2-3-36 翼对体的干扰参数，进行插值
D_LW_=[0,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.6,0.7,0.8,0.9];%D/LW_
f_D_LW_=[0,0.095,0.155,0.190,0.220,0.240,0.260,0.280,0.300,0.295,0.330,0.341,0.350,0.356,0.360];%f(D/LW_)
f_D_LW_1=interp1(D_LW_,f_D_LW_,DB/LW_1);%翼面参数
f_D_LW_2=interp1(D_LW_,f_D_LW_,DB/LW_2);%舵面参数

if M>=1
xdBW=xbr1+br1*(xdW_+0.02*zxb1*Xt1)+f21*(M^2-1)^0.5-A1/2*sgb1*LW1*(Xt1+(M^2-1)^0.5);
xdBT=xbr2+br2*(xdT_+0.02*zxb2*Xt2)+f22*(M^2-1)^0.5-A1/2*sgb2*LW2*(Xt2+(M^2-1)^0.5);
else
    xdBW=xbr1+br1*((0.125*zxb1*(1+sgb1)*X01-0.25-LW1/(2*br1)*X1_4_1*f_D_LW_1)*((zxb1*(1-M^2)^0.5-4)/4)^2+(0.25+LW1/(2*br1)*X1_4_1*f_D_LW_1))-A1/2*sgb1*LW1*Xt1;
    xdBT=xbr2+br2*((0.125*zxb2*(1+sgb2)*X02-0.25-LW2/(2*br2)*X1_4_2*f_D_LW_2)*((zxb2*(1-M^2)^0.5-4)/4)^2+(0.25+LW2/(2*br2)*X1_4_2*f_D_LW_2))-A1/2*sgb2*LW2*Xt2;
end
    
xdWB_=xdWB/Lb;%有弹体干扰时在外露弹翼的压力中心系数
xdTB_=xdTB/Lb;%有弹体干扰时在外露尾翼的压力中心系数
xdBW_=xdBW/Lb;%弹翼干扰在弹体上引起的升力的压力中心系数
xdBT_=xdBT/Lb;%尾翼干扰在弹体上引起的升力的压力中心系数

xDWB_=(CyB*xdB_+(CyW*KWB1*xdWB_+CyW*KBW1*xdBW_)*SW1/SB)/(CyB+CyW*(KWB1+KBW1)*SW1/SB);%翼体组合体压力中心系数

xd_=(CyWB*xDWB_+CyBT*xdBT_+CyTB*xdTB_)/Cy%全弹压力中心系数