DRAFT: This module has unpublished changes.
DRAFT: This module has unpublished changes.

DRAFT: This module has unpublished changes.

CP calculations based on Barrowman's Equations done in MATLAB

 

%  barrowman.m

%  Created by Dan Paulsen on 2/16/15.

%  Copyright (c) 2015 Dan Paulsen. All rights reserved.

 

close all;

clear all;

 

%nose cone 

Ln = 408.84; 

TypeFactor = 0.5;

CN_alpha_n = 2.0;

 

%body

D  = 102.21;

 

%fin

CR = 215.9;

LF = 109.009;

S  = 106.68;

CT = 132.979;

XS = 82.9213;

XF = 2490.22;

N  = 3;

 

%Rocket CG

Xcg = 1710;

 

%nose eq

XN = TypeFactor * Ln;

 

%body

R = D/2;

 

%fins eq

Kfb = 1 + (R/(S+R));

 

%LF = (S^2 + (XS +CT/2 - CR/2)^2)^0.5;

CN_alpha_f = (4 * N * (S/D)^2) / (1 + (1 +((2 * LF)/(CR + CT))^2)^0.5);

CN_alpha_fb = Kfb * CN_alpha_f;

xf =XF + (XS * (CR + 2*CT))/(3 * (CR + CT)) + ((1/6)*(CR + CT - ((CR*CT)/(CR + CT))));

 

%final Cp(Xcp)

CN_alpha_T = CN_alpha_n  + CN_alpha_fb;

Xcp =(CN_alpha_n * XN + CN_alpha_fb * xf) / CN_alpha_T;

StaicMargin = (Xcp-Xcg)/D;

 

 

DRAFT: This module has unpublished changes.

Apogee Simulator done in MATLAB.

 

Features.

- uses documented motor thrust curves from thrustcurves.org

- uses 1976 Standard Atmosphere Modeling

- Simultaneous multi-motor testing

- Numerates and plots AGL, Thrust, Acceleration, Drag, Mach, G-force, Velocity with a one milisecond interval.

 

   
   
   
   

 

%  Rocket.m

%  Created by Dan Paulsen on 2/16/15.

%  Copyright (c) 2015 Dan Paulsen. All rights reserved.

 

close all; clear all;

 

g = 9.81;%m/s^2

Interval = 0.001;

LaunchAltitude = 1327.4;%m 

LaunchTemp = 32;%90F 

 

pc = ['r','g','b'];

Gon = 1;

 

wi = 0;

%for w = 4.54:0.1:15.87

wi = wi+1;

%Weight_lb(wi) = w * 2.20462;

 

 

for m = 1:3

 

%m= 3;    

    

%Motors

if (m== 1) MotorFile = 'AeroTech_L1365.eng'; end

if (m== 2) MotorFile = 'AeroTech_L1420.eng'; end

if (m== 3) MotorFile = 'AeroTech_L2200.eng'; end

 

[mf_Time,mf_Thrust,mf_cMass,mf_pMass] = MotorThrustFromFile( MotorFile );

BurnTime = mf_Time(length(mf_Time));%S

 

%rocket parameters

Cd = 0.4; 

%SolidWorkMass = 14.866;

PayLoadWeight = 5.806;%kg

%AirFrameWeight = SolidWorkMass - (mf_cMass + mf_pMass + PayLoadWeight);

AirFrameWeight = 6.441;%kg

%AirFrameWeight = w;%kg

 

%AirFrameWeight = w

 

FinLength = 0.10668;

FinThickness =  0.003175;

FinCount = 3;

FinArea = FinCount * FinThickness * FinLength;

 

rMass = AirFrameWeight+ PayLoadWeight; %kg

Diameter = 4.024 * 25.4 /1000;%m;

TubeArea = (Diameter^2 * pi() /4);

Area = TubeArea + FinArea; 

 

 

%LaunchTowerExitVelocity

LTEV = 100 * 0.3048; 

