ESS_Flight_Dynamics
/
2-DoF-TrajectorySim
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Migration to ESS Flight Dynamics.

master
Jason Ong 2 years ago
parent
d8ffb9c304
commit
8636b276c6
18 changed files with 6362 additions and 1 deletions
Unified View
Diff Options
Show Stats Download Patch File Download Diff File
  1. BIN
      2-DoF vs Astos 6-Dof.png
  2. 25
      Code/DragInterpolations_ICANSAT.m
  3. 76
      Code/DragInterpolations_MAD.m
  4. 12
      Code/DragInterpolations_Spark.m
  5. 6
      Code/ICANSAT/ICANSAT_DragData.txt
  6. 64
      Code/ICANSAT/V0/ICANSAT.m
  7. 202
      Code/ICANSAT/V0/ICANSAT_Thrust_v_Time.txt
  8. 31
      Code/ICANSAT/V1/ICANSAT_V1.m
  9. 82
      Code/ICANSAT/V1/ICANSAT_V1_Thrust_v_Time.txt
  10. 43
      Code/MAD/MAD v 4.2B_Thrust_v_Time.txt
  11. 11
      Code/MAD/MAD_DragData.txt
  12. 2695
      Code/MAD/MAD_Nominal_Inputs.txt
  13. 2603
      Code/MAD/MAD_Trajectory_Output.txt
  14. 52
      Code/MAD/MAD_Validation.m
  15. 441
      Code/SingleStageRocket.m
  16. 8
      Code/Spark/Spark_DragData.txt
  17. 12
      README.md
  18. BIN
      zjd11_Thesis_Paper_Final.docx

BIN
2-DoF vs Astos 6-Dof.png
View File

Before After
Width: 1457  |  Height: 762  |  Size: 61 KiB

25
Code/DragInterpolations_ICANSAT.m
View File

@ -0,0 +1,25 @@
function [coeffdrag] = DragInterpolations_ICANSAT(velocity, altitude, time, tb)
% global tb
drag_data = dlmread("ICANSAT/ICANSAT_DragData.txt", "\t",1,0);
% Mach No
mach_no = drag_data(:,1);
% Full base drag
cd_full = drag_data(:,3);
% No base drag
cd_zero = drag_data(:,2);
[temp, sound_speed, P, rho] = atmosisa(altitude);
mach_speed = velocity/sound_speed;
if (time > tb)
coeffdrag = interp1(mach_no, cd_full, mach_speed);
else
coeffdrag = interp1(mach_no, cd_zero, mach_speed);
end
end

76
Code/DragInterpolations_MAD.m
View File

@ -0,0 +1,76 @@
function [coeffdrag] = DragInterpolations_MAD(velocity, altitude, time, tb)
drag_data = dlmread("MAD/MAD_DragData.txt", "\t",1,0);
% Mach No
mach_no = drag_data(:,1);
% Altitude 0m, full base drag
cdfull_0 = drag_data(:,4);
% Altitude 10000m, full base drag
cdfull_10000 = drag_data(:,5);
% Altitude 30000m, full base drag
cdfull_30000 = drag_data(:,6);
% Altitude 45000m, full base drag
cdfull_45000 = drag_data(:,7);
% Altitude 60000m, full base drag
cdfull_60000 = drag_data(:,8);
% Altitude 0m, no base drag
cdzero_0 = drag_data(:,2);
% Altitude 10000m, no base drag
cdzero_10000 = drag_data(:,3);
% Interpolate for Cd values at two stated altitudes, then interpolate for
% Cd value at a specific altitude
[temp, sound_speed, P, rho] = atmosisa(altitude);
mach_speed = velocity/sound_speed;
% If propellant is still available
if (time > tb)
if (altitude <= 10000)
coeffdrag = interpolateCd(0, 10000, mach_no, cdfull_0, ...
cdfull_10000, mach_speed, altitude);
end
if (altitude <= 30000 && altitude > 10000)
coeffdrag = interpolateCd(10000, 30000, mach_no, cdfull_10000, ...
cdfull_30000, mach_speed, altitude);
end
if (altitude <= 45000 && altitude > 30000)
coeffdrag = interpolateCd(30000, 45000, mach_no, cdfull_30000, ...
cdfull_45000, mach_speed, altitude);
end
if (altitude <= 60000 && altitude > 45000)
coeffdrag = interpolateCd(45000, 60000, mach_no, cdfull_45000, ...
cdfull_60000, mach_speed, altitude);
end
if (altitude > 60000)
coeffdrag = interp1(mach_no, cdfull_60000, mach_speed);
end
else
if (altitude <= 10000)
coeffdrag = interpolateCd(0, 10000, mach_no, cdzero_0, ...
cdzero_10000, mach_speed, altitude);
end
if (altitude > 10000)
coeffdrag = interp1(mach_no, cdzero_10000, mach_speed);
end
end
end
function Cd_specific_alt = interpolateCd(low_alt, high_alt, mach_num, low_alt_Cd, ...
high_alt_Cd, mach_speed, specific_alt)
low_alt_interpolation = interp1(mach_num, low_alt_Cd, mach_speed);
high_alt_interpolation = interp1(mach_num, high_alt_Cd, mach_speed);
Cd_specific_alt = interp1([low_alt high_alt], [low_alt_interpolation high_alt_interpolation], specific_alt);
end

