Solving system of motion equations. Need help

[Phil Smith]

Hey guys!
I try to figure out how to make ballistic calculations for a real launch vehicle (a good ol' backyard rocketry project with dreams of put 15 kilos in LEO :lol: ).
I've integrated basic equations of motion (just in 3 degree of freedom system for now) and have a problems to solve that system of non-leaner equations.
Here what I got:
a) The initial system of deferential equations (taken from some Russian book called "Design of launch vehicle and satellite trajectories" - see Fig.01 in attachments.
b) The same system after solving diff/eqs - see Fig.02,
where:
v = vehicle velocity in certain time;
u'_v - engine's exhaust velocity in vacuum, m/s;
u'_0 - engine's exhaust velocity at sea level, m/s;
p(h) - known function atmosphere pressure vs altitude, Pa;
p(0) - atmosphere pressure at launch site, Pa;
μ - vehicle weight ratio = (vehicle weight at certain time)/(total vehicle start mass) or μ = (m_0 - mdot*t)/m_0 - basically it's our iteration step;
m_0 - vehicle start mass, kg;
mdot - engine's flow rate, kg/sec;
S - Vehicle max cress-section area, m2;
c_x(M) - known function of vehicle drag coef. vs Mach number;
ρ(h) - known function of air density vs altitude, kg/m3;
g(h) - known function of g-force vs altitude, m/sec^2;
θ - velocity angle, rad;
r - vehicle radius-vector 1;
φ_pr(t) - known pitch angle program, rad;
l_1 - vehicle length, m;
x_c.p.(α) - known function of vehicle center of pressure vs angle of attack, m;
x_c.m.(t) - known function of vehicle center of gravity vs time, m;
с'_y1 - known function of vehicle lift coef. vs Mach number;
x & y - vehicle coordinates from launch site, m;
M - mach number;
a(h) - known function of speed of sound vs altitude, m/s;
r' - radius-vector 2, m;
R - mean Earth radius, m;
h - altitude, m;
α - angle of attack, rad.

oh yeah, that's it..probably..)
So my question is - what is the best way of compute these equations? (via excel for example) here are 9 unknowns and 9 equations - just like in a text book, but what makes me confused is known functions vs unknown values ( p(h) for ex. - we got a table with altitude vs pressure, density, speed of sound, temperature, etc, but we dont know altitude itself yet) - it means auto interpolations among these tables and a lot of coding..
So what is your suggestion?
Very appreciate any help.
PS - my mind is a mess right know:lol:

 

[Thunder Chicken]

Look up Euler's method - it is an explicit method that is pretty intuitive and easy to enter inta spreadsheet. Only disadvantage is that it is only first-order accurate, but if you make your time step small enough it won't matter.

At a given time t with initial position x(t) and velocity v(t) you can calculate all of the forces F(t) on the vessel at that time. Then you basically just solve a linearized version of F = ma to update to a new position and velocity assuming F is approximately constant at that time step.

a(t) = F(t)/m(t)

v(t+dt) = v(t) + a(t)*dt

x(t+dt) = x(t) + v(t)*dt + (1/2)*a(t)*dt*dt

Actually, table interpolations are relatively efficient vs evaluating a formula. You can get Excel to do interpolations.

There are more accurate and complicated methods, but this is a good place to start and will get you where you want to go.

---------- Post added at 07:32 PM ---------- Previous post was at 07:18 PM ----------

