Andrew L. and Andrew Z.
Computational Physics
Problem Statement
There are numerous (or innumerable) methods for modeling and considering characteristics of rocket trajectory. These models present advantages in simplicity or complexity depending on the intended model application. A trajectory model intended for a commercial customer may be exceptionally resource intensive in set-up cost or processing resources but may be necessary depending on the computational design basis’s accuracy, precision, and stability. A less process-intensive model may be simpler to compute and present less initial setup investment but may not meet the requirements for much beyond a theoretical application.
We intend to model rocket trajectories considering various computational models. We intend to compare variations in Deterministic and Probabilistic computational models for rocket trajectories and identify models that present value in processing efficiency, flexibility, and exhaustiveness in analysis.
Modeling the Atmosphere
Build a model to identify atmospheric temperature as a function of z-altitude, add randomized variance
Build a model to identify atmospheric pressure as a function of z-altitude
Build a model to identify atmospheric density as a function of z-altitude
Consider the expected sparse array dimensions from initial conditions.
References
http://www.braeunig.us/space/atmmodel.htm
clear all;
% init_param specifies the initial parameters for everything in two sets,
% so that we can compare the sets to each other on a plot.
%
% Format: v_init (m/sec), rho_0 (kg/m^3), A (m^2), m_shot (dens*vol...kg), theta (rad),
% dt (sec)
% v_init = initial cannonball speed
% rho_0 = reference air density
% A = cross-sectional area of cannonball
% m_shot = mass of cannonball
% theta = initial angle of trajectory (rad)
% phi = angle to y as measured from x (rad)
% dt = time window
% FYI: 1 m/sec ~ 2.2 mph
% ... and these are the parameters (you can use one or more...ncompare
% allows you to choose how many of these sets of parameters you will plot.
for ii=1:8;
init_param(1,:)=[0.0001, 1.225, pi*0.1^2, 7870*(4/3)*pi*(0.1)^3, 45*pi/180, (ii-1)*90*pi/180, 0.01];
cannon_3d_1(init_param);
end
Notes and Functions:
Describe the Function
function cannon_3d_1(init_param)
%
% 3D Cannonball with drag, variable density and drag coeffient...but no
% wind.
% This program allows you to compare trajectories given ncompare sets of parameters.
%
%%%%%%%%% Initialization part %%%%%%%%%%%%%%%
% This sets the number of parameter sets that you can compare to each other
ncompare=size(init_param,1);
clr=['k' 'r' 'b' 'm' 'g' 'c'];
% Physical parameters
%
g=9.8; % g=9.8 m/sec^2
zref=1e4; %Altitude reference for air density (m)
% Initial positions...keep these the same
x_init=0;
y_init=0;
z_init=0;
t_init=0;
vecc=[];
vel=@(v) sqrt(v(1).^2+v(2).^2+v(3).^2);
for icomp=1:ncompare %Loop over <ncompare> parameter changes
v_init =init_param(icomp,1);
rho_0 =init_param(icomp,2);
A =init_param(icomp,3);
m_shot =init_param(icomp,4);
theta =init_param(icomp,5);
phi =init_param(icomp,6);
dt =init_param(icomp,7);
Fthr = 10000; %The Thrust Force in Newtons
tburn = 10; % Time it takes to burn the fuel in seconds
stges = 1; %Number of stages rockets
x=[];
y=[];
z=[];
vx=[];
vy=[];
vz=[];
tim=[];
vtot=[];
x(1) = x_init;
y(1) = y_init;
z(1) = z_init;
vx(1) = v_init*cos(theta)*cos(phi);
vy(1) = v_init*cos(theta)*sin(phi);
vz(1) = v_init*sin(theta);
Fthrx(1) = Fthr*cos(theta)*cos(phi);
Fthry(1) = Fthr*cos(theta)*sin(phi);
Fthrz(1) = Fthr*sin(theta);
tim(1) = t_init;
vtot(1)=vel([vx(1) vy(1) vz(1)]);
C=drag_coeff(rho_0*exp(-z(1)/zref),0.1,vtot(1));
%%%%%%%%% Euler part %%%%%%%%%%%%%%%
zz=100; % just to make the loop work the first time around
step=1;
while(zz>0) % stop loop if altitude goes below 0 again
x(step+1)=x(step)+vx(step)*dt;
y(step+1)=y(step)+vy(step)*dt;
z(step+1)=z(step)+vz(step)*dt;
zz=z(step+1);
fx=-0.5*C*rho_0*exp(-z(step)/zref)*A*vx(step)*vtot(step)/m_shot + Fthrx/m_shot;
fy=-0.5*C*rho_0*exp(-z(step)/zref)*A*vy(step)*vtot(step)/m_shot + Fthry/m_shot;
fz=-0.5*C*rho_0*exp(-z(step)/zref)*A*vz(step)*vtot(step)/m_shot - g + Fthrz/m_shot;
vx(step+1) = vx(step) + fx*dt;
vy(step+1) = vy(step) + fy*dt;
vz(step+1) = vz(step) + fz*dt;
vtot(step+1)=vel([vx(step+1) vy(step+1) vz(step+1)]);
if tim(step) < tburn
m_shot = m_shot - ((0.15)/tburn)*dt*m_shot;
else
Fthrx = 0;
Fthry = 0;
Fthrz = 0;
end
tim(step+1) = tim(step) + dt;
step=step+1;
C=drag_coeff(rho_0*exp(-z(step)/zref),0.1,vtot(step));
%m_shot=mass_adj(m_shot,)
end
% Plot both numerical and analytical solution
figure(1);
hold on;
if(icomp==ncompare)
for idd=1:ncompare
plot(0,0,clr(idd),'LineWidth',5);
strr{idd}=sprintf('params #%d',idd);
end
gg=legend(strr{:},'Location','SW');
set(gg,'FontSize',16);
end
plot(x,z,clr(icomp),'LineWidth',2);
xlabel('x (m)','FontSize',16);
ylabel('z (m)','FontSize',16);
title('X-Z Trajectory');
hold off;
% Y-trajectory plot
figure(2);
if(icomp>1)
hold on;
end
plot(y,z,clr(icomp),'LineWidth',2);
xlabel('y (m)','FontSize',16);
ylabel('z (m)','FontSize',16);
title(' Y-Z Trajectory');
hold off;
% Z-trajectory plot
figure(3);
if(icomp>1)
hold on;
end
% Plot both numerical and analytical solution
semilogy(tim,vx,[clr(icomp) ':'], ...
tim,vy,[clr(icomp) '--'], ...
tim,abs(vz),clr(icomp),'LineWidth',2);
gg=legend('vx','vy','|vz|','Location','SE');
set(gg,'FontSize',16);
ylim([1e-2 1e4]);
xlabel('time (sec)','FontSize',16);
ylabel('velocity (m/sec)','FontSize',16);
title('Velocity Profiles vs. Time');
text(500,50,sprintf('Final x=%.2f m',x(end)), ...
'FontSize',16);
text(500,200,sprintf('Final time=%.2f sec',tim(end)), ...
'FontSize',16);
hold off;
% 3D-trajectory plot
figure(4);
if(icomp>1)
hold on;
end
% Plot both numerical and analytical solution
plot3(x,y,z,clr(icomp),'LineWidth',2);
xlabel('x (m)','FontSize',16);
ylabel('y (m)','FontSize',16);
zlabel('z (m)','FontSize',16);
vecc=[vecc x(end) y(end) z(end)];
% if(icomp==ncompare)
% maxx=max(vecc);
% xlim([-maxx maxx]);
% ylim([-maxx maxx]);
% zlim([0 maxx]);
% end
title('Trajectory');
hold off;
end
end
Drag
function cd = drag_coeff( dens, L, v )
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
re_num=dens*v*L/1.75e-5;
cd(1:length(re_num))=1.5;
idx=find(re_num<1e5 & re_num>2);
if(length(idx)>0)
cd(idx)=0.5+exp(-(re_num(idx))/300);
end
idx=find(re_num<6e5 & re_num>=1e5);
if(length(idx)>0)
cd(idx)=0.5-0.4*(1./(1+exp(-3e-5*(re_num(idx)-3e5))));
end
idx=find(re_num>=6e5);
if(length(idx)>0)
cd(idx)=0.1+0.35*(1-exp(-(re_num(idx)-6e5)/3e6));
end
end
Augusta Physics 