towerHeightFlag = 0; 

 

Time(m) = 0;

i = 0;

%Powered Ascent   

while(Time(m) < BurnTime) 

  Time(m) = Interval * i;

  Mass(i+1,m) = RocketMass(Time(m), mf_cMass + rMass +mf_pMass,mf_cMass + rMass, BurnTime);

  Thrust(i+1,m) = RocketThrust(Time(m), BurnTime, mf_Time,mf_Thrust);

  

  

  if(i == 0)

     AccZ(i+1,m) = ((Thrust(i+1,m)/Mass(i+1,m)) - g);

     VelZ(i+1,m) = 0;

     Alt(i+1,m) = LaunchAltitude;

     Drag(i+1,m) = 0;

     Mach(i+1,m) = 0;

  else

     AccZ(i+1,m) = ((Thrust(i+1,m)- Drag(i,m))/ Mass(i+1,m)) - g;

     VelZ(i+1,m) = VelZ(i,m) + AccZ(i+1,m) * Interval;

     Alt(i+1,m)  = Alt(i,m) + VelZ(i+1,m) * Interval;

     [rho,a,T,P,nu,H] = stdatmo(Alt(i+1,m),LaunchTemp,'SI',true);

     Drag(i+1,m) = RocketDrag(rho,VelZ(i+1,m), Area, Cd);

     Mach(i+1,m) = VelZ(i+1,m) / a;

  end

  

  if VelZ(i+1,m) < 0 

    VelZ(i+1,m) = 0;

  end 

  

  if ((VelZ(i+1,m)> LTEV) && (towerHeightFlag == 0)) 

    LaunchTowerHeight(m) = Alt(i+1,m) - LaunchAltitude;

    towerHeightFlag = 1;

  end

  

  

  i= i+1;

end

 

 

BurnOut_AGL(wi,m) = (Alt(i-1,m)- LaunchAltitude) * 3.28084

 

 

%Coasting  

while(VelZ(i-1,m) > 0)

    Time(m) = Interval * i;

    Mass(i+1,m) = RocketMass(Time(m), mf_cMass + rMass + mf_pMass,mf_cMass + rMass, BurnTime);

    Thrust(i+1,m) = 0;

    AccZ(i+1,m) = (-Drag(i,m)/ Mass(i+1,m)) - g;

    VelZ(i+1,m) = VelZ(i,m) + AccZ(i+1,m) * Interval;

    Alt(i+1,m)  = Alt(i,m) + VelZ(i+1,m) * Interval;

    [rho,a,T,P,nu,H] = stdatmo(Alt(i+1,m),LaunchTemp,'SI',true);

    Drag(i+1,m) = RocketDrag(rho,VelZ(i+1,m), Area, Cd);

    Mach(i+1,m) = VelZ(i+1,m) / a;

    i= i+1;

end

 

 

 

 

%free fall 

%droug fall  

%Chute deploy

 

MotorFile

Max_AGL_Feet(m) = (Alt(i-1,m)- LaunchAltitude) * 3.28084;

 

MaxG(m) = max(AccZ(:,m))/g;

Max_VelZ_Feet(m) = max(VelZ(:,m)) * 3.28084;

%Launch_Tower_Height_Feet = LaunchTowerHeight * 3.28084

Max_MACH(m)= max(Mach(:,m));

 

 

end

 

 

FS='FontSize';

SimInfo = strcat('Cd:',num2str(Cd),' Dry Weight:',num2str(rMass * 2.20462,4), ' lb.')

 

 

if (Gon == 1)

 

h = figure(1);

hold on

grid on

for m = 1:3

    t =0:Interval:Time(m);

    plot(t,(Alt(1:length(t),m)-LaunchAltitude) * 3.28084,pc(m))

end

xlabel('Time(s)',FS,14)

ylabel('AGL(ft)',FS,14)

legend('\it L1365','\it L1420','\it L2200','Location','southeast')