12
Code/DragInterpolations_Spark.m
View File

@ -0,0 +1,12 @@
function [coeffdrag] = DragInterpolations_Spark(velocity, altitude)
spark_drag_data = dlmread("Spark/Spark_DragData.txt", "\t",1,0);
[temp, sound_speed, P, rho] = atmosisa(altitude);
mach_speed = velocity/sound_speed;
mach = spark_drag_data(:,1);
drag = spark_drag_data(:,2);
coeffdrag = interp1(mach, drag, mach_speed);
end

6
Code/ICANSAT/ICANSAT_DragData.txt
View File

@ -0,0 +1,6 @@
Mach NBD FBD
0.00 0.259 0.35
0.30 0.26 0.35
0.60 0.229 0.324
0.90 0.233 0.328
0.99 0.261 0.4

64
Code/ICANSAT/V0/ICANSAT.m
View File

@ -0,0 +1,64 @@
clc
addpath('../')
addpath('../..')
% Df, Lr, Ls, ms, mp, tb
icansat = SingleStageRocket(0.1567, 2.391, 4, 14, 1.92, 8.0);
% Timestep, no. of timesteps, launch angle
icansat.SetSimulationParams(0.01, 10000, 88);
icansat.SetRocketID("ICANSAT");
icansat.InitializeVars();
% Uncomment to use non-constant mass flow rate
% icansat.SetVariableMassFlow();
icansat.CalculateThrust('ICANSAT_Thrust_v_Time.txt');
icansat.CalculateTrajectory();
icansat.WriteDataFile();
subplot(2,3,1);
plot(icansat.t(1:icansat.NStepEnd), icansat.y(1:icansat.NStepEnd)/1000);
hold on
xlabel('Time (s)');
ylabel('Altitude (km)');
ylim([0 2]);
grid on
subplot(2,3,2);
plot(icansat.t(1:icansat.NStepEnd), icansat.M(1:icansat.NStepEnd));
hold on
xlabel('Time (s)');
ylabel('Mach No.');
grid on
subplot(2,3,3);
plot(icansat.t(1:icansat.NStepEnd), icansat.thrust(1:icansat.NStepEnd));
hold on
xlabel('Time (s)');
ylabel('Thrust (N)');
grid on
subplot(2,3,4);
plot(icansat.x(1:icansat.NStepEnd)/1000, icansat.y(1:icansat.NStepEnd)/1000);
hold on
xlabel('Downrange(km)');
ylabel('Altitude (km)');
ylim([0 2]);
grid on
subplot(2,3,5);
plot(icansat.t(1:icansat.NStepEnd), icansat.DvDt(1:icansat.NStepEnd));
hold on
plot([0,50], [0, 0], 'k');
xlabel('Time (s)');
ylabel('Accleration (m/s^2)');
grid on
subplot(2,3,6);
plot(icansat.t(1:icansat.NStepEnd), icansat.m(1:icansat.NStepEnd));
hold on
xlabel('Time (s)');
ylabel('Mass (kg)');
grid on
% ---------

202
Code/ICANSAT/V0/ICANSAT_Thrust_v_Time.txt
View File

