% Vehicle Propellant Budget and Tank Sizing Code

Appendix A.3

A.3.1

tanksize.m

%

% AAE 450 Fall 2001

% Vehicle propellant budget and tank sizing code

% Written by Casey Kirchner for Spring 2001, modified by Jon Edwards for Fall 2001

clear all;

%%%%%%%%%%%%%

% Constants %

%%%%%%%%%%%%%

g_0 = 9.81; % m/s^2 This is a scaling factor to convert between force and mass, NOT necessarily to denote Earth's gravitational influence g_mars = 3.71; % m/s^2 Acceleration of gravity on Mars

Ru = 8314.3;% J/kmol-K, Universal Gas Constant

gamma = 1.66; % Isentropic parameter for helium

Mw_He = 4.002602; % kg/kmol, molecular weight of helium

D_prop = 1.0;% m, propellant tank diameter

rho_g = 1550; % kg/m^3 density of graphite (Table 5.15, Humble)

Ftu_g = 895000000; % Pa, material strength of graphite (Table 5.15, Humble)

fs = 2;% safety factor (Humble p. 269)

rho_al = 2800; % kg/m^3 (alloy 2219, Table 5.15, Humble)

rho_MMH = 0.878*1000; % density of fuel - convert from g/cc (TEP output) to kg/m^3

rho_NTO = 1.444*1000; % density of oxidizer

tw_al = 0.5/1000; % thickness of propellant tank aluminum liner 1/2 mm --> m

% Entry Angle Change Parameters

gamma_deg = 0; % deg, Angle Change

gamma = gamma_deg*(pi/180); % rad, change to radians

Ve = 8390; % m/s Speed at initial entry

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% Vehicle Propulsion Calculations %%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Engine Specifications

Isp_RCS = 306;% R-40A Marquardt RCS Engines

Isp_retro = 313;% Shuttle OMS Engine (Hammond p.215)

landing_mass = 47.324;% tonnes Final mass as vehicle touches down on Mars surface

Ae_At = 55; % Shuttle OMS expansion ratio

O_F = 1.6; % O/F ratio Shuttle OMS Engine (Hammond p.215)

% Delta V required for spin-up maneuver

DV_tether_sul = 99.51; % m/s (Table 4.4.2 Project PERForM Final Report)

% Delta V required for course corrections

DV_RCS = 100*1.02; % m/s (Spring 2001 Data, Table 4.4.2)

% Delta V required for entry angle change

DV_entry = Ve*sin(gamma); % m/s

% Delta V required for perigee lowering/raising in parking orbit

DV_orbit = 0; % m/s

% Apply rocket equation to solve for Hab mass fractions: DV = g_0*Isp*ln(Mfinal/Minitial) - g_local*tb (neglect drag term);

% Parking orbit perigee burn phase

Mfraction_orbit = exp(-DV_orbit/(g_0*Isp_retro));

Minitial_orbit = landing_mass/Mfraction_orbit; % tonnes

prop_orbit_mass = Minitial_orbit - landing_mass; % tonnes

Mfinal_entry = Minitial_orbit; %tonnes

% Entry angle change burn phase

Mfraction_entry = exp(-DV_entry/(g_0*Isp_retro));

Minitial_entry = Mfinal_entry/Mfraction_entry; % tonnes

prop_entry_mass = Minitial_entry - Mfinal_entry; % tonnes

Mfinal_RCS = Minitial_entry; %tonnes

% Enroute maneuvering/attitude phase

Mfraction_RCS = exp(-DV_RCS/(g_0*Isp_RCS));

Minitial_RCS = Mfinal_RCS/Mfraction_RCS;

spun_mass = Mfinal_RCS;

prop_RCS_mass = Minitial_RCS - spun_mass;

Mfinal_sul = Minitial_RCS;

% Mass fractions for Spin-up phase

Mfraction_sul = exp(-DV_tether_sul/(g_0*Isp_RCS));

Minitial_sul = Mfinal_sul/Mfraction_sul;

prop_spin_up_mass = Minitial_sul - Mfinal_sul;

launched_mass = Minitial_sul;

%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Delta V and Prop Masses %

%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Total delta V required for Hab RCS/Retro

HAB_DV = DV_tether_sul + DV_RCS + DV_entry + DV_orbit; % m/s

% Total Hab propellant mass required

prop_total_mass = prop_orbit_mass + prop_entry_mass + prop_RCS_mass + prop_spin_up_mass; % metric tons

%%%%%%%%%%%%%%%%%%

% Engine Sizing %

%%%%%%%%%%%%%%%%%%

% Thrust level

F_HAB = 26689; % N, Thrust of Shuttle OMS engines

