探空火箭優化設計

箭頭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