how to: rocket launch (mass consumption, equation of motion... good literature?)
12 views (last 30 days)
Show older comments
Hi :)
Sry for the maybe obvious questions, but i'm pretty new to Matlab.
Background:
We should model and simulate a rocketflight into stratosphere and back to earth. Goal is to launch from coordinates (lat,lon,a) to different coordinates on earth.
1.Question:
How is it possible to formulate a time dependent mass equation for a rocket launch in Matlab?
Stats:
The rocket should have 2 stages of 'propellant-burning'.
I've considered:
P.I_sp_1 = 180; % [s], specific impusle stage 1
P.I_sp_2 = 100; % [s], -"- stage 2
P.m_dot = 515; % [kg/s], constant massflow
P.mfuel_1 = P.I_sp_1*P.m_dot; % [kg]
P.mfuel_2 = P.I_sp_2*P.m_dot; % [kg]
P.m_0 = 10000; % [kg] empty rocketmass
P.m_rkt_0 = P.m_0 + P.mfuel_1 + P.mfuel_2; % initial mass
How do I formulate m(t) for the whole flight? (do i have to know the whole 'flighttime'?)
2.Question:
How can i formulate this 'equation of motion' best in Matlab? -> i want to get the velocity (v(t)) & position vector (r(t)) out of it!
F_net= P.F_gravity + P.F_drag + P.F_thrust; % vectors of form [1x3]
% F_net = m(t) * (ddot*r)/(dt)^2 ... equation of motion
% with:
P.F_gravity = gravitywgs84(height,lat,lon,'Exact','None')
% P.F_thrust = dependent on azimuth, elevation, specific impusle, massflow
% P.F_drag = 0.5 * P.rho * P.velocity^2 * P.reference_area * P.drag_coefficient
% -> P.rho is dependent of altitude (ref. NASA)
% -> P.reference_area = const.
% -> P.drag_coefficient = const.
% -> P.velocity = still required! --> wanted to get it out of the mass equation above
Thanks in advance!
If anyone of you has some good helping literature or code to our 'Rocket-Project' I would be really happy!
Thanks!
0 Comments
Answers (1)
Hynod
on 12 Sep 2023
Edited: Hynod
on 12 Sep 2023
Hi @moinmoinnoob69, Hoping you are practic with ordinary differential equations, I think that it is appropriate to approach the script with the two body problem structure with a perturbation ,which in this case is the thrust. I will assume that you have the thrust curve of each of the phases of your project , from which you will compute the acceleration dividing the thrust by the mass, each instant. You are lucky because I am writing a thesis about a spacecraft which orbits around Mercury, and also if it seems very different, the structure of the script and the functions are not:
This is the main script
clc,clear
format long
mu_m = ...
s_0 = ...
%%%%%%%%%% Integrazione del de-orbiting
options = odeset('RelTol', 1.e-16, 'AbsTol', 1.e-16);
[t,S_1] = ode113(@ODE_1_LanderOrbit, t_span1 , s_0' , options) ;
%%%%%%%%%% Integrazione della fase propulsa a solido e coasting %%%%%%%%%%
[t,m_2] = ode113(@ODE_SRMSupport, 0:0.01:(t_f2-t_i2) , 3490 ) ;
t_span_SRM_int = 0:2.5:65 ;
F_SRM_int = 4.44822*[0 17500 18000 20000 22000 23500 26000 25000 26000 26500 27500 28500 30000 31000 31000 32000 32500 34000 34000 34000 34000 32500 32000 32000 31000 15000 0 ]; % La spinta in libbre , bha !
F_SRM_eva = interp1(t_span_SRM_int,F_SRM_int,0:0.01:(t_f2-t_i2),'linear','extrap') ;
for i2 = 1:length(0:0.01:(t_f2-t_i2))
a_SRM(i2)= F_SRM_eva(i2)/m_2(i2) ;
end
[t, S_2] = ode113(@(t, x) ODE_2_LanderSRM( t, x, a_SRM , t_span2 ), t_span2, S_1(end, 1:6));
[t, S_3] = ode113(@ODE_3_LanderCoasting, t_span3, [S_2(end, 1:6) ]);
%%%%%%%%%% Integrazione della fase propulsa a liquido %%%%%%%%%%
[t,m_3] = ode113(@ODE_LRESupport, 0:0.01:(t_f4-t_i4) , m_2(end)-560*0.453592 ) ;
t_span_LRE_int = [0 0.1 0.5 0.7 5 7.5 12 13 13.5 15 18 19 19.001] ;
F_LRE_int = 1.03*4.44822*[0 1400 1450 1400 2100 2400 2800 2825 2850 2850 2850 1000 0 ];;
F_LRE_eva = interp1(t_span_LRE_int,F_LRE_int,0:0.01:(t_f4-t_i4),'linear','extrap');
for i_2 = 1:length(0:0.01:(t_f4-t_i4))
a_LRE(i_2)= F_LRE_eva(i_2)/m_3(i_2);
end
[t, S_4] = ode113(@(t, x) ODE_4_LanderLRE(t, x, a_LRE,t_span4), t_span4, S_3(end, 1:6));
Here it is an example of the function that will you provide the mass each instant
function dA = ODE_SRMSupport(t,A1)
g_0 = 3.7;
m_p = 6205 * 0.453592 ;
I_t = 4.44822 * 1911070 ;
x= 0:2.5:65 ;
T = 4.44822*[0 17500 18000 20000 22000 23500 26000 25000 26000 26500 27500 28500 30000 31000 31000 32000 32500 34000 34000 34000 34000 32500 32000 32000 31000 15000 0 ]; % La spinta in libbre , bha !
i = interp1(x,T,t,'linear','extrap');
A5 = i*m_p/I_t ;
dA = [ -A5 ]' ;
end
And here an example of the SRM stage function, while the orbiting one can be easily provided removing the a_p. (the LRE is obviously very similar to the SRM one)
function dx = ODE_SRMSupport(t, x, ts,t_span_f)
istat = interp1(t_span_f,ts,t,'linear','extrap');
mu_m = ...
grav = - mu_m / (norm(x(1:3),2))^3;
versore = x(4:6)/norm(x(4:6)) ;
a_p = versore * istat / 1000;
A = [ 0,0,0, 1,0,0
0,0,0, 0,1,0
0,0,0, 0,0,1
grav,0,0, 0,0,0
0,grav,0, 0,0,0
0,0,grav, 0,0,0];
dx = A*x+[0 0 0 -a_p']';
Obviously it will not run because i rudimentally modificated it , deleting the lines that are not necessary for your problem, but I think it is a very good starting point (I also didn't delete my parameters, I am sorry XD). Hope it will help! I am sorry if i will be not very active to answer your questions this week but I am new to this forum and I only got registered because I had an issue about the thesis I am writing. Therefore, I will do my best to solve other problems you will propose. I am not familiar with some good literature about what you have to do, I am sorry about that.
I didn't add the perturbation of the drag (because Mercury doesen't have an atmosphere hahaha), but it can be easily computed as
a_drag = - velocity_versor * (module of the drag, function of the position of the rocket, provided by the integrator itself)
I am sorry for my bad language but I am italian :) .
Good luck!
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!