@ -0,0 +1,202 @@
0.00 0.000
0.01 857.647
0.04 954.785
0.08 913.544
0.12 927.823
0.16 919.100
0.20 919.595
0.24 916.682
0.28 915.395
0.32 913.525
0.36 912.058
0.40 910.551
0.44 909.180
0.48 907.763
0.52 906.496
0.56 905.294
0.60 904.034
0.64 902.880
0.68 901.733
0.72 900.644
0.76 899.489
0.80 898.448
0.84 897.394
0.88 896.390
0.92 895.312
0.96 894.339
1.00 893.343
1.04 892.377
1.08 891.435
1.12 890.419
1.16 889.500
1.20 888.499
1.24 887.654
1.28 886.720
1.32 885.754
1.36 884.862
1.40 883.893
1.44 883.073
1.48 882.110
1.52 881.241
1.56 880.337
1.60 879.399
1.64 878.531
1.68 877.641
1.72 876.782
1.76 875.889
1.80 874.959
1.84 874.098
1.88 873.210
1.92 872.281
1.96 871.422
2.00 870.539
2.04 869.685
2.08 868.795
2.12 867.867
2.16 867.006
2.20 866.117
2.24 865.186
2.28 864.323
2.32 863.430
2.36 862.496
2.40 861.630
2.44 860.789
2.48 859.850
2.52 858.979
2.56 858.077
2.60 857.133
2.64 856.256
2.68 855.348
2.72 854.450
2.76 853.492
2.80 852.609
2.84 851.690
2.88 850.731
2.92 849.837
2.96 848.912
3.00 847.996
3.04 397.662
3.08 270.454
3.12 262.949
3.16 254.969
3.20 247.296
3.24 239.841
3.28 232.600
3.32 225.594
3.36 218.788
3.40 212.177
3.44 205.794
3.48 199.581
3.52 193.569
3.56 187.741
3.60 182.078
3.64 176.589
3.68 171.267
3.72 166.117
3.76 161.104
3.80 156.244
3.84 151.542
3.88 146.973
3.92 142.552
3.96 138.249
4.00 134.093
4.04 130.048
4.08 126.133
4.12 122.330
4.16 118.650
4.20 115.075
4.24 111.608
4.28 108.254
4.32 104.988
4.36 101.828
4.40 98.758
4.44 95.787
4.48 92.895
4.52 90.102
4.56 87.389
4.60 84.752
4.64 82.201
4.68 79.722
4.72 77.324
4.76 74.998
4.80 72.734
4.84 70.543
4.88 68.418
4.92 66.359
4.96 64.357
5.00 62.417
5.04 60.539
5.08 58.714
5.12 56.945
5.16 55.226
5.20 53.563
5.24 51.947
5.28 50.383
5.32 48.864
5.36 47.391
5.40 45.962
5.44 44.577
5.48 43.231
5.52 41.929
5.56 40.664
5.60 39.437
5.64 38.247
5.68 37.092
5.72 35.973
5.76 34.887
5.80 33.836
5.84 32.814
5.88 31.826
5.92 30.863
5.96 29.930
6.00 29.027
6.04 28.150
6.08 27.302
6.12 26.477
6.16 25.676
6.20 24.900
6.24 24.147
6.28 23.418
6.32 22.710
6.36 22.022
6.40 21.358
6.44 20.712
6.48 20.086
6.52 19.477
6.56 18.889
6.60 18.316
6.64 17.762
6.68 17.224
6.72 16.703
6.76 16.197
6.80 15.706
6.84 15.230
6.88 14.769
6.92 14.320
6.96 13.887
7.00 13.466
7.04 13.058
7.08 12.661
7.12 12.277
7.16 11.904
7.20 11.543
7.24 11.193
7.28 10.852
7.32 10.523
7.36 10.203
7.40 9.890
7.44 9.587
7.48 9.301
7.52 9.020
7.56 8.748
7.60 8.486
7.64 8.213
7.68 7.969
7.72 7.736
7.76 7.509
7.80 7.289
7.84 7.074
7.88 6.864
7.92 6.661
7.96 6.464
8.00 0.000

31
Code/ICANSAT/V1/ICANSAT_V1.m
View File

@ -0,0 +1,31 @@
clc
addpath('../')
addpath('../..')
% Df, Lr, Ls, ms, mp, tb
icansat = SingleStageRocket(0.1567, 2.391, 4, 14, 1.160534, 4.0);
% RunTrajectorySim(icansat, 88);
% RunTrajectorySim(icansat, 86);
% RunTrajectorySim(icansat, 84);
% RunTrajectorySim(icansat, 82);
% RunTrajectorySim(icansat, 80);
% RunTrajectorySim(icansat, 78);
% RunTrajectorySim(icansat, 76);
% RunTrajectorySim(icansat, 74);
% RunTrajectorySim(icansat, 72);
RunTrajectorySim(icansat, 70);
function RunTrajectorySim(rocket, launch_angle)
% Timestep, no. of timesteps, launch angle
rocket.SetSimulationParams(0.01, 10000, launch_angle);
rocket.SetRocketID("ICANSAT");
rocket.InitializeVars();
% Uncomment to use non-constant mass flow rate
rocket.SetVariableMassFlow();
rocket.CalculateThrust('ICANSAT_V1_Thrust_v_Time.txt');
rocket.CalculateTrajectory();
rocket.WriteDataFile();
end

82
Code/ICANSAT/V1/ICANSAT_V1_Thrust_v_Time.txt
View File