I actually just saw your attachments (didn't show up before). It looks like everything indexes on mu at some level, so you should be able to calculated h, p(h), x, y, etc. The velocity equations looks iterative. Hmm.

 

[Phil Smith]

thanks for your support!
a(t) = F(t)/m(t) - it's a basic. add gravity loss and we get the same first eq. on Figure 01 at my attachments -
a(t) or dV/dt = (F(t)/m(t)) - g*cosθ - g*(x/r)*sinθ, so
V - V0 = ((F(t)/m(t)) - g*cosθ - g*(x/r)*sinθ)(t-t0)
or
delta_V = ((F(t)/m(t)) - g*sinθ - g*(x/r)*cosθ)*delta_t,
where delta_t is out step-size.
F(t) = F*cosδ - X = F*cosδ - c_x*q*S
where F - engine(s) thrust, δ - engine(s) gimbal angle, X - drag force.
but dear Lord, it smells like a numerical analysis here. I read some books, and found a Runge–Kutta method is used in most ballistic calculations:
https://en.wikipedia.org/wiki/Runge–Kutta_methods
I try to dig in that direction..
PS. I don't know why, but in all examples I read they integrate this system for weight ratio μ, rather than t.

 

[Thunder Chicken]

Runge-Kutta is a higher order numerical method that is appropriate for this application. Euler is simpler and perhaps easier to get into a spreadsheet, but not quite as accurate.

The mass ratio is a proxy for time t as both appear to be linear (constant thrust?).

There is a pitch program - is that provided as a function of time or altitude, or something else?

I *think* you need to implement an iterative routine to solve this system.

Assuming initial mass, time, and velocity, solve for an initial "guess" solution, say assuming no pitch program or drag - rocket flies straight upward. You can check your numerical solution against the [ame="http://en.wikipedia.org/wiki/Tsiolkovsky_rocket_equation"]Tsiolkovsky equation[/ame].

At this point you'll have a preliminary solution for velocity and altitude verified against an analytical equation (always a good starting point when working with numerical analysis). With these initial guess values you can then evaluate the altitude and velocity dependent parameters like density and drag and then include those in the equations for the next iteration. Update the position and velocity solution, then update the altitude dependent stuff, Lather, rinse and repeat the process until the solution converges. Depending on how the pitch program is provided, you can roll that into the math and include that.

Getting this into a spreadsheet will be challenging but possible. However, this is really something better done in a program that can handle loops, especially if you have a lot of table lookups.

Some of the folks more versed in Visual C++ stuff in Windows might be able to steer you to how to program this. Matlab is a good option if you have access to it. I work on Linux and use a program called Octave that is very similar to Matlab for this sort of calculation.

EDIT: Look into 'numerical shooting methods'.

 

[Thunder Chicken]

I pulled up Octave and started playing with a shooting method using parameters for a V2 rocket fired directly upward.

This code will horrify the professional coders in the audience with its inefficiency but hopefully it illustrates the shooting method. You can work on optimization once you understand what it is doing. This should run in Matlab if you have access to it. It should also be reasonably easy to convert to another programming language.

clear
clc
%
% Numerical Rocket Trajectory Calculator
% Specs based on V-2 rocket http://www.astronautix.com/lvs/v2.htm
%
% Rocket empty mass (kg)
mass_empty = 4008;
% Rocket fuel mass on pad (kg)
mass_fuel = 8797;
% Rocket diameter (m)
d = 1.65; 
% Frontal area of rocket (m2)
A = 0.25*pi*d*d;
% Vacuum Thrust (N)
thrust_vac = 311800;
% ISP at sea level (s)
ISP_sl = 203;
% ISP in vacuum (s)
ISP_va = 239;
% Burn time (s)
t_burn = 68;
%
% Physical constants
%
% gravitational acceleration at Earth's surface (m/s^2)
g_0 = 9.80665;
% standard atmospheric pressure at sea level (N/m^2)
p_0 = 101325;
% standard atmospheric temperature at sea level (K)
T_0 = 288.15;
% lapse rate (K/m)
L = 0.0065;
% molar mass of dry air (kg/mol)
M = 0.0289644; 
% universal gas constant (J/(mol•K))
R = 8.31447;
%
% Numerical parameters
%
% Number of time steps per shot
nsteps = 500;
% Time step size (s) %smaller step size ~ 0.1 is a bit more accurate but memory costs go up by 10X
dt = 1.0;
% Number of shooting iterations
nshots = 100;

% Allocate memory

a = zeros(nsteps);
v = zeros(nsteps);
x = zeros(nsteps);
t = zeros(nsteps);

vts = zeros(nsteps); %check against the Tsiolkovsky Eqn.
xts = zeros(nsteps); %check against the Tsiolkovsky Eqn.

% populate time vectors here (doesn't change shot to shot)

for i = 2:nsteps % t(1) = 0, already initialized.

	t(i) = t(i-1) + dt;

end

% Perform initial "shot" assuming no drag, implicit Euler 

for i = 2:nsteps % i = 1 is the initial condition, doesn't change

	%check to see if motor is still burning
	
	if t(i) <= t_burn
		mass_flow = mass_fuel/t_burn;
	else
		mass_flow = 0;
	end

	% calculate average mass and weight at middle of time step
	
	mass = mass_empty + mass_fuel*max(0,(t_burn - (i-0.5)*dt))/t_burn;
	weight = mass * g_0;
	
	% no pressure correction for ISP with altitude for first shot.
	
	ISP = ISP_va;
	
	% calculate thrust
	thrust = ISP * mass_flow * g_0;
	
	% sum forces on rocket (going straight up, no drag)
	force = thrust - weight;

	%update kinematics
	
	a(i) = force / mass;
	v(i) = v(i-1) + a(i)*dt;
	x(i) = x(i-1) + v(i-1)*dt + 0.5*a(i)*dt*dt;
	
	%check velocity against Tsiolkovsky Eqn.
	%won't agree if drag and non-vacuum ISP are included or if burn is terminated before nsteps is reached.
	
	vts(i) = ISP_va*g_0*log((mass_empty + mass_fuel)/mass) - g_0*t(i);
	xts(i) = x(i);

end

% OK, now we have a preliminary altitiude and velocity "guess" from the first "shot".
% Now we can calculate all the velocity and altitude- related forces explicitly using the guessed data
% and repeatedly update the solution for a number of shots.

for j = 1:nshots

for i = 2:nsteps % i = 1 is the initial condition, doesn't change

	%check to see if motor is still burning
	
	if t(i) <= t_burn
		mass_flow = mass_fuel/t_burn;
	else
		mass_flow = 0;
	end

	% calculate average mass and weight at middle of time step
	
	mass = mass_empty + mass_fuel*max(0,(t_burn - (i-0.5)*dt))/t_burn;
	weight = mass * g_0;
	
	% apply pressure correction for ISP with average altitude at middle of time step.
	
	ISP = ISP_va - ((1-(L*0.5*(x(i-1)+x(i))/T_0))^((g_0*M)/(R*L)))*(ISP_va - ISP_sl);
	
	% calculate thrust 
	% Check if routine works by setting ISP = ISP_va (should get Tsiolkovsky Eqn)
	
	thrust = ISP * mass_flow * g_0;
	
	% calculate drag
	% This is a rough hack, really should get a drag curve as a function of Mach number
	% I'm just calling the drag coefficient Cd = 0.01 for simplicity.
	% I am also just using air density rho at sea level, really should plug in a correction for altitude.
	% I'm calculating drag at the average velocity between the start and end of the time step.
	
	Cd = 0.01; 	% drag coefficient (dimensionless). 
			% Check if routine works by setting Cd = 0 (should get Tsiolkovsky Eqn)
	
	rho = 1.223; 	%air density at sea level (kg/m3)
	
	drag = Cd*0.5*rho*A*0.25*(v(i)+v(i-1))*(v(i)+v(i-1));
	
	% sum forces on rocket (going straight up, drag points down)
	
	force = thrust - weight - drag;

	%update kinematics
	
	a(i) = force / mass;
	v(i) = v(i-1) + a(i)*dt;
	x(i) = x(i-1) + v(i-1)*dt + 0.5*a(i)*dt*dt;
	
end

end

plot(t, xts, t, x)
xlabel ('Flight Time (s)')
ylabel ('Altitude (m)')

I've attached the plot output from this program. Upper blue curve is the solution with no drag. The lower curve is the one with the crude drag model applied. The engine thrusts for 68 seconds, afterwards the upper curve is a parabolic arc.

You should also build in some sort of convergence check to make sure the trajectory stops changing shot to shot. I simply told the program to shoot 100 times, updating the solution each time.

Hope this helps. It's been fun tinkering with this stuff again. Cheers

 

[Phil Smith]

Damn, you're a man!
I'll try this code soon! Thanks a lot!
PS right now i've been creating math model of Earth in matlab (pressure, density, temperature, speed of sound, g vs altitude) based on "US standard atmosphere" & Russian GOST 4401-81 "standard atmosphere". This matrices ARE huge (for range of 0-250 km there are 1173 values of each parameter..). I'll upload it when I've done. Hope it'll be useful.

 

[Phil Smith]

So I've been started digging into the system and stuck..))
MatLab doesn't show any results - all parameters are still zeros..
As far as the vehicle has pitch program, now there are two coordinates - x and y and also angle of velocity THETA.
The system looks like this:

a(i) = (Thrust - Drag)/m(t) - g_i*sin(theta(i)) - (x(i)/r)*g_i*cos(theta(i));
theta(i) = theta(i-1)- (dt/v(i))*(((fi_pr_i - theta(i))/cur_mass)*(Thrust + Lift_1) - g_i*(cos(theta(i)) + (x(i)/r)*sin(theta(i))));
v(i) = v(i-1) + a(i)*dt;

% x = x(i-1) + vx(i)*dt + 0.5*ax(i)*dt*dt,
% if vx(i) = v(i)*cos(theta(i))
% and vy(i) = v(i)*sin(theta(i))
% then:

x(i) = x(i-1) + v(i)*cos(theta(i))*dt + 0.5*a(i)*cos(theta(i))*dt*dt = x(i-1) + a(i)*dt*cos(theta(i))*dt + 0.5*a(i)*cos(theta(i))*dt*dt = x(i-1) + 1.5*a(i)*cos(theta(i));
% same for y:
y = y(i-1) + vy(i)*dt + 0.5*ay(i)*dt*dt = y(i-1) + 1.5*a(i)*sin(theta(i));

Also have a problem setting initial theta(i) value = 90 degree (1.57 rad) (cause we launch the rocket vertically)

I've typed interpolation tables for finding pressure/density/Speed of sound vs altitude, lift and drag coefficients vs mach number and pitch angle vs time and placed em in separate file called "All_parameters.m".

I attach .zip archive this two matlab files (main.m and All_parameters.m) - they both need to be located in the same folder, and for calculation you should start main.m. Also there is a word file to see how the system of equations looks like in a proper form))

