Equation of motion with linear acceleration and velocity-squared drag

[darrencubitt]

Homework Statement

This isn't exactly homework, but this seemed like the right place for this question. I'm working on an add-on for Orbiter - the space flight simulator - and would like to be able to determine the equation of motion for an object traveling at an initial velocity (v₀), with a linear acceleration term (a), a velocity-squared drag proportionality constant (k), and elapsed time (t).

Homework Equations

v = v + a * t - k * v^2 * t

The Attempt at a Solution

Being of limited mathematical ability, I scoured the web for some solutions. I found many dealing with velocity-proportional drag, and also some velocity-squared functions which, unfortunately, assume v0=0. Reluctantly, I tried my own hand and I managed to hammer out an equation which can deal with velocity and drag but, alas, no acceleration.

Looking back at how I came to this solution, I am sure that my logic is utterly, utterly wrong. This is almost certainly a case of being right for the wrong reasons. For your entertainment (and also in the hope of getting an answer), my working follows:

I started by writing the basic equation:

1. v_t = v₀ - k * v₀^2 * t

Obviously, since velocity is constantly changing through time, the (k * v₀^2) term will only be valid at t=0, so I divided both sides by v^2 to remove it from this part of the equation:

2. 1/v_t = (1/v₀ - k * t)

I then reasoned that since I am now dealing with the reciprocal of the velocity, I should add the drag term rather than subtract it.

3. 1/v_t = (1/v₀ + k * t)

Finally, I multiplied the reciprocals out:

4. v_t = v₀ / (1 + k * t * v)

And, realising velocity may be negative, used the magnitude of the velocity on the right hand side.

5. v_t = v₀ / (1 + k * t * |v|)

After studying the relationships between the initial and final velocity, I concluded that displacement could be found by:

6. LOG(1 + k * t * |v|) / k (for k != 0)

I guess you could call eqn. 6 "solving by imperical observation".

I am sure that my workings are hilariously wrong but, to my surprise, these equations actually work. I compared the equation to the output of a numerical approximation and it agrees almost 100%, with some slight deviations which I put down to the latter being an approximation.

I tried repeating this trick when including a linear acceleration term (a) but always end up with a nasty (a/v^2) term stuck inside the equation (replacing the k*v^2 problem I was initially trying to solve).

Can anyone provide me with any hints on how to perform this calculation with linear acceleration included? I have seen some solutions bordering on what I am looking for, but either assume v₀=0 or break down in situations with negative acceleration, etc. Is what I am looking for even possible?

Thanks for any help.

 

[Filip Larsen]

In general you will have to solve a differential equations such as

$$m\frac{dv}{dt} = -m g + c v^2$$

which describe a particle of mass m falling vertical down in a fluid with constant gravity and constant drag coefficient. Simple differential equations like this can be solved analytically, but in general you should not expect this to be possible.

For instance, if you were to model the dynamics of an vehicle entering the atmosphere with orbital speed, you would need a more detailed model that included atmospheric density variations and inverse square gravity (as opposed to simple uniform gravity), and you would possible also have to consider vehicle lift too.

Perhaps you could explain a bit more what you want to model?

If you are trying to model atmospheric re-entry and want to understand the models yourself and do not mind a bit of calculus, I can recommend [1]. Simpler models, like the one mentioned above, are typically part of most books on mechanics, like for instance [2], but that also required a bit of calculus.[1] Dynamics of Atmospheric Re-Entry, Regan, AIAA 1993.
[2] Analytical Mechanics, 4th, Fowles, CBS College Publishing, 1983.

 

[darrencubitt]

For instance, if you were to model the dynamics of an vehicle entering the atmosphere with orbital speed, you would need a more detailed model that included atmospheric density variations and inverse square gravity (as opposed to simple uniform gravity), and you would possible also have to consider vehicle lift too.

Perhaps you could explain a bit more what you want to model?

Hi, Filip. Thanks for the reply.