@ -0,0 +1,82 @@
0.000000000 28.61809561
0.049382716 865.7172252
0.098765432 828.788367
0.148148148 868.0198335
0.197530864 855.1079959
0.24691358 857.8965696
0.296296296 854.5947699
0.345679012 854.7495807
0.395061728 853.0271145
0.444444444 851.8792688
0.49382716 850.5777483
0.543209877 849.3378201
0.592592593 848.0391454
0.641975309 846.8108491
0.691358025 845.5911444
0.740740741 844.3810837
0.790123457 843.1795528
0.839506173 841.9861293
0.888888889 790.6663048
0.938271605 814.7418358
0.987654321 803.4557486
1.037037037 808.316949
1.086419753 805.7695017
1.135802469 806.6979621
1.185185185 805.9333018
1.234567901 806.0156238
1.283950617 805.7302387
1.333333333 805.6325586
1.382716049 805.403352
1.432098765 805.2754677
1.481481481 805.1345303
1.530864198 804.9574627
1.580246914 804.8347208
1.62962963 804.7205219
1.679012346 804.5587272
1.728395062 804.506189
1.777777778 804.3839195
1.827160494 804.2516444
1.87654321 804.1588364
1.925925926 804.0801232
1.975308642 803.9493204
2.024691358 803.8779173
2.074074074 803.7574503
2.12345679 803.6943463
2.172839506 803.6351266
2.222222222 803.5275729
2.271604938 803.5265503
2.320987654 803.4009857
2.37037037 803.3696582
2.419753086 803.3223497
2.469135802 803.2353159
2.518518519 803.1992204
2.567901235 803.1158937
2.617283951 803.0865164
2.666666667 376.035258
2.716049383 253.688266
2.765432099 248.46047
2.814814815 242.6136227
2.864197531 236.9364371
2.913580247 231.3951724
2.962962963 225.9851368
3.012345679 220.6900131
3.061728395 215.53518
3.111111111 210.4890237
3.160493827 205.5629687
3.209876543 200.7540814
3.259259259 196.0723345
3.308641975 191.4770118
3.358024691 187.002752
3.407407407 182.6232343
3.456790123 178.3595068
3.50617284 174.1854098
3.555555556 170.1104297
3.604938272 166.1321595
3.654320988 162.2588809
3.703703704 158.4668032
3.75308642 154.7548871
3.802469136 151.1405636
3.851851852 147.6122193
3.901234568 144.157891
3.950617284 140.7854521
4.000000000 0.0

43
Code/MAD/MAD v 4.2B_Thrust_v_Time.txt
View File

@ -0,0 +1,43 @@
0.00 0.00
0.10 9359.85
0.50 9263.03
1.00 9184.49
1.50 9118.59
2.00 9061.93
2.50 9012.31
3.00 8968.24
3.50 8928.65
4.00 8892.74
4.50 8859.92
5.00 8829.72
5.50 8801.76
6.00 8775.76
6.50 8751.47
7.00 8728.69
7.50 8707.25
8.00 8687.01
8.50 8667.85
9.00 8649.67
9.50 8632.37
10.00 8615.88
10.50 8600.13
11.00 8585.06
11.50 8570.61
12.00 8556.75
12.50 8543.42
13.00 8530.59
13.50 8518.23
14.00 8506.30
14.50 8494.77
15.00 8483.63
15.50 8472.85
16.00 8462.40
16.50 8452.27
17.00 8442.44
17.50 8432.90
18.00 8423.62
18.50 8414.60
19.00 8405.82
19.50 8397.27
20.00 8388.94
20.10 0.00

11
Code/MAD/MAD_DragData.txt
View File

@ -0,0 +1,11 @@
Mach NBD_0k NBD_10k FBD_0k FBD_10k FBD_30k FBD_45k FBD_60k
0.0 0.38 0.428 0.42 0.468 0.725 1.133 1.819
0.5 0.38 0.428 0.42 0.468 0.725 1.133 1.819
0.8 0.37 0.41 0.4 0.441 0.67 1.026 1.611
1.1 0.51 0.547 0.58 0.617 0.813 1.115 1.607
1.2 0.5 0.54 0.573 0.609 0.793 1.078 1.54
1.5 0.46 0.5 0.53 0.563 0.734 1 1.424
2.0 0.38 0.42 0.438 0.467 0.619 0.854 1.239
3.0 0.28 0.31 0.315 0.34 0.463 0.661 1
4.0 0.217 0.235 0.242 0.261 0.367 0.54 0.842
5.0 0.217 0.235 0.242 0.261 0.367 0.54 0.842

2695
Code/MAD/MAD_Nominal_Inputs.txt
File diff suppressed because it is too large
View File

2603
Code/MAD/MAD_Trajectory_Output.txt
File diff suppressed because it is too large
View File

52
Code/MAD/MAD_Validation.m
View File