num_RCS = 24; % Number of RCS Marquardt R-40A engines

% Engine Physical Dimensions

ME_HAB = 118; % kg (All OMS Engine Data from Astronautix.com

ME_RCS = 10 * num_RCS; % kg (All RCS Engine Data from Astronautix.com

ME_HAB_total = ME_HAB + ME_RCS; % kg

LE_HAB = 1.96; % m

DE_HAB = 1.17; % m

LE_RCS = 1.0; % m

DE_RCS = 0.5; % m

total_engine_weight = ME_HAB_total; % kg

HAB_engine_weight = ME_HAB_total; % kg

%%%%%%%%%%%%%%%

% Tank Sizing %

%%%%%%%%%%%%%%%

% Propellant Masses

ox_total = (O_F/(O_F+1))*(prop_total_mass*1000); % kg (Humble Eq 4.58)

fu_total = (1/(O_F+1))*(prop_total_mass*1000); % kg (Humble Eq 4.58)

% Propellant Volumes

ox_volume = ox_total/rho_NTO*1.03; % m^3 Extra 3% for ullage (Humble p.268)

fu_volume = fu_total/rho_MMH*1.03; % m^3

% Propellant Tank Pressures

V_flow = 10; % m/s typical flow velocity (Humble p.207)

Pc = 862000; % Pa, Chamber pressure of OMS Engine (from

DP_feed = 50000; % Pa, Max pressure drop through feedsystem (Humble Eq 5.18)

DP_cool = 0.15*Pc; % Pa, Pressure loss through cooling system (Humble Eq 5.19)

DP_dynamic_ox = 0.5*rho_NTO*(V_flow)^2; % Oxidizer dynamic pressure drop (Humble Eq 5.16)

DP_dynamic_fu = 0.5*rho_MMH*(V_flow)^2; % Fuel dynamic pressure drop (Humble Eq 5.16)

DP_injector = 0.3*Pc; % Delta P through injector for throttled system (Humble Eq 5.21)

DP_total = DP_feed + DP_dynamic_ox + DP_dynamic_fu + DP_injector + DP_cool; % Pa, Total Pressure loss through feed system

Ptank = Pc + DP_total; % Pa, Tank Pressure needed for OMS Chamber Pressure to be reached (Humble Eq 5.24)

% Pressurant tank pressures

Ti = 300; % K, temp of pressurant tank before blowdown (Humble p.278)

Tf = 200; % K, temp of pressurant tank after blowdown (Humble p.278)

HAB_press_vol = ox_volume + fu_volume; % Initial guess of final volume of pressurant required (Humble p. 279)

HAB_press_mass = Ptank*HAB_press_vol*Mw_He/(Ru*Tf); % kg, mass of the pressurant required (Humble Eq 5.83)

press_on_HAB_vol = HAB_press_mass*Ru*Ti/(Ptank*Mw_He); % m^3, Required pressurant tank volume (Ideal Gas Law)

burst_pressure = fs * Ptank; % Structure's design burst pressure (Humble Eq 5.73)

% Tank Lengths

Loxt_HAB = (ox_volume - (4/3)*pi*(D_prop/2)^3)/(pi*(D_prop/2)^2); % m, Cylindrical length of ox tank (Humble Eq 5.78)

Lfut_HAB = (fu_volume - (4/3)*pi*(D_prop/2)^3)/(pi*(D_prop/2)^2); % m, Cylindrical lenght of fuel tank (Humble Eq 5.78)

Dprt_HAB = 2*(3*press_on_HAB_vol/(4*pi))^(1/3); % m, Diameter of Pressurant Tank (Humble Eq 5.74)

% Tank Surface Areas

Soxt_HAB = pi*D_prop^2 + pi*D_prop*Loxt_HAB; % m^2, Surface area of ox tank (Humble Eq 5.79 and Eq 5.75)

Sfut_HAB = pi*D_prop^2 + pi*D_prop*Lfut_HAB; % m^2, Surface area of fuel tank (Humble Eq 5.79 and Eq 5.75)

Sprt_HAB = pi*Dprt_HAB^2; % m^2, Surface area of pressurant tank (Humble Eq 5.75)

% Tank Wall Thicknesses

tw_oxt_HAB = (D_prop*burst_pressure)/(2*Ftu_g); % m, ox tank thickness (Humble Eq 5.80)

tw_fut_HAB = (D_prop*burst_pressure)/(2*Ftu_g); % m, fuel tank thickness (Humble Eq 5.80)

tw_prt_HAB = (Dprt_HAB*burst_pressure)/(4*Ftu_g); % m, pressurant tank thickness (Humble Eq 5.76)

% Tank Masses

Moxt_HAB = rho_g*Soxt_HAB*tw_oxt_HAB + rho_al*Soxt_HAB*tw_al; % kg, ox tank mass (includes aluminum liner)

Mfut_HAB = rho_g*Sfut_HAB*tw_fut_HAB + rho_al*Sfut_HAB*tw_al; % kg, fuel tank mass (includes aluminum liner)

Mprt_HAB = rho_g*Sprt_HAB*tw_prt_HAB; % kg, pressurant tank mass

total_liner_mass = rho_al*tw_al*(Soxt_HAB + Sfut_HAB); % kg, total mass of tank liners

HAB_tank_mass = Moxt_HAB + Mfut_HAB + Mprt_HAB; % kg, total mass of tanks

%%%%%%%%%%%%%%

% Inert mass %

%%%%%%%%%%%%%%

% According to the AAE 439/539 text, we add 10% of the total inert mass for support structure (p. 281 Humble)

HAB_support_mass = 0.10*(ME_HAB_total + HAB_tank_mass);

HAB_inert_mass = ME_HAB_total + HAB_tank_mass + HAB_support_mass;

%%%%%%%%%%

% Output %

%%%%%%%%%%

fprintf(1,'%1s\n',' ');

fprintf(1,'%19s\n','Delta V Budget, m/s');

fprintf(1,'%1s\n',' ');

fprintf(1,'%32s %8.2f\n','Mars Orbit DV:',DV_orbit);

fprintf(1,'%32s %8.2f\n','Entry Angle Change DV:',DV_entry);

fprintf(1,'%32s %8.2f\n','Spin-up DV',DV_tether_sul);

fprintf(1,'%32s %8.2f\n','Enroute RCS/Maneuvering DV:',DV_RCS);

fprintf(1,'%32s %8.2f\n','Total DV:',HAB_DV);

fprintf(1,'%1s\n',' ');

fprintf(1,'%6s\n','Masses');

fprintf(1,'%32s %15s\n',' ',' Hab Subsystem');

fprintf(1,'%32s %15.2f\n','OMS Engine mass (kg): ',ME_HAB);

fprintf(1,'%32s %15.2f\n','RCS Engine mass (total) (kg): ',ME_RCS);

fprintf(1,'%32s %15.2f\n','Oxidizer tank masses (kg): ',Moxt_HAB);

fprintf(1,'%32s %15.2f\n','Fuel tank masses (kg): ',Mfut_HAB);

fprintf(1,'%32s %15.2f\n','Pressurant tank masses (kg): ',Mprt_HAB);

fprintf(1,'%32s %15.2f\n','Struct. support mass (kg): ',HAB_support_mass);

fprintf(1,'%1s\n',' ');

fprintf(1,'%32s %15.2f\n','Total inert mass (kg): ',HAB_inert_mass);

fprintf(1,'%1s\n',' ');

fprintf(1,'%32s %15.2f\n','Oxidizer masses (kg): ',ox_total);

fprintf(1,'%32s %15.2f\n','Fuel masses (kg): ',fu_total);

fprintf(1,'%32s %15.2f\n','Pressurant masses (kg): ',HAB_press_mass);

fprintf(1,'%32s %15.2f\n','Total prop/press masses (kg): ',ox_total+fu_total+HAB_press_mass);

fprintf(1,'%1s\n',' ');

fprintf(1,'%32s %15.2f\n','Total mass (kg): ',ox_total+fu_total+HAB_press_mass+HAB_inert_mass);

fprintf(1,'%1s\n',' ');

fprintf(1,'%15s\n','Tank and Engine Geometry');

fprintf(1,'%32s %15s %15s\n',' ',' Length (m)',' Diameter (m)');

fprintf(1,'%32s %15.2f %15.2f\n','OMS Engine: ',LE_HAB,DE_HAB);

fprintf(1,'%32s %15s %15.2f\n','OMS Engine Throat: ',' N/A',(pi*DE_HAB^2/4)/Ae_At);

fprintf(1,'%32s %15.2f %15.2f\n','RCS Engines on Hab :',LE_RCS,DE_RCS);

fprintf(1,'%32s %15.2f %15.2f\n','Oxidizer tank on Hab: ',Loxt_HAB+D_prop,D_prop);

fprintf(1,'%32s %15.2f %15.2f\n','Fuel tank on Hab: ',Lfut_HAB+D_prop,D_prop);

fprintf(1,'%32s %15s %15.2f\n','Pressurant tank on Hab: ',' Spherical',Dprt_HAB);

fprintf(1,'%1s\n',' ');

A.3.2

supersubhm.m

% AAE 450 Fall 2001

%

% Mars Descent Simulation

% Written by Jeremy Davis for Spring 2001 AAE 450, modified by Jon Edwards for Fall 2001

clear all

close all

warning off

global msc areapsuper numsuper cdsub areapsub numsub gamma R gmars

%Constants

msc = 47042.7; % kg Mass of the spacecraft

dsc = 9; % m Diameter of the spacecraft

areasc = pi*((dsc/2)^2); %m^2 Area of the spacecraft

gamma = 1.359; % Ratio of specific heats for CO2

R = 191.18; % J/kg*K % Universal Gas Constant for Martian atmosphere

gmars = 3.71; % Mars gravitational acceleration

superdiam

% INITIAL CONDITIONS

init_alt = 12.138; % km Initial Altitude where Mach 3 is reached

init_M = 3;

T = 0.0018*init_alt^3 - 0.0915*init_alt^2 - 0.2238*init_alt + 214.93; % deg K Temperature (Curve fit from Table 1 of Post-Viking Atmospheric Models Handout)

sos = sqrt(gamma*R*T); % m/s speed of sound for Martian atmosphere at initial altitude

init_vel = init_M*sos; % m/s

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%Initial Conditions

t(1) = 0; % seconds

x(1) = 0; % meters

y(1) = init_alt*1000; % meters

vi = init_vel; % m/s

thetaideg = -9.41; % deg Flight path angle

thetai = thetaideg*pi/180; %radians Flight path angle

vx(1) = vi*cos(thetai); %m/s Velocity in x direction

vy(1) = vi*sin(thetai); %m/s Velocity in y direction

time = max(roots([-gmars/2 vy(1) y(1)])); % Time it would take to hit surface if no parachute

p0 = [x(1) vx(1) y(1) vy(1)]; % All in SI units

ti = [0 time];

options = odeset('Reltol',1e-12);

[tt,xx] = ode45('supereom',ti,p0,options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% PERTINENT INFORMATION

disp('PERTINENT INFORMATION');

disp('------');

disp('Initial Mass of the HAB in tonnes: ');

disp(msc/1000);

disp('Altitude at Parachute Deployment in km: ');

disp(init_alt);

disp('Velocity at Parachute Deployment in m/s: ');

disp(init_vel);

disp('------');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% TRAJECTORY SOLUTIONS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x = xx(:,1); % ALL IN METERS

vx = xx(:,2);

y = xx(:,3);

vy = xx(:,4);

t = tt;

v = sqrt((vx.^2) + (vy.^2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% TRIM THE VECTORS TO REMOVE NEGATIVE ALTITUDES

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sz = size(x);

i = 1;

while i < sz(1,1);

if abs(v(i)) < 0.9*(-6E-05*y(i)^4 + 0.0051*y(i)^3 - 0.1357*y(i)^2 + 0.4464*y(i) + 234.78); % Speed of Sound Curve fit from Post Viking Atmospheric Models Handout)

brk = i - 1;

t = t(1:brk);

x = x(1:brk);

y = y(1:brk);

vx = vx(1:brk);

vy = vy(1:brk);

v = sqrt((vx.^2) + (vy.^2));

break

end

i = i + 1;

end

endofsuper=i-1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% FINDING THE ACCELERATION AND RESULTING G-LOADS DURING

% PARACHUTE DEPLOYMENT

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

gload = -1.*acceldiff(t,v)./9.81;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% PERTINENT INFORMATION

disp('Supersonic Parachute Pertinent Information');

disp('------');

disp('Number of Supersonic Parachutes Being Used: ');

disp(numsuper);

disp('Diameter of each parachute in meters is [m]: ');

disp(diamsuper);

disp('The mass of the supersonic parachutes combined are [kg]: ');

disp(Mparasuper);

disp('Supersonic Parachute Deployment Time [sec]: ');

disp(t(end));

disp('------');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% ======

% ======

% ======STAGE BREAK ======

% ======

% ======

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% SUBSONIC REEFED STAGE

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subreefdiam

reefref = .3; % reefing factor as the percentage of nominal area

surfacearea_sub=surfacearea_sub*reefref;

%Initial Conditions

reeftime = 10; % Deployment time of reefed subsonic chutes

treef(1) = t(end);

xreef(1) = x(end); % meters

yreef(1) = y(end); % meters

vxreef(1) = vx(end); % m/s

vyreef(1) = vy(end); % m/s

p0reef = [xreef(1) vxreef(1) yreef(1) vyreef(1)]; % All in SI units

tireef = [treef(1) reeftime+treef(1)];

options = odeset('Reltol',1e-12);

[ttreef,xxreef] = ode45('subreefeom',tireef,p0reef,options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% REEFED TRAJECTORY RESULTS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xreef = xxreef(:,1);

vxreef = xxreef(:,2);

yreef = xxreef(:,3);

vyreef = xxreef(:,4);

treef = ttreef;

vreef = sqrt((vxreef.^2) + (vyreef.^2));

gloadreef = -1.*acceldiff(treef,vreef)./9.81;

endofreef = length(ttreef) + endofsuper - 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% COMBINING THE SUPERSONIC AND REEFED RESULTS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x = [x' xreef']';

y = [y' yreef']';

vx = [vx' vxreef']';

vy = [vy' vyreef']';

t = [t' treef']';

v = [v' vreef']';

gload = [gload' gloadreef'];

sz = size(x);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% SUBSONIC FULLY DEPLOYED RINGSLOT PARACHUTES

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subreefdiam

%Initial Conditions

deptime = 200; % deployment time of fully deployed ringslot subsonic chutes

areapsub = numsub*(pi*(diamsub/2)^2);

tsub(1) = t(end);

xsub(1) = x(end); % meters

ysub(1) = y(end); % meters

vxsub(1) = vx(end); %m/s

vysub(1) = vy(end); %m/s

p0sub = [xsub(1) vxsub(1) ysub(1) vysub(1)]; % All in SI units

tisub = [tsub(1) deptime+tsub(1)];

options = odeset('Reltol',1e-12);

[ttsub,xxsub] = ode45('subreefeom',tisub,p0sub,options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% FULLY DEPLOYED CHUTE TRAJECTORY RESULTS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xsub = xxsub(:,1);

vxsub = xxsub(:,2);

ysub = xxsub(:,3);

vysub = xxsub(:,4);

tsub = ttsub;

vsub = sqrt((vxsub.^2) + (vysub.^2));

gloadsub = -1.*acceldiff(tsub,vsub)./9.81;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% CUT OFF TIMES AFTER TERMINAL VELOCITY IS REACHED

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sz = size(tsub);

z = 2;

while z < sz(1,1)

if gloadsub(z) <= 0.05 % Ringslots are cut when the acceleration of the vehicle falls below 0.05 m/s^2 which is close to terminal velocity

termpoint = z;

xsub = xsub(1:termpoint);

ysub = ysub(1:termpoint);

vxsub = vxsub(1:termpoint);

vysub = vysub(1:termpoint);

vsub = vsub(1:termpoint);

tsub = tsub(1:termpoint);

gloadsub = gloadsub(1:termpoint);

sos_deploy = (-6E-05*ysub(end)^4 + 0.0051*ysub(end)^3 - 0.1357*ysub(end)^2 + 0.4464*ysub(end) + 234.78);

break

end

z = z + 1;

end

endofsub = length(tsub) + endofreef - 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% COMBINING THE SUPERSONIC, REEFED, AND SUBSONIC RESULTS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x = [x' xsub']';

y = [y' ysub']';

vx = [vx' vxsub']';

vy = [vy' vysub']';

t = [t' tsub']';

v = [v' vsub']';

gload = [gload gloadsub'];

sz = size(x);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% PERTINENT INFORMATION

disp('Fully Deployed Subsonic Chute Pertinent Information');

disp('------');

disp('Number of Subsonic Parachutes: ');

disp(numsub);

disp('Diameter of Each Subsonic Parachute [m]: ');

disp(diamsub);

disp('The mass of the subsonic parachutes combined are [kg]: ');

disp(Mparasub);

disp('Deployment Time of Subsonic Parachutes [sec]: ');

disp(tsub(end)-tsub(1));

disp('Velocity at 1500 meters altitude when parafoil is deployed [m/s]: ');

disp(v(size(v,1)));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% PARAFOIL PART OF DESCENT

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

parafoildiam

%Initial Conditions

deptime_para = 200; % deployment time of fully deployed subsonic chutes

tpara(1) = t(end); % seconds

xpara(1) = x(end); % meters

ypara(1) = y(end); % meters

vxpara(1) = vx(end); % m/s

vypara(1) = vy(end); % m/s

p0para = [xpara(1) vxpara(1) ypara(1) vypara(1)]; % All in SI units

tipara = [tpara(1) deptime_para+tpara(1)];

options = odeset('Reltol',1e-12);

[ttpara,xxpara] = ode45('parafoileom',tipara,p0para,options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% PARAFOIL TRAJECTORY RESULTS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xpara = xxpara(:,1);

vxpara = xxpara(:,2);

ypara = xxpara(:,3);

vypara = xxpara(:,4);

tpara = ttpara;

vpara = sqrt((vxpara.^2) + (vypara.^2));

gloadpara = -1.*acceldiff(tpara,vpara)./9.81;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% CUT OFF TIMES AFTER TERMINAL VELOCITY IS REACHED

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sz = size(tpara);

z = 2;

while z < sz(1,1)

if ypara(z) <= 0 % At 0 altitude the spacecraft has landed

termpoint = z;

xpara = xpara(1:termpoint);

ypara = ypara(1:termpoint);

vxpara = vxpara(1:termpoint);

vypara = vypara(1:termpoint);

vpara = vpara(1:termpoint);

tpara = tpara(1:termpoint);

gloadpara = gloadpara(1:termpoint);

break

end

z = z + 1;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% COMBINING THE SUPERSONIC, REEFED, SUBSONIC, AND PARAFOIL RESULTS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x = [x' xpara']';

y = [y' ypara']';

vx = [vx' vxpara']';

vy = [vy' vypara']';

t = [t' tpara']';

v = [v' vpara']';

gload = [gload gloadpara'];

sz = size(x);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% PARAFOIL PERTINENT INFORMATION

disp('Parafoil Pertinent Information');

disp('------');

disp('Span of the Parafoil [m]: ');

disp(b);

disp('The mass of the parafoil as ratio [kg]: ');

disp(Mparafoil_ratio);

disp('The mass of the parafoil as calculated [kg]: ');

disp(Mparafoil);

disp('Deployment Time of Parafoil [sec]: ');

disp(tpara(end)-tpara(1));

disp('Time to land [sec]:');

disp(t(end));

disp('Max G-loading [m/s^2]: ');

disp(max(gload));

disp('Velocity at Touchdown [m/s]: ');

disp(v(size(v,1)));

disp('Vertical Velocity at Touchdown [m/s]: ');

disp(vy(size(v,1)));

disp('Horizontal Velocity at Touchdown [m/s]: ');

disp(vx(size(v,1)));

disp('Total Mass of Parachutes [kg]: ');

total_mass = Mparasuper+Mparasub+Mparafoil_ratio;

disp(total_mass);

fprintf('%15s %15s %15s\n','Supersonic','Subsonic','Parafoil');

fprintf('%15.2f %15.2f %15.2f\n',Mparasuper,Mparasub,Mparafoil_ratio);

sz = size(t);

figure(1); hold on;

plot(t,y,'b'); axis([0 max(t)+10 0 init_alt*1000+1000]);plot(t(endofsuper),y(endofsuper),'kx'); plot(t(endofreef),y(endofreef),'kx'); plot(t(endofsub),y(endofsub),'kx');

xlabel('Time [sec]');

ylabel('Altitude [meters]');title('Altitude vs. Time');

figure(2); hold on;

plot(t,v,'b');plot(t(endofsuper),v(endofsuper),'kx'); plot(t(endofreef),v(endofreef),'kx'); plot(t(endofsub),v(endofsub),'kx');

xlabel('Time [sec]');

ylabel('Velocity [m/s]');title('Velocity vs. Time')

figure(3); hold on;

plot(t(1:sz(1,1)-2),gload,'b'); plot(t(endofsuper),gload(endofsuper),'kx'); plot(t(endofreef),gload(endofreef),'kx'); plot(t(endofsub),gload(endofsub),'kx');

xlabel('Time [sec]');

ylabel('Gload');title('G-load vs. Time');

figure(4); hold on;

plot(x,y,'b'); axis([0 max(x)+1000 0 18000]); plot(x(endofsuper),y(endofsuper),'kx'); plot(x(endofreef),y(endofreef),'kx'); plot(x(endofsub),y(endofsub),'kx');

xlabel('horizontal position (m)')

ylabel('altitude (m)')

title('Parachute Trajectory')

figure(5); hold on;

plot(t,abs(vy),'b'); plot(t(endofsuper),abs(vy(endofsuper)),'kx'); plot(t(endofreef),abs(vy(endofreef)),'kx'); plot(t(endofsub),abs(vy(endofsub)),'kx');

xlabel('time (s)')

ylabel('vertical velocity (m/s)')

title('Vertical Velocity vs. Time')

A.3.3

superdiam.m

%

% AAE 450 Fall 2001

% Written by Jeremy Davis for Spring '01, modified by Jon Edwards for Fall '01

% Hemsiflo Ribbon Parachute parameters and mass calculation

global msc numsuper dsc surfacearea

diamsuper = 30; % meters Constructed diameter of each chute

numsuper = 1; % Number of parachutes being used

Dc_Do = 0.62; % Constructed diameter to nominal diameter ratio of Hemisflo Ribbon Parachute (Knacke Table 5-2)

nom_diameter = diamsuper/Dc_Do; % m Nominal diameter of Hemisflo Ribbon Parachute

surfacearea = (nom_diameter/1.1284)^2; % m^2 Total canopy surface area (Knacke p. 5-2)

circumsuper = pi*diamsuper; % m Circumference of each parachute

areapsuper = pi*(diamsuper/2)^2; % m^2 Area of each chute

% CANOPY MASS

% ------

dens = 0.0373; % kg/m^2 Specific canopy weight (Knacke p. 6-98)

Mcsuper = dens*surfacearea; % kg Mass of the parachute canopy

% SUSPENSION LINE MASS

% ------

lodsuper = 6;% Rigging length over diameter of parachute (Knacke 5-99)

lsussuper = lodsuper*dsc; % meters Length of suspension lines (Knacke 5-99)

nsuper = round(circumsuper/4); % Number of rigging lines

densl = .01488; % kg/m Linear density of nylon suspension lines (Knacke 6-99)

Mslsuper = numsuper*nsuper*lsussuper*densl; % kg Mass of the suspension lines

Mparasuper = Mcsuper + Mslsuper; % kg Mass of the total parachute system

A.3.4

subreefdiam.m

%

% AAE 450 Fall 2001

% Written by Jeremy Davis for Spring '01, modified by Jon Edwards for Fall '01

% We are using ringslot parachutes in this analysis since according to Knacke (Table 5-2)

% ringslot parachutes are used in applications where 0.1 < M < 0.9

global msc cdsub surfacearea_sub numsub dsc

cdsub = .60; % Drag coefficient of Ringslot Parachute (Knacke p 5-29)

diamsub = 70; % m Diameter of each chute

Dc_Do_sub = 1.00; % Constructed diameter to nominal diameter ratio of ringslot parachute (Knacke Table 5-2)

nom_diameter_sub = diamsub/Dc_Do_sub; % m Nominal diameter of ringslot parachute

surfacearea_sub = (nom_diameter_sub/1.1284)^2; % m^2 Total canopy surface area (Knacke p. 5-2)

numsub = 3; % Number of parachutes being used

circumsub = 2*pi*(diamsub/2); % meters Circumference of each chute

areapsub = pi*(diamsub/2)^2; % meters^2 Area of the parachute

% CANOPY MASS

% ------

dens = 0.0373; % kg/m^2 Specific canopy weight (Knacke p. 6-98)

Mcsub = numsub*dens*surfacearea_sub; % kg Mass of all the ringslot parachute canopies

% SUSPENSION LINE MASS

% ------

lod = 6;% Rigging length over diameter of parachute (Six-forebody-diameter rule, Knacke 5-21)

lsussub = lod*dsc; % meters, length of suspension lines

nsub = round(circumsub/4); % Number of rigging lines

densl = .01488; % kg/m, Linear density of nylon (Knacke 6-99)

Mslsub = numsub*nsub*lsussub*densl; % kg Mass of the suspension lines

Mparasub = Mcsub + Mslsub; % kg Mass of the total parachute system

A.3.5

parafoildiam.m

%

% AAE 450 Fall 2001

% Written by Jon Edwards for Fall '01

global areaparafoil Cd_parafoil Cl_parafoil msc Mparafoil alpha_rad

Lod = 3.2; % Parafoil maneuverable parachute (Knacke fig 5-106 Curve C at Best Glide Point)

Cl_parafoil = 0.43; % (Knacke fig 5-106 Best Glide Point)

Cd_parafoil = Cl_parafoil/Lod; % (Knacke fig 5-106 Best Glide Point)

alpha_rad = atan(Cd_parafoil/Cl_parafoil); % rad, fixed angle of attack of parafoil (Knacke fig. 5-105)

alpha_deg = alpha_rad*(180/pi); % deg, fixed angle of attack of parafoil

AR = 2.7; % Aspect Ratio of X-38 Parafoil Landing System (AIAA-99-1730 p. 206)

b = 176; % meters, Span of Parafoil (variable parameter set by user)

c = b/AR; % meters, Chord of Parafoil

areaparafoil = b*c; % m^2 Area of Parafoil

numcells = round(areaparafoil*(31/510.97)); % Number of cells (proportional to X-38 Parafoil Landing System)

thickness = .148*c; % m Thickness of Parafoil (proportional to X-38 parafoil landing system)

totalarea = areaparafoil*2 + (numcells+1)*thickness*c; % m^2 Total surface area of parafoil

Mparafoil_ratio = msc*(390/11340); % kg Mass estimate (ratio from X-38 parafoil landing system)

% CANOPY MASS

% ------

dens = .0373; % in kg/m^2 ("specific weight" of nylon)

Mcparafoil = totalarea*dens; % in kg

% SUSPENSION LINE MASS

% ------

lob = 0.62; % Line length to span ratio (from AIAA-99-1730 p. 207)

linelength = lob*b; % m Line length

numlines = (12/13.72)*c*(numcells+1); % Number of suspension lines (based on X-38 parafoil)

kevlar = 0.00521; % kg/m Kevlar specific weight (Knacke p 6-99)

Mslparafoil = numlines*linelength*kevlar; % Mass of suspension lines

Mparafoil = Mcparafoil + Mslparafoil; % kg Total mass of the parafoil

A.3.6

supereom.m

%

% A&AE 450 Fall 2001

% Written by Jeremy Davis Spring '01, modified by Jon Edwards Fall '01

% Code used to solve the ODE's of the vehicle during the

% supersonic (M > 0.9) stages of deployment.

% Hemisflo Ribbon Supersonic Parachutes are used in this analysis

function xdot = supreefeom(t,x)

global msc numsuper gamma R gmars surfacearea

v = sqrt(x(2)^2 + x(4)^2); % m/s Velocity of the spacecraft

ro = 0.016*exp(-0.092*(x(3)/1000)); % kg/m^3 Density (Table 1 of Post-Viking Atmospheric Models handout)

q = .5*ro*(v.^2); % N/m^2 Dynamic Pressure

temp = 0.0018*x(3)^3 - 0.0915*x(3)^2 - 0.2238*x(3) + 214.93; % deg K Temperature (Table 1 of Post-Viking Atmospheric Models Handout)

mach = v/sqrt(gamma*R*temp); % Mach No.

if mach > 1.8

cdsuper = -0.1983*mach+0.645; % Curve fit from Figure 5-96 in Knacke for M > 1.8

else

cdsuper = 0.45; % Figure 5-96 from Knacke for 0.9 < M < 1.8

end

drag = numsuper*cdsuper*surfacearea*q; % N Drag (Knacke 4-9)

xdot(1) = x(2);

xdot(2) = -((drag*x(2))/(v*msc));

xdot(3) = x(4);

xdot(4) = -((drag*x(4))/(v*msc)) - gmars;

xdot = xdot';

return

A.3.7

subreefeom.m

%

% A&AE 450 Fall 2001

% Written by Jeremy Davis Spring '01, modified by Jon Edwards Fall '01

% Code used to solve the ODE's of the vehicle during the

% subsonic stages of deployment.

function xdot = subreefeom(t,x)

global msc reefref cdsub surfacearea_sub numsub gmars

v = sqrt(x(2)^2 + x(4)^2); % m/s Velocity of the spacecraft

ro = 0.016*exp(-0.092*(x(3)/1000)); % kg/m^3 Density (Curve fit from Table 1 of Post-Viking Atmospheric Models handout)

q = .5*ro*(v^2); % N/m^2 Dynamic Pressure

dragsub = numsub*cdsub*surfacearea_sub*q; % N Drag (Knacke 4-9)

xdot(1) = x(2);

xdot(2) = -(dragsub*x(2))/(v*msc);

xdot(3) = x(4);

xdot(4) = -(dragsub*x(4))/(v*msc) - gmars;

xdot = xdot';

return

A.3.8

parafoileom.m

%

% A&AE 450 Fall 2001

% Written by by Jon Edwards Fall '01

% Code used to solve the ODE's of the vehicle during the

% parafoil stage of deployment.

function xdot = parafoileom(t,x)

global Cd_parafoil areaparafoil msc Cl_parafoil gmars alpha_rad

v = sqrt(x(2)^2 + x(4)^2); % m/s Velocity of the spacecraft

ro = 0.016*exp(-0.092*(x(3)/1000)); % kg/m^3 Density (Table 1 of Post-Viking Atmospheric Models handout)

q = .5*ro*(v.^2); % N/m^2 Dynamic Pressure

drag = Cd_parafoil*areaparafoil*q; % N Drag (Knacke 4-9)

lift = Cl_parafoil*areaparafoil*q; % N Lift

flight_angle = acos(x(2)/v);

xdot(1) = x(2);

xdot(2) = lift*sin(flight_angle)/msc-drag*cos(flight_angle)/msc;

xdot(3) = x(4);

xdot(4) = lift*cos(flight_angle)/msc+drag*sin(flight_angle)/msc - gmars;

xdot = xdot';

return

A.3.9

acceldiff.m

%

% A&AE 450 Fall 2001

% Code used to numerically differentiate velocity vectors

% in order to obtain acceleration histories.

% Written by Jeremy Davis for Spring 2001, modified by Jon Edwards for Fall 2001

function accel = acceldiff(t,v)

sz = size(t);

i = 2;

accel(1,1) = 0;

while i < sz(1,1)

accel(i) = (v(i) - v(i-1))/(t(i) - t(i-1));

i = i + 1;

end

size(accel);

accel(1,1) = accel(1,2);

accel = accel';

return

A-3-1