PS. Little bit about flight plan:
T + 000.00 sec - lift off,
T + 120.00 sec - payload fairing jettison;
T + 185.00 sec - 1st stage engine cut-off;
----free ballistic flight----
T + 295.00 sec - 2nd stage engine ignition;
T + 449.00 sec - 2nd stage engine cut-off, orbit insertion.

so I'll appreciate any help to solve this problem..

 

[Phil Smith]

Gosh, found a lot major issues at previous code..
so finally I've got my system working well.. At least I guess))
this code computes just 1st stage by now.
I've attach current version of the program and screen of velocity and pitch angles of my vehicle first stage.
I а you wanna, you can put your stages performance you know well just to confirm it works or doesn't.
As for pitch angle program - I got short vertical segment (from 0 to 10 sec), a little turn and gravity turn almost to the end of the first stage burn (180 sec) and little correction setting pitch angle to 30 degrees at the stage cut-off (185 sec).

PS Cx(M) and Cy(M) functions are taken from TITAN II booster data.
Atmospheric tables (p(h), rho(h) and speed of sound(h)) are from 0 to 269 km by now, I gotta increase it up to 645 km as soon as I can. Also you can use em if you wanna - I don't mind))
and again, huge thanks to Thunder Chicken - your basic code structure was a great help to solve my problem!