title('AGL(ft)',FS,18)

axis([0 27 0 12500])

annotation('textbox',[0.42,0.076,0.1,0.1],'string',SimInfo)

hgexport(h, 'AGL.jpg', hgexport('factorystyle'), 'Format', 'jpeg');

 

 

h = figure(2);

hold on

grid on

for m = 1:3

    t =0:Interval:Time(m);

    plot(t,(VelZ(1:length(t),m)) * 3.28084,pc(m))

end

xlabel('Time(s)',FS,14)

ylabel('Velocity(ft/s)',FS,14)

legend('\it L1365','\it L1420','\it L2200','Location','northeast')

title('Velocity(ft/s)',FS,18)

axis([0 27 0 1200])

annotation('textbox',[0.42,0.813,0.1,0.1],'string',SimInfo)

hgexport(h, 'Velocity.jpg', hgexport('factorystyle'), 'Format', 'jpeg');

 

h = figure(3);

hold on

grid on

for m = 1:3

    t =0:Interval:Time(m);

    plot(t,(AccZ(1:length(t),m)/g),pc(m))

end

xlabel('Time(s)',FS,14)

ylabel('G''s',FS,14)

legend('\it L1365','\it L1420','\it L2200','Location','northeast')

title('G''s',FS,18)

axis([0 27 -4 25])

annotation('textbox',[0.42,0.813,0.1,0.1],'string',SimInfo)

hgexport(h, 'Gs.jpg', hgexport('factorystyle'), 'Format', 'jpeg');

 

h = figure(4);

hold on

grid on

for m = 1:3

    t =0:Interval:Time(m);

    plot(t,(Mach(1:length(t),m)),pc(m))

end

xlabel('Time(s)',FS,14)

ylabel('Mach',FS,14)

legend('\it L1365','\it L1420','\it L2200','Location','northeast')

title('Mach',FS,18)

axis([0 27 0 1])

annotation('textbox',[0.42,0.813,0.1,0.1],'string',SimInfo)

hgexport(h, 'Mach.jpg', hgexport('factorystyle'), 'Format', 'jpeg');

 

h = figure(5);

hold on

grid on

for m = 1:3

    t =0:Interval:Time(m);

    plot(t,(Drag(1:length(t),m)),pc(m))

end

xlabel('Time(s)',FS,14)

ylabel('Drag(N)',FS,14)

legend('\it L1365','\it L1420','\it L2200','Location','northeast')

title('Drag',FS,18)

%axis([0 27 0 1])

annotation('textbox',[0.42,0.813,0.1,0.1],'string',SimInfo)

hgexport(h, 'Drag.jpg', hgexport('factorystyle'), 'Format', 'jpeg');

 

 

h = figure(6);

hold on

grid on

for m = 1:3

    t =0:Interval:Time(m);

    plot(t,(AccZ(1:length(t),m)),pc(m))

end

xlabel('Time(s)',FS,14)

ylabel('Acceleration(ft/s^2)',FS,14)

legend('\it L1365','\it L1420','\it L2200','Location','northeast')

title('Acceleration',FS,18)

annotation('textbox',[0.42,0.813,0.1,0.1],'string',SimInfo)

hgexport(h, 'Acceleration.jpg', hgexport('factorystyle'), 'Format', 'jpeg');

 

h = figure(7);

hold on

grid on

 

for m = 1:3

    t =0:Interval:Time(m);

    plot(t,(Thrust(1:length(t),m)),pc(m))

end

xlabel('Time(s)',FS,14)

ylabel('Thrust(N)',FS,14)

legend('\it L1365','\it L1420','\it L2200','Location','northeast')

title('Thrust',FS,18)

axis([0 5 0 3500])

annotation('textbox',[0.42,0.813,0.1,0.1],'string',SimInfo)

hgexport(h, 'Thrust.jpg', hgexport('factorystyle'), 'Format', 'jpeg');

 

 

end

%end

 