@ -0,0 +1,52 @@
clc
addpath('../')
mad = SingleStageRocket(0.205, 4.96, 10, 68, 78, 20.1);
mad.SetSimulationParams(0.05, 10000, 85);
mad.SetRocketID("MAD");
mad.InitializeVars();
mad.SetThrustCorrection();
mad.CalculateThrust('MAD v 4.2B_Thrust_v_Time.txt');
mad.CalculateTrajectory();
mad.WriteDataFile();
% MAD Data
% -----------------------------------------------------------------------------
nominal_input = dlmread("MAD/MAD_Nominal_Inputs.txt", "\t",1,0);
nominal_time = nominal_input(:, 1);
astos_thrust = nominal_input(:, 2);
trajectory_output = dlmread("MAD/MAD_Trajectory_Output.txt", "\t",1,0);
trajectory_time = trajectory_output(:, 1);
altitude = trajectory_output(:, 2);
mach = trajectory_output(:, 6);
subplot(1,3,1);
plot(mad.t(1:mad.NStepEnd), mad.y(1:mad.NStepEnd)/1000);
hold on
plot(trajectory_time, altitude, '--o','MarkerIndices',1:100:length(trajectory_time));
xlabel('Time (s)');
ylabel('Altitude (km)');
legend('2-DoF', 'ASTOS 6-DoF');
ylim([0 60]);
grid on
subplot(1,3,2);
plot(mad.t(1:mad.NStepEnd), mad.M(1:mad.NStepEnd));
hold on
plot(trajectory_time, mach, '--o','MarkerIndices',1:100:length(trajectory_time));
xlabel('Time (s)');
ylabel('Mach No.');
legend('2-DoF', 'ASTOS 6-DoF');
grid on
subplot(1,3,3);
plot(mad.t(1:mad.NStepEnd), mad.thrust(1:mad.NStepEnd)/1000);
hold on
plot(nominal_time, astos_thrust, '--o','MarkerIndices',1:100:length(nominal_time));
xlabel('Time (s)');
ylabel('Corrected Thrust (kN)');
legend('2-DoF', 'ASTOS 6-DoF');
xlim([0 25]);
grid on
% -----------------------------------------------------------------------------

441
Code/SingleStageRocket.m
View File