In this particular case, I am attempting to model a crafts horizontal velocity over the surface of a planet (with atmosphere) as it is on a landing approach. I would like to be able to determine the exact (or close enough) time the craft will have to fire its braking rockets to come to a complete stop over the landing site (ignoring the vertical descent portion).

I can use the 'coasting drag' equation (v / (1 + v * t * k)) to predict when the craft will reach the landing site (since it is not accelerating horizontally), and at what velocity, with quite a good degree of precision. Unfortunately, when I try to solve for the required rocket burn-time without taking into account drag (which is all I can do at the moment), it severely over-estimates the time required.

Thanks again.

 

[Filip Larsen]

The equation you gave is the solution to the differential equation $\dot{v} = -kv^2$. By adding a constant brake acceleration of magnitude $a$ this becomes

(1) $$\dot{v} = -kv^2 - a$$

which can be integrated to

(2) $$\tan^{-1}(v_0\sqrt{k/a})-\tan^{-1}(v\sqrt{k/a}) = t \sqrt{ak}$$

The time it takes to brake from $v_0$ to zero speed can then be found as

(3) $$t_b = \frac{1}{\sqrt{ak}} \tan^{-1}(v_0 \sqrt{k/a})$$

To be really useful I guess you would also need to know the distance covered during the brake segment. However, my integration skills are not sufficient to handle equation (2) without a more serious effort (I'm not even sure if it can be integrated). I will give it a try, but let's hope someone else can jump in with an integration of (2) to get a relation between distance and time.

 

[Filip Larsen]

Looks like equation (2) can be integrated to give the relation

(4) $$ac^2 s = \log(\cos(ac(t_b-t))) - \log(\cos(act_b))$$

where $c = \sqrt{k/a}$. Inserting $t = t_b$ gives the breaking distance

(5) $$s_b = -\frac{1}{ac^2}\log(\cos(act_b))$$

If the vehicle is coasting with speed $v_0$ and has distance to target $s_0 = s_c + s_b$ consisting of a coasting distance and a breaking distance, and the breaking starts when $v = v_b$, then the coasting distance is

(6) $$s_c = -\frac{1}{k}\log(v_b/v_0)$$

which added to (5) gives

(7) $$-\frac{1}{ac^2}\log(\cos(\tan^{-1}(v_bc)))-\frac{1}{k}\log(v_b/v_0) = s_0$$

Looks like this one has to be solved numerically for $v_b$. Assuming a solution can be found, this can then be inserted in the original coast speed equation

(8) $$v = v_0 \frac{1}{1+k v_0 t}$$

to give the time to coast as

(9) $$t_c = \frac{v_0 - v_b}{v_b}\frac{1}{kv_0}$$

All of the above completely untested and may easily contain errors on my part, but this is the best I can do this evening. The late hour is also why the explanations are somewhat lacking, but feel free to ask and I may pick it up later tomorrow.

 

[Filip Larsen]

Finally got time to go all the way through. Let me make a summary post including solving the last equation.

We have a vehicle in horizontal flight with the instantaneous speed $v_0$ heading directly towards a landing point at the distance $s_0$. The acceleration of the vehicle is given by the differential equation

(1) $$\dot{v} = -k v^2 - a$$

where $k$ is the total constant drag coefficient and $a$ is the constant deceleration. Note, that $k$ is considered independent of $v$ which means the vehicle does not use increased lift to remain in horizontal flight while slowing down but instead utilize a vertical thruster or similar.

In the first segment of approach to the landing site the vehicle coast while only being decelerated by drag. Setting $a = 0$ and integrating (1) twice gives us the distance covered for a coast segment where the speed slows from $v_1$ to $v_2$ as

(2) $$s_{12} = -\frac{1}{k}log(v_2/v_1)$$

or, equivalently,

(2b) $$v_2 = v_1 \exp(-k s_{12})$$

The relation between speed and time in a coast segment is

(2c) $$v_2 = \frac{v_1}{1+k v_1 t}$$

Note, that since $s_{12} \rightarrow \infty$ for $v_2 \rightarrow 0$ the vehicle will never reach zero speed, but if $v_m$ is the lowest speed we want to consider during coasting, then the maximum range of coasting is

(3) $$s_{max} = -\frac{1}{k}log(v_m/v_0)$$

When the vehicle enters the second segment of its approach to the landing site it adds a breaking deceleration $a > 0$. Integrating (1) once with speed varying from $v_1$ to $v_2$ gives the relation

(4) $$ac t = \arctan(v_1 c) - \arctan(v_2 c)$$

where $c = \sqrt{k / a}$. For instance, if the vehicle breaks from speed $v_1 = v_b$ down to $v_2 = 0$, this will be related to the breaking time $t_b$ by

(5) $$ac t_b = \arctan(v_b c)$$

Integrating (4) to get distance with the speed varying from $v_b$ to zero gives the breaking distance

(6) $$s_b = -\frac{1}{k}\log(\cos(\arctan(v_b c)))$$

If we insert the vehicles current speed $v_0$ instead of $v_b$ in (6) we get the minimum distance the vehicle need to break to a full stop

(7) $$s_{min} = -\frac{1}{k}\log(\cos(\arctan(v_0 c)))$$

(continued in next post)

 

[Filip Larsen]

(continued)

Now, at any given time with the vehicle having speed $v_0$ and distance $s_0$ from landing, with $s_{min} < s_0 < s_{max}$, we would like to know how much time the vehicle should coast before starting breaking at speed $v_b$ such that the vehicle ends with zero speed over the landing site. The distance to the landing site is the sum of the distance to coast and distance to break, that is

(8) $$s_0 = s_c + s_b$$

which by using (2) and (6) and some transformations gives the relation

(9) $$v_b \cos(\arctan(v_b c)) = v_0 \exp(-k s_0) = v_c$$

This relation determines the speed $v_b$ at which the vehicle should start breaking given either the current speed and distance to landing site or, alternatively by comparing with (2b), given the residual speed $v_c$ which is the speed the vehicle would have had if it started with speed $v_0$ and coasted the distance $s_0$. The equation (9) can be solved for $v_b$ to give

(10) $$v_b = \frac{v_c}{\sqrt{1-c^2v_c^2}}$$

which has a real number solution when $c^2v_c^2 < 1$ or, in terms of the acceleration involved, when

(11) $$a > k v_c^2 = k v_0^2 \exp(-2 k s_0)$$

Finally, for a given $v_b$ the time to coast is then given by (2c) as

(12) $$t_c = \frac{1}{k v_0}\frac{v_0 - v_b}{v_b}$$

and the time to break by (5) as

(13) $$t_b = \frac{1}{ac}\arctan(v_b c)$$

 

[darrencubitt]

Holy cow, Filip! Thanks for providing all this information! I can't even begin to explain how much I appreciate this.

The best I had done so far was to find http://www.grc.nasa.gov/WWW/K-12/airplane/flteqs.html" which gives some equations for drag-impeded flight. I ended up using the equation for displacement and solving numerically for acceleration, but I will certainly take a look at your equations in the morning when I am not so sleep deprived!

There are some other things I would like to implement (such as calculating crossrange distance), which I am sure this information will be invaluable in solving. One thing that you have really helped me with is that in all the equations I have found on the web, all of them use the term $$\sqrt[]{a/k}$$. This turns out to be the terminal velocity for a given a and k, but it is interesting to see how this comes about. I feel more comfortable using equations like this when I at least have some basic understanding of why they work, rather than just taking it on faith.

If you are at all interested in seeing what you have helped me out with, I have uploaded a http://www.youtube.com/watch?v=1CyW8aVBwB8" demonstrating a landing on Mars using pure mathematics to work out when to fire the retro rockets. As you can see, once drag is taken into account, the timing is almost perfect. It is hard to see in this video, but the retro burn starts almost a full 2 minutes before reaching the landing site. I was quite amazed that the equations worked out right to within a few meters, given the other variables which haven't been taken into account (such as the drag coefficient changing with velocity - which I know Orbiter does actually calculate).

Thanks again, Filip. I can not even begin to express my gratitude.