how to: rocket launch (mass consumption, equation of motion... good literature?)

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!

Answer 1

[Hynod]

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 = ...

Invalid expression. Check for missing or extra characters.

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!

[moinmoinnoob69]

Thanks a lot!! I'm really looking forward to going through your code tomorrow morning :D