Solving diff equations

[mardlamock]

Hello everyone, i ve been trying to get the max height for my rocket taking drag into account (irl), variations in the angle of the rocket with respect to the airflow and a whole bunch of other things and i ended up with several differential equations. Im not going to explain much of what the equations look like themselves but just trying to get a general idea. The way i ve been solving them so far is by simply using excel spreadsheets where for a small time interval i treat all of the functions of time as a constant and then do an approximation as if the acceleration was constant. This could sort of be expressed like so,

from t0 to t1:a=r(t1)-n(x'(t0))*(m(x(t0)). Then i get the displacement in t1: 1/2*(r(t1)-n(x'(t1))*(m(x(t0)))*t1^2, and velocity: t1*(r(t1)-n(x'(t1))*(m(x(t0))) From t1 to t2 acceleration would be:

r(t2)-n(x'(t1))*m(x(t1)). As you can see, to solve for t=n i am using previous values from t=n-k, now would there be a way of expressing this whole thing with maybe sigma notation and then take the limit as t1-t2 approaches zero and the ammount of small t s i need approaches infinity.

Sorry if you dont understand anything of what im saying, im not all that good at explaining myself and do not know all that much about these topics.

Thank you very much!

[K^2]

If you ignore change of rocket's mass as it burns fuel and assume constant thrust, there is an analytic solution.

Given an equation:

v' = a - kv²

You can find terminal velocity, v_t = Sqrt[a/k] for which v' = 0. Naturally, a = F/m - g, where F is thrust and m is mass of the rocket, which you can also think of as acceleration of the rocket without drag. In that case, velocity changes as the function of time.

v(t) = v_t Tanh(a t/v_t)

Tanh is hyperbolic tangent. Tanh(x) = sinh(x)/cosh(x). If your software does not have hyperbolic functions, sinh(x) = (exp(x) - exp(-x))/2 and cosh(x) = (exp(x) + exp(-x))/2. This equation gives v(0) = 0, starts out accelerating at rate a, and then settles on v = v_t after a while.

Of course, you want h(t). The solution to h' = v from above yields:

h(t) = v_t²/a ln(cosh(a t/v_t))

So all you need to know are a, k, and the amount of time t that the engine burns to figure out altitude and velocity the rocket attains once the engine cuts out. From there on, rocket coasts.

v' = -g - kv²

The parameter k is still exactly the same, but without thrust, the only other force is gravity, which also slows down your rocket. Re-defining v_t as Sqrt(g/k), we arrive at a solution which is very similar to the earlier equation.

v(t) = -v_t Tan(g t/v_t - c)

This is essentially the same equation, but with Tan instead of Tanh and a bit of a twist that without parameter c, v(0) would be 0. And we want v(0) to be whatever velocity the rocket had when engine cut out. So solving for v(0) = v₀ we get the value for c.

c = Tan^-1(v₀/v_t)

This also tells us when the rocket will reach the apex. It will happen the moment gt/v_t = c. So all that's left is figuring out altitude. Again, I'm going to use h' = v, so this will need some corrections in a moment.

h(t) = v_t²/g log(cos(g t/v_t - c))

You should note immediately that at the apex, I get h = 0. That's because h(0) is negative. Just a side effect of solving a differential equation. But if all you want is the altitude to which rocket climbs, then -h(0) is exactly what you are looking for.

Putting it all together, the full equation for the actual maximum altitude the rocket will reach is the following.

H = v_t1²/a ln(cosh(a T/v_t1)) - v_t2²/g log(cos(-c))

Where:

H : Maximum height.

v_t1 = Sqrt(a/k) : Ascent terminal velocity.

a = F/m - g : Initial acceleration of the rocket

g : Acceleration due to gravity, or 9.8m/s²

T : Length of time that engine runs.

F : Average thrust of the engine. (Assumed to be constant.)

v_t2 = Sqrt(g/k) : Coasting terminal velocity.

c = Tan^-1(v₀/v_t2) : Time offset.

v₀ = v_t1 Tanh(a T/v_t1) : Rocket's velocity when engine cuts out.

k : Drag coefficient divided by mass. It's probably easier to estimate coasting terminal velocity and get k from that. Or use drag formulas available for rocketry.

You should be able to estimate all of these parameters and get a somewhat descent estimate of max height.

Of course, if you want to be more precise, you need to take into account the thrust profile and the change in mass. (Which also depends on thrust profile.) That will require numerical integration.

If you want to do the numerical integration, I would strongly recommend forgetting about Excell and learning to use Matlab/Octave. They use mostly the same language, but you can get Octave for free. I can walk you through setting up a numerical integrator for this problem in Octave, if you want.

Oh, and somebody should check my math above. The functional forms are correct, but I might have messed up the coefficients here and there. Or forgot a minus sign. Or something equally silly. Please, let me know if there is a mistake.

[mardlamock]

I know that solution but that is not taking the density of air changing as a function of height which is quite important when it comes to the max height of my rocket, considering it has an impulse of about 800ns.

The general equation i am looking to solve is this: k/(m(t))-r(x´(t))*n(×(t))/m(t)-x´´(t)=0. In my case, m(t) would be the rocket s mass as a function of time, k the force it exerts,

r(x´(t) the velocity s part of the drag equation, and n(x(t)) the density as a function of height .I ve been trying to solve it by using the regular e^rt method used in homogenous diff equations but it seems to get quite complicated. I would be more than pleased if you could help me set the equations up in Octave and solve it. This will be useful because if I want to build the flight computer i ve been thinking of it d have to solve these equations on the fly to estimate apogee and then even more complex equations for the max deflection angle of the fins and such. Thank you very much!

[Jimbimbibble]

It looks like the initial method you described is essentially Euler's Method (am I correct?). If you want to improve your approximation, you can use the Improved Euler Method and iterate a few times per data point until it converges. This would be too tedious to do by hand, but you could write a program in VBA to solve this. In case you don't know, VBA is a programming language used to write macros is Excel. You can access it by going to the toolbar settings and turn on "developer". You can use the spreadsheet to read in data and output results, which is convenient and there are a lot of built-in functions in VBA that make it easy to learn.

[K^2]

It looks like the initial method you described is essentially Euler's Method (am I correct?).

Looks like it. And while implicit Euler is an improvement, there are better, simpler methods. Even though it's not designed to deal with drag, Verlet should work extremely well here.

The general equation i am looking to solve is this: k/(m(t))-r(x´(t))*n(×(t))/m(t)-x´´(t)=0. In my case, m(t) would be the rocket s mass as a function of time, k the force it exerts,

r(x´(t) the velocity s part of the drag equation, and n(x(t)) the density as a function of height .I ve been trying to solve it by using the regular e^rt method used in homogenous diff equations but it seems to get quite complicated. I would be more than pleased if you could help me set the equations up in Octave and solve it.

Here is the basic code using Velocity Verlet as the integration method. rocket.m

To run it in octave, just find create a file with this code somewhere, say, "C:\code\rocket.m". Start up octave, navigate to this directory using command "cd c:\code" and then just run "rocket". (Extension isn't needed. It will look for a .m file.)

It's pretty basic, but you can probably figure out how to modify it to fit your parameters.

[mardlamock]

Thank you very much, i had no idea what i was doing had a name and everything, i just thought of it on my own and ended up with that method. I also thought of doing a more accurate method where i do the integrals of the rocket s acceleration due to its own thrust (f/m(t)) and get a better result. Thanks a LOT for the code, comments made it much easier to understand. My math teacher got me a book on diff equations of these sorts and i will be reading that to try and find an exact solution to the problem, that would make it i think a lot easier for when i put all the code onto an arduino and have it log data and correct its path, that will happen much later though. Do you guys happen to have any idea on how it is done in real spacecraft, do they use really good approximations or do they try to solve the equations exactly? I ve been looking into what might make in the futurue a good business, small disposable launchers capable of taking multiple cubesats into orbit and for a decent price, also giving flexibility to the companies wishing to use them. Once again,

Thank you very much!

[K^2]

Small rockets are much harder to do efficiently than large ones, unfortunately. That's why you don't see a lot of dedicated cubesat launchers. They are usually taken "along for a ride" with heavier cargo.

In terms of what they really use for real rockets, they are essentially the same equations, but the methods for solving them are more complicated. A general book on diff equations probably won't be much help on numerical methods, but there are books on numerical methods for solving differential equations. There is good info on Wikipedia, too. If you look up Runge-Kutta Methods, you'll find some good info to get you started.

The other big difference is drag model. I'm pretty sure your rocket will stay strictly sub-sonic. Once you get into transonic regions, drag model changes a lot. Working out exactly how your rocket is going to behave as it goes from Mach .9 to Mach 1.1 is going to be very difficult. But that can be a fairly quick section, especially if you throttle up to punch through transonic faster. But then the supersonic flight will have a different drag coefficient, and as you start getting into hypersonic, things start to change again. This is something that a lot of numerical work would get into. The good news is that if your rocket is perfectly cylindrical, especially if your upper stages have no exposed fins, you can build a fairly simple numerical solver that will give you very good results. But you do need to understand a lot of hydrodynamics to write the code for it. There might be some libraries or even complete programs for it out there, though.

And, of course, you'll need more realistic fuel flow solution. In that sample code, I've assumed constant fuel flow for constant thrust and constant I_SP. That's not going to be the case for a real rocket. Be it solid fuel, hybrid, or liquid rocket, you will have variations in fuel flow throughout ascent. So you'll have a separate differential equation for mass of the rocket, which you'll be solving together with equations of motion as a system. It's not that much more difficult, but one is the first order and the other is second order, and it's something you need to understand how to deal with. If you understand completely how Runge-Kutta methods work, though, and how to apply them to higher order equations, you'll be able to figure it out, even if you end up using a more complex method as a basis.