Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- % Program:
- % rocket-trajectory-simulation.m
- % Multi-stage rocket dynamics and trajectory simulation.
- %
- % Description:
- % Predicts multi-stage rocket dynamics and trajectories based on the given
- % rocket mass, engine thrust, launch parameters, and drag coefficient.
- %
- % Variable List:
- % Delta = Time step (s)
- % t = Time (s)
- % Thrust = Thrust (N)
- % Mass = Mass (kg)
- % Mass_Rocket_With_Motor = Mass with motor (kg)
- % Mass_Rocket_Without_Motor = Mass without motor (kg)
- % Theta = Angle (deg)
- % C = Drag coefficient
- % Rho = Air density (kg/m^3)
- % A = Rocket projected area (m^2)
- % Gravity = Gravity (m/s^2)
- % Launch_Rod_Length = Length of launch rod (m)
- % n = Counter
- % Fn = Normal force (N)
- % Drag = Drag force (N)
- % Fx = Sum of forces in the horizontal direction (N)
- % Fy = Sum of forces in the vertical direction (N)
- % Vx = Velocity in the horizontal direction (m/s)
- % Vy = Velocity in the vertical direction (m/s)
- % Ax = Acceleration in the horizontal direction (m/s^2)
- % Ay = Acceleration in the vertical direction (m/s^2)
- % x = Horizontal position (m)
- % y = Vertical position (m)
- % Distance_x = Horizontal distance travelled (m)
- % Distance_y = Vertical travelled (m)
- % Distance = Total distance travelled (m)
- % Memory_Allocation = Maximum number of time steps expected
- clear, clc % Clear command window and workspace
- % Parameters
- Delta = 0.001; % Time step
- Memory_Allocation = 30000; % Maximum number of time steps expected
- % Preallocate memory for arrays
- t = zeros(1, Memory_Allocation);
- Thrust = zeros(1, Memory_Allocation);
- Mass = zeros(1, Memory_Allocation);
- Theta = zeros(1, Memory_Allocation);
- Fn = zeros(1, Memory_Allocation);
- Drag = zeros(1, Memory_Allocation);
- Fx = zeros(1, Memory_Allocation);
- Fy = zeros(1, Memory_Allocation);
- Ax = zeros(1, Memory_Allocation);
- Ay = zeros(1, Memory_Allocation);
- Vx = zeros(1, Memory_Allocation);
- Vy = zeros(1, Memory_Allocation);
- x = zeros(1, Memory_Allocation);
- y = zeros(1, Memory_Allocation);
- Distance_x = zeros(1, Memory_Allocation);
- Distance_y = zeros(1, Memory_Allocation);
- Distance = zeros(1, Memory_Allocation);
- C = 0.556; % Drag coefficient
- Rho = 1.204; % Air density (kg/m^3)
- A = 6.76*10^-2; % Rocket projected area (m^2)
- Gravity = 9.81; % Gravity (m/s^2)
- Launch_Rod_Length = 5; % Length of launch rod (m)
- Mass_Rocket = 239.034; % Mass with propllant (kg)
- Mass_Rocket_Empty = 138.384; % Mass without propllant (kg)
- Mass_Propellant = 100.65; % Mass of propellant (kg)
- Mass_Flow_Rate = 14.4; % Mass Loss of propellant (kg/s)
- Burn_Time = 6.985; % Rocket Burn Time (s)
- Engine_Max_Thrust = 29130; % Max Thrust of rocket (N)
- Theta(1) = 70; % Initial launch angle (deg)
- Vx(1) = 0; % Initial horizontal speed (m/s)
- Vy(1) = 0; % Initial vertical speed (m/s)
- x(1) = 0; % Initial horizontal position (m)
- y(1) = 0.1; % Initial vertical position (m)
- Distance_x(1) = 0; % Initial horizontal distance travelled (m)
- Distance_y(1) = 0; % Initial vertical distance travelled (m)
- Distance(1) = 0; % Initial distance travelled (m)
- Mass(1) = Mass_Rocket; % Initial rocket mass (kg)
- n = 1; % Initial time step
- while y(n) > 0 % Run until rocket hits the ground
- n = n+1; % Increment time step
- t(n)= (n-1)*Delta; % Elapsed time
- % Determine rocket thrust and mass based on launch phase
- if t(n) <= 0.1 % Rocket Startup
- Thrust(n) = (Engine_Max_Thrust*10)*t(n);
- Mass(n) = Mass_Rocket;
- elseif t(n) >= 0.1 && t(n) < Burn_Time % Rocket Boost Phase
- Thrust(n) = Engine_Max_Thrust;
- Mass(n) = Mass_Rocket - (t(n)*Mass_Flow_Rate);
- elseif t(n) >= Burn_Time % Rocket Shutdown
- Thrust(n) = 0;
- Mass(n) = Mass_Rocket_Empty;
- end
- % Normal force calculations
- if Distance(n-1) <= Launch_Rod_Length % Launch rod normal force
- Fn(n) = Mass(n)*Gravity*cosd(Theta(1));
- else
- Fn(n) = 0; % No longer on launch rod
- end
- % Drag force calculation
- Drag(n)= 0.5*C*Rho*A*(Vx(n-1)^2+Vy(n-1)^2); % Calculate drag force
- % Sum of forces calculations
- Fx(n)= Thrust(n)*cosd(Theta(n-1))-Drag(n)*cosd(Theta(n-1))-Fn(n)*sind(Theta(n-1)); % Sum x forces
- Fy(n)= Thrust(n)*sind(Theta(n-1))-(Mass(n)*Gravity)-Drag(n)*sind(Theta(n-1))+Fn(n)*cosd(Theta(n-1)); % Sum y forces
- % Acceleration calculations
- Ax(n)= Fx(n)/Mass(n); % Net accel in x direction
- Ay(n)= Fy(n)/Mass(n); % Net accel in y direction
- % Velocity calculations
- Vx(n)= Vx(n-1)+Ax(n)*Delta; % Velocity in x direction
- Vy(n)= Vy(n-1)+Ay(n)*Delta; % Velocity in y direction
- % Position calculations
- x(n)= x(n-1)+Vx(n)*Delta; % Position in x direction
- y(n)= y(n-1)+Vy(n)*Delta; % Position in y direction
- % Distance calculations
- Distance_x(n) = Distance_x(n-1)+abs(Vx(n)*Delta); % Distance in x
- Distance_y(n) = Distance_y(n-1)+abs(Vy(n)*Delta); % Distance in y
- Distance(n) = (Distance_x(n)^2+Distance_y(n)^2)^(1/2); % Total distance
- % Rocket angle calculation
- Theta(n)= atand(Vy(n)/Vx(n)); % Angle defined by velocity vector
- end
- figure('units','normalized','outerposition',[0 0 1 1]) % Maximize plot window
- % Figure 1
- subplot(3,2,1)
- plot(x(1:n),y(1:n));
- xlim([0 inf]);
- ylim([0 inf]);
- xlabel({'Range (m)'});
- ylabel({'Altitude (m)'});
- title({'Height-Range Graph (Trajectory)'});
- % Figure 2
- subplot(3,2,2)
- plot(t(1:n),Vx(1:n));
- xlabel({'Time (s)'});
- ylabel({'Vx (m/s)'});
- title({'Vertical Velocity Graph'});
- % Figure 3
- subplot(3,2,3)
- plot(t(1:n),Vy(1:n));
- xlabel({'Time (s)'});
- ylabel({'Vy (m/s)'});
- title({'Horizontal Velocity Graph'});
- % Figure 4
- subplot(3,2,4)
- plot(t(1:n),Mass(1:n));
- xlim([0 60]);
- ylim([0 240]);
- xlabel({'Time (s)'});
- ylabel({'Mass (kg)'});
- title({'Rocket Mass Graph'});
- % Figure 5
- subplot(3,2,5)
- plot(t(1:n),((Ax(1:n) + Ay(1:n))/2));
- xlabel({'Time (s)'});
- ylabel({'Acceleration (m/s^2)'});
- title({'Acceleration-Time Graph'});
- % Figure 6
- subplot(3,2,6)
- plot(t(1:n),y(1:n));
- xlim([0 inf]);
- ylim([0 inf]);
- xlabel({'Time (s)'});
- ylabel({'Altitude (m)'});
- title({'Altitude-Time Graph'});
- % Create a table with the data and variable names
- T = table(Theta, 'VariableNames', {'Rocket Angle'} )
- % Write data to text file
- writetable(T, 'MyFile.txt')
- disp(Theta(n))
Add Comment
Please, Sign In to add comment