PS PS - I got a big concern - as you see from diagram aerodynamic velocity loss is 707 m/s. For example Saturn V S-Ic stage had 50 m/s aerodynamic loss only, and my vehicle has almost the same shape. Maybe it's because using drag and lift coefficients of bulky Titan II..


oops, I got a problem with x coordinate - it doesn't change almost up to 90 sec, but horizontal motion should occur after 8-10 sec when vehicle starts pitching.. strange..

 

[Thunder Chicken]

You've been busy!

Your velocity loss may well be significantly different between the Titan and Saturn V as the Titan accelerated more quickly than the Saturn (i.e. it was pushing through the high-density parts of the atmosphere faster).

Did you try this with V-2 rocket data? Being a single stage with a suborbital trajectory, it is an easier test case and things like the apogee and downrange distance are known. I also believe there is some NACA/NASA data on captured V-2 drag coefficients somewhere.

 

[Phil Smith]

You've been busy!

Your velocity loss may well be significantly different between the Titan and Saturn V as the Titan accelerated more quickly than the Saturn (i.e. it was pushing through the high-density parts of the atmosphere faster).

Did you try this with V-2 rocket data? Being a single stage with a suborbital trajectory, it is an easier test case and things like the apogee and downrange distance are known. I also believe there is some NACA/NASA data on captured V-2 drag coefficients somewhere.

