探空火箭优化设计

箭头60kg，总冲量400KN·s，总推进剂180kg，总质量500kg，请编程设计二级探空火箭并优化设计一二级分配方案，使火箭飞行高度最大。

程序编的比较粗糙，仅供参考。

clear

clc

Cd=0.2;%阻力系数

Cl=0.1;%升力系数

alpha=0; %攻角

S=0.1; %参考面积，m^2

g=9.8; %重力加速度

row=1; %大气密度

Dt=0.001;

n1=1;

for m1=40:20:110  %一级推进剂质量，kg

    m2=180-m1; %二级推进剂质量，kg

    for I1=100000:20000:300000   %一级总比冲

        I2=400000-I1;  %二级总比冲

        for F1=(500*g+100):500:20000

            t1=I1/F1;

            for F2=((500-m1-20)*g+100):500:20000

                t2=I2/F2;

                x=0;

                y=0;

                v=35;

                t=0;

                theta=90*pi/180; %弹道倾角

                flag=0;

                dy=0;

                while 1

                    if  v<=0 %速度为零时是弹道最高点

                        flag=0;

                        break;

                    end

                    if dy<0

                        flag=0; %

                        break;

                    end

                    if t<t1

                        F=F1;

                        M=500-m1/t1*t;

                        if M<= 500-m1-20

                            flag=1;%非正常中止，置flag为1

                            break;

                        end

                    elseif t<t1+t2

                        F=F2;

                        M=500-m1-20-(180-m1)/t2*t;

                        if M<=60 %60kg为弹头质量

                            flag=1;%非正常中止，置flag为1

                            break;

                        end

                    else

                        F=0;

                        M=60;

                    end

                    dv=(F*cos(alpha)-1/2*row*v^2*Cd*S-M*g*sin(theta))/M;

                    dtheta=(F*sin(alpha)+1/2*row*v^2*Cl*S-M*g*cos(theta))/(M*v);

                    dy=v*sin(theta);

                    dx=v*cos(theta);

                    v=v+dv*Dt;

                    theta=theta+dtheta*Dt;

                    y=y+dy*Dt;

                    x=x+dx*Dt;

                    t=t+Dt;

                end

            if flag==1

                continue;

            end

            state(:,n1)=[m1;m2;I1;I2;F1;F2;t1;t2]; %一二级参数配置

            y_store(n1)=y; %高度

            n1=n1+1;

            end

        end

    end

end

[H_max,n]=max(y_store);

disp('最大高度为：')

H_max

disp('------相应参数配置-----')

disp('一级配置：')

disp('一级质量：')

state(1,n)

disp('一级总比冲：')

state(3,n)

disp('一级发动机推力：')

state(5,n)

disp('一级飞行时间：')

state(7,n)

disp('二级配置：')

disp('二级质量：')

state(2,n)

disp('二级总比冲：')

state(4,n)

disp('二级发动机推力：')

state(6,n)

disp('二级飞行时间：')

state(8,n)

运行结果：

最大高度为：

H_max =

  2.0077e+004

------相应参数配置-----

一级配置：

一级质量：

ans =

   100

一级总比冲：

ans =

      240000

一级发动机推力：

ans =

       20000

一级飞行时间：

ans =

    12

二级配置：

二级质量：

ans =

    80

二级总比冲：

ans =

      160000

二级发动机推力：

ans =

  3.8240e+003

二级飞行时间：

ans =

   41.8410