@ -0,0 +1,441 @@
classdef SingleStageRocket < handle
properties
% Rocket parameters
Df % Diameter of fuselage (meters)
Ls % Length of structure (meters)
Lr % Length of launch rod (meters)
ms % Mass of structure (kg)
mp % Mass of propellant (kg)
tb % burn time (sec)
% Boolean - constant mass flow or variable mass flow
mdot_constant
% Boolean - thrust correction
thrust_correction
% Derived parameters
mt % Total mass
mdot % Mass flow
% Simulation parameters
Timestep
NSteps % No. of timesteps
n_tb % Burnout timestep
% Thrust data file, thrust, isp
thrust_dat_file
thrust
isp
% Simulation variables
t % Time
V, Vcr, M % Velocity magnitude, critical velocity, mach number;
v, u, ux % Velocity components in y and x, x-component velocity calc.4
utheta % X-component of velocity for launch rod
m, Cd % Mass, drag coefficient
x, y, theta % X and Y positions, angle w.r.t. horizontal
DvDt, dvdt % Absolute and relative accelerations
kd, kt % Coefficients of drag and thrust
rho, g, c, T % Density, gravity, sound speed, temperature
area % Cross sectional area
q % Dynamic pressure
launch_angle
% Rocket identifier (MAD, ICANSAT, Spark etc.)
% For purposes of drag identification
ID
% End of simulation
NStepEnd
end
methods
% Initialization
function self = SingleStageRocket(Df, Ls, Lr, ms, mp, tb)
% Set rocket parameters
self.Df = Df;
self.Ls = Ls;
self.Lr = Lr;
self.ms = ms;
self.mp = mp;
self.tb = tb;
% Set booleans
self.mdot_constant = true;
self.thrust_correction = false;
% Derived parameters
self.mt = self.ms + self.mp;
end
% Set rocket ID
function SetRocketID(obj, id)
obj.ID = id;
end
% Set booleans
function SetThrustCorrection(obj)
obj.thrust_correction = true;
end
function SetVariableMassFlow(obj)
obj.mdot_constant = false;
end
% Set simulation parameters
function SetSimulationParams(obj, dt, nsteps, launch_angle)
obj.Timestep = dt;
obj.NSteps = nsteps;
obj.launch_angle = launch_angle;
fprintf('Launch angle: %4.2f degrees.\n', launch_angle);
fprintf('Timestep: %5.3f s. Maximum no. of timesteps: %d.\n', dt, nsteps);
end
% Initialize variables
function InitializeVars(obj)
obj.t = zeros(obj.NSteps, 1);
obj.V = zeros(obj.NSteps, 1);
obj.Vcr = zeros(obj.NSteps, 1);
obj.M = zeros(obj.NSteps, 1);
obj.v = zeros(obj.NSteps, 1);
obj.u = zeros(obj.NSteps, 1);
obj.ux = zeros(obj.NSteps, 1);
obj.utheta = zeros(obj.NSteps, 1);
obj.m = zeros(obj.NSteps, 1);
obj.Cd = zeros(obj.NSteps, 1);
obj.x = zeros(obj.NSteps, 1);
obj.y = zeros(obj.NSteps, 1);
obj.theta = zeros(obj.NSteps, 1);
obj.DvDt = zeros(obj.NSteps, 1);
obj.dvdt = zeros(obj.NSteps, 1);
obj.mdot = zeros(obj.NSteps, 1);
obj.thrust = zeros(obj.NSteps, 1);
obj.isp = zeros(obj.NSteps, 1);
obj.kd = zeros(obj.NSteps, 1);
obj.kt = zeros(obj.NSteps, 1);
obj.rho = zeros(obj.NSteps, 1);
obj.g = zeros(obj.NSteps, 1);
obj.c = zeros(obj.NSteps, 1);
obj.T = zeros(obj.NSteps, 1);
obj.q = zeros(obj.NSteps, 1);
obj.area = (pi/4)*(obj.Df.^2);
% Set initial values
obj.V(1) = 0.447; % Velocity - m/s
obj.theta(1) = deg2rad(obj.launch_angle);
obj.u(1) = obj.V(1)*cos(obj.theta(1)); % x-component of velocity
obj.v(1) = obj.V(1)*sin(obj.theta(1)); % y-component of velocity
obj.y(1) = 0.0; % Altitude - m (XY)
obj.t(1) = 0; % Flight Time - s
obj.M(1) = 0; % Mach Number
obj.c(1) = 340.3961; % Speed of Sound - m/s^2
obj.rho(1) = 1.23; % Sea Level Density - kg/m^3
obj.T(1) = 288.16; % Temperature (K)
obj.g(1) = 9.8066;
obj.m(1) = obj.mt;
end
% Load thrust data (include thrust correction)
function CalculateThrust(obj, file)
data = dlmread(file);
force = data(:,2);
time = data(:,1);
dt = obj.Timestep;
resampled_time = 0:dt:obj.tb;
resampled_force = interp1(time,force,resampled_time,'linear');
for i = 1:length(resampled_time)
obj.thrust(i) = resampled_force(i);
end
obj.n_tb = length(resampled_time);
% Find mass flow rate up till end of burnout
obj.CalculateMassFlow(length(resampled_time));
end
function [t,p,rho] = CalcAtmosQuantities(obj, h)
hg = h; % geometric alt in meters
R = 287; % Gas Constant
re = 6.356766e6; % radius of earth
g0 = 9.8066; % gravity
h = re/(re+hg)*hg;
% Altitude in meters
h0 = 0.0;
h1 = 11000.0;
h2 = 25000.0;
h3 = 47000;
h4 = 53000;
h5 = 79000;
h6 = 90000;
h7 = 105000;
% Temperature in Kelvin.
t0 = 288.16;
t1 = 216.66;
t2 = 216.66;
t3 = 282.66;
t4 = 282.66;
t5 = 165.66;
t6 = 165.66;
t7 = 225.66;
%%equations
a01 = (t1-t0)/(h1-h0); % slope of T/h in K/m for 0 to 11 km
a23 = (t3-t2)/(h3-h2); % slope of T/h in K/m for 25 to 47 km
a45 = (t5-t4)/(h5-h4); % slope of T/h in K/m for 53 to 79 km
a67 = (t7-t6)/(h7-h6); % slope of T/h in K/m for 90 to 105 km
% Pressure in N/m^2
p0 = 101325.0;
p1 = p0*(t1/t0)^(-g0/(a01*R));
p2 = p1*exp(-(g0/(R*t1))*(h2-h1));
p3 = p2*(t3/t2)^(-g0/(a23*R));
p4 = p3*exp(-(g0/(R*t3))*(h4-h3));
p5 = p4*(t5/t4)^(-g0/(a45*R));
p6 = p5*exp(-(g0/(R*t5))*(h6-h5));
p7 = p6*(t7/t6)^(-g0/(a67*R));
if h < h1
t = t0+a01*h;
p = p0*(t/t0)^(-g0/(a01*R));
elseif h < h2
t = t1;
p= p1*exp(-(g0/(R*t))*(h-h1));
elseif h < h3
t = t2+a23*(h-h2);
p = p2*(t/t2)^(-g0/(a23*R));
elseif h < h4
t = t3;
p = p3*exp(-(g0/(R*t))*(h-h3));
elseif h < h5
t = t4+a45*(h-h4);
p = p4*(t/t4)^(-g0/(a45*R));
elseif h < h6
t = t5;
p = p5*exp(-(g0/(R*t))*(h-h5));
elseif h < h7
t = t6+a67*(h-h6);
p = p6*(t/t6)^(-g0/(a67*R));
else
t= t7;
p = p7*exp(-(g0/(R*t))*(h-h7));
end
rho=p/(R*t); % Density in kg/m^3
end
% Thrust correction based on ASTOS manual
function T = ThrustCorrection(obj, h, thrust)
[t,p,rho] = obj.CalcAtmosQuantities(h);
T = thrust*1.095 - 0.014 * p;
end
% Calculate mass at each iteration
function CalculateMassFlow(obj,n_tb)
g0 = 9.8066;
if obj.mdot_constant == true
% Use constant mass flow rate
obj.mdot(2:n_tb) = obj.mp/obj.tb;
% Calculate instantaneous Isp
for i = 2:n_tb
if obj.mdot(i) == 0
obj.isp(i) = 0;
else
obj.isp(i) = obj.thrust(i)/(obj.mdot(i) * g0);
end
end
else
% Assuming constant Isp, use thrust values as ratios over
% the combined sum to obtain the mass flow rate
mdot_sum = obj.mp/obj.Timestep;
force_sum = sum(obj.thrust(1:n_tb));
for i = 2:n_tb
obj.mdot(i) = obj.thrust(i)/force_sum * mdot_sum;
if obj.mdot(i) == 0
obj.isp(i) = 0;
else
obj.isp(i) = obj.thrust(i)/(obj.mdot(i) * g0);
end
end
end
end
% Density calculation
function [T, rho] = CalculateRhoTemp(obj, h)
[T, p, rho] = obj.CalcAtmosQuantities(h);
return;
end
% Gravity calculation
function g = CalculateGravity(obj, h)
g0 = 9.80665; % Sea Level Gravity - m/s^2
RE = 6376000; % Radius of the Earth - m
g = g0*(RE/(RE+h))^2;
return;
end
% Drag interpolation
% Drag interpolation functions need to be specified separately in
% other files due to non-similarity in data.
function [Cd] = DragInterpolation(obj, V, y, t)
switch obj.ID
case "ICANSAT"
Cd = DragInterpolations_ICANSAT(V, y, t, obj.tb);
case "MAD"
Cd = DragInterpolations_MAD(V, y, t, obj.tb);
% Spark data not validated
case "Spark"
Cd = DragInterpolations_Spark(V, y);
end
return;
end
% Trajectory calculation
function CalculateTrajectory(obj)
Rg = 286.9;
nSteps = obj.NSteps;
dt = obj.Timestep;
% Calculate trajectory
for i = 2:nSteps
obj.t(i) = obj.t(i-1) + dt;
V_prev = obj.V(i-1);
y_prev = obj.y(i-1);
t_current = obj.t(i);
obj.Cd(i) = obj.DragInterpolation(V_prev, y_prev, t_current);
% Get thrust from resampled thrust data
if obj.t(i) >= obj.tb
thrust_i = 0;
else
% Thrust correction if boolean true
if obj.thrust_correction == true
obj.thrust(i) = obj.ThrustCorrection(obj.y(i-1), obj.thrust(i));
if obj.thrust(i) < 0
obj.thrust(i) = 0;
end
end
thrust_i = obj.thrust(i);
end
% Update mass
obj.m(i) = obj.m(i-1) - obj.mdot(i-1)*dt;
% Find critical velocity for leaving launch structure
obj.Vcr(i) = sqrt((thrust_i/obj.m(i))*obj.Ls*2);
% X Y Velocity Component Equations
obj.kd(i) = (obj.rho(i-1)*obj.Cd(i)*obj.area)/(2*obj.m(i));
obj.kt(i) = thrust_i/obj.m(i);
obj.V(i) = sqrt(obj.u(i-1)^2 + obj.v(i-1)^2);
obj.DvDt(i) = obj.kt(i)*obj.v(i-1)/obj.V(i) - ...
obj.kd(i)*obj.v(i-1)*obj.V(i) - obj.g(i-1);
% Use non-zero acceleration during burn phase
if obj.DvDt(i) < 0 && obj.t(i-1) <= obj.tb
obj.dvdt(i) = 0;
else
obj.dvdt(i) = obj.DvDt(i);
end
obj.v(i) = obj.v(i-1) + obj.dvdt(i)*dt;
obj.utheta(i) = obj.v(i) * obj.u(i-1)/obj.v(i-1);
obj.ux(i) = obj.u(i-1)+((obj.kt(i)*obj.u(i-1)/obj.V(i-1)) - ...
obj.kd(i)*obj.u(i-1)*obj.V(i-1))*dt;
% Zero acceleration on the ground
if obj.DvDt(i) < 0 && obj.t(i-1) <= obj.tb
obj.u(i) = obj.u(i-1);
% Within the launch structure
elseif obj.x(i-1) <= obj.Lr*cos(obj.theta(1))
obj.u(i) = obj.utheta(i);
% Under critical velocity for leaving launch structure
elseif obj.V(i) < obj.Vcr(i)
obj.u(i) = obj.utheta(i);
% In the air
else
obj.u(i) = obj.ux(i);
end
obj.theta(i) = atan(obj.v(i)/obj.u(i));
obj.x(i) = obj.x(i-1) + obj.u(i)*dt;
obj.y(i) = obj.y(i-1) + obj.v(i)*dt;
% Calculate dynamic pressure
obj.q(i) = 0.5*obj.rho(i-1)*(obj.V(i-1)^2);
% Update density
[obj.T(i), obj.rho(i)] = obj.CalculateRhoTemp(obj.y(i));
% Update sound speed
obj.c(i) = sqrt(1.4*Rg*obj.T(i));
% Update gravity
obj.g(i) = obj.CalculateGravity(obj.y(i));
% Update Mach no.
obj.M(i) = abs(obj.V(i))/obj.c(i);
if obj.y(i) < 0
obj.NStepEnd = i;
break
end
end
obj.theta = obj.theta*(180/pi);
obj.q = obj.q./1000;
%% PREDICTED PERFORMANCE
Vx = max(obj.V);
Mx = max(obj.M);
Yx = max(obj.y);
Xx = max(obj.x);
tx = max(obj.t);
Qx = max(obj.q);
[valx,idxx] = max(obj.x);
[valy,idxy] = max(obj.y);
tX = obj.t(idxx);
tY = obj.t(idxy);
disp(' ');
formatSpec = 'Apogee is %6.3f m occuring at %5.3f seconds.\n';
fprintf(formatSpec,Yx,tY)
formatSpec = 'Landing Distance is %6.3f m occuring at %5.3f seconds.\n\n';
fprintf(formatSpec,Xx,tX)
end
% Write data file
function WriteDataFile(obj)
if obj.mdot_constant == true
result_dat_file = obj.ID + "_Result_" + obj.launch_angle + "deg.txt";
else
text = "VariableMassFlow";
result_dat_file = obj.ID + "_Result_" + obj.launch_angle + "_deg_" + text + ".txt";
end
result_matrix = [obj.t obj.x/1000 obj.y/1000 obj.M obj.thrust obj.DvDt obj.m obj.mdot obj.isp];
result_data = reshape(result_matrix, [], 9);
dlmwrite(result_dat_file, 'Time(s),x(km),y(km),Mach,Thrust(N),Acceleration(m/s2),Mass(kg),mdot(kg/s),Isp (s)','delimiter','')
dlmwrite(result_dat_file, result_data, '-append', 'delimiter', ',');
end
end
end