oh yeah, I have)) Fortunately, my last week was a vacation week))
perhaps you're right, I'll try to dig into V-2 simulations until I got it's working well. recently I saw some of its aerodynamic parameters

 

[Phil Smith]

Hey everybody!)
Oh boy, finally gathered enough info about V-2 missile and optimize calculations, but with Excel - MatLab didn't wanna compute em right desperately so I decided to create an excel spread sheet(s) for this purpose.
So results really look nice. I've attack my .xlsx file if you'd like to check it out - there are 6 sheets:
1. Theory - it's an explanation of what formulas I use for computing the problem;
2. Input data sheet - All input parameters of a vehicle and launch conditions, for the Earth model the WGS-84 ellipsoid is used;
3. Main - actually the main sheet of the file with A-AA columns, stepsize (now it's 0.01 sec) and few diagrams at the top right table corner - velocity vs time, altitude and range vs time, acceleration vs time and pitch angle, velocity vector and attack angles vs time;
4. Pitch program - Interpolation table for finding a current pitch angle;
5, Cd & Cl vs M - Interpolation table for finding current Drag and Lift coefficients versus Mach number;
6. Atmospheric data - Interpolation table for finding current air pressure, density and speed of sound.

Comparing results table (with real v-2 flight data) is also attached.
Btw - altitude and range values are little bit bigger than real V-2 had, bud velocities are pretty much the same

PS I will appreciate much if somebody will test my table with different vehicles. My next step - to simulate mercury-redstone Freedom 7 flight.

PS PS Masses in column "D" should be tweaking manually. I dont know how to create multi-IF command in one cell w/o typing looooong excel-style code.
For example, a problem:
I need to calculate current mass (column D) in each period of time (column A), but I also need to verify engine's still burning and also at 120 sec we got payload fairing jettison. In matlab it would be looking like this (cutoff time= 160 sec, payload fairing mass = M_pf):

if t(i) < 160
    M_current = M(i-1) - mdot * dt
elseif t(i) = 120
    M_current = M(i-1) - mdot * dt - M_pf;
else
    M_current = M_dry;
end

Probably I'll try create some macros for this task.

2010 Excel file:
https://app.box.com/s/b242o0sms6qz5f2oozxu

 

[Phil Smith]

yeah baby! I've adjusted specific impulses little bit and BINGO!
all parameters are within V-2 flight data range.

Fixed version of the excel file:
https://app.box.com/s/b242o0sms6qz5f2oozxu