%h = figure(8);

%hold on

%grid on

%for m = 1:3

    %Max_AGL_Feet(wi,m)

    %Weight_lb(wi,m)

%    plot(Weight_lb,Max_AGL_Feet(:,m),pc(m))

%end

%xlabel('Weight(lb)',FS,14)

%ylabel('AGL(ft)',FS,14)

%legend('\it L1365','\it L1420','\it L2200','Location','northeast')

%title('Airframe and Recovery Weight Vs. AGL',FS,18)

%axis([  0 12500])

%annotation('textbox',[0.42,0.813,0.1,0.1],'string',strcat('Cd:',num2str(Cd)))

%hgexport(h, 'WeightVsAGL.jpg', hgexport('factorystyle'), 'Format', 'jpeg');

 

 

 

 

 

DRAFT: This module has unpublished changes.

%  MotorThrustFromFile.m MATLAB module for reading trust curve files

%  Created by Dan Paulsen on 2/16/15.

%  Copyright (c) 2015 Dan Paulsen. All rights reserved

 

function [ mf_Time,mf_Thrust,mf_cMass,mf_pMass] = MotorThrustFromFile( MotorFile )

    mData = importdata(MotorFile);

    vars = fieldnames(mData);

    for i = 1:length(vars)

        assignin('base', vars{i}, mData.(vars{i}));

    end

 

    mf_Thrust(1) = 0; 

    mf_Time(1) = 0;

 

    for i = 1:length(mData.data(:,1))

        mf_Thrust(i+1) = mData.data(i,2);

        mf_Time(i+1) = mData.data(i,1);

    end

 

    mf_cMass = 0;

    mf_pMass = 0;

    

    td = textscan(mData.textdata{4,:},'%s %f %f %s %f %f %s');

    mf_cMass = td{1,6} - td{1,5};

    mf_pMass = td{1,5};

end

 

DRAFT: This module has unpublished changes.

%  RocketMass.m MATLAB module for burn fuel mass lose assume linear lose

%  Created by Dan Paulsen on 2/16/15.

%  Copyright (c) 2015 Dan Paulsen. All rights reserved

 

function Mass = RocketMass( Time, WetMass, DryMass, BurnTime )

%UNTITLED5 Summary of this function goes here

%   Detailed explanation goes here

    if(Time <= BurnTime)

        Mass = DryMass + ((WetMass - DryMass) * (1 - (Time /BurnTime)));  

    else

       Mass = DryMass;

    end

end

DRAFT: This module has unpublished changes.

%  RocketDrag.m MATLAB module for drag calculations

%  Created by Dan Paulsen on 2/16/15.

%  Copyright (c) 2015 Dan Paulsen. All rights reserved

 

function Drag = RocketDrag( rho, Velocity, Area, Cd)

    Drag = 0.5 * rho * Velocity^2 * Area * Cd; 

end

 

DRAFT: This module has unpublished changes.

%  RocketThrust.m MATLAB module for linear interpolation of thrust data

%  Created by Dan Paulsen on 2/16/15.

%  Copyright (c) 2015 Dan Paulsen. All rights reserved

 

function Thrust = RocketThrust(Time, BurnTime, mf_Time,mf_Thrust)

%UNTITLED6 Summary of this function goes here

%   Detailed explanation goes here

 

    Thrust = 0;

    if(Time <= BurnTime)

        for i = 2:length(mf_Time) 

            if((Time > mf_Time(i-1)) && (Time < mf_Time(i)))

                dt = mf_Time(i) - mf_Time(i-1);

                ts = (Time - mf_Time(i-1)) / dt;

                dT = mf_Thrust(i) - mf_Thrust(i-1);

                Thrust = mf_Thrust(i-1) + ts*dT; 

            end

            

            if(Time == mf_Time(i))

                Thrust = mf_Thrust(i);

            end

            

        end

    end

end

 

DRAFT: This module has unpublished changes.
Show Comments (0)