8
Code/Spark/Spark_DragData.txt
View File

@ -0,0 +1,8 @@
Mach CA
0 0.41
0.1 0.41
0.5 0.347
0.8 0.328
0.95 0.347
1.1 0.441
1.2 0.434

12
README.md
View File

@ -1,2 +1,12 @@
# 2-DoF-TrajectorySim
## 2-DoF Rocket Trajectory Simulation Program
This 2-DoF rocket trajectory simulation is reproduced from the Master's thesis by Zachary Doucet. The original code has been rewritten in a OOP form so that it can be reused in a simple manner for future trajectories of single-stage rockets. Drag interpolations are included as standalone files for each rocket since data format is non-similar in nature.
A validation test has been made and compared against the 6-DoF software ASTOS, from which the results of the MAD rocket have been obtained. Figures shown below are the comparisons of the altitude, Mach number and the corrected thrust curve for both 2-DoF and 6-DoF simulations.
![](https://git.equatorial.space/jason.ong/2-DoF-TrajectorySim/raw/branch/master/2-DoF%20vs%20Astos%206-Dof.png)
Input settings used for the 2-DoF simulation for validation are found in the MAD_Validation.m file. A similar workflow has now been created for the ICANSAT data to demonstrate code-reuseability.
Comments:
1. Numerical calculations are adapted based on the Angle_Tests_XY.m program used in the original code by Zachery Doucet.
2. Launch tower deflection is not included.

BIN
zjd11_Thesis_Paper_Final.docx
View File

Write Preview
Loading…
Cancel
Save