Predicting model rocket height (with and without drag)

[brainpushups]

Hello all,

I recently made some model rockets and thought it would be fun to predict their height before setting them off this weekend. I thought it would also be interesting to compare the difference in the predictions of a model with no drag and a model with quadratic drag. Surprisingly, the model without drag predicted a lower height. I probably made some error somewhere, but I don't see it and I was wondering if anyone could offer some insight. I'll provide all the information and hope there isn't some simple answer that I've overlooked. Thanks in advance for anybody even willing to look at this rather lengthy post.

The simplifying assumptions I've made for both models during the boost phase are:
1) Use the average thrust of the engine
2) Use the average mass of the rocket+engine+fuel

The thrust of the engine is 4.2N
The average mass of the rocket is 0.08kg
The burn time of the engine is 1.1s

No Drag

This is a simple exercise in kinematics. The equation of motion is T-mg=ma where T is the average thrust and m is the average mass. Solving the differential equation predicts a velocity of 46m/s and an altitude of 26m when the engine has used up its fuel.

The equation of motion during the coasting phase is ma=-mg. Solving this for the total height gives 136m.


Quadratic Drag

My drag coefficient (1/2 rho * Cd *Area) is equal to 0.00035 kg/m

The equation of motion is ma=T-mg-cv^2. To simplify the integration I solved for the terminal velocity using this equation which is about 97m/s. The solution to the integral solving for the velocity of the rocket at burnout is
Subscript[v, b]=Subscript[v, ter]((1-E^-k\[Tau])/(1+E^-k\[Tau])) where tau is the burn time of the engine and k is a constant that cleaned up the mess of the variables a bit. The height when the engine burns out was solved similarly and came out as
Subscript[y, b]=
\!\(\*OverscriptBox["m", "_"]\)/(2c) ln (\!
\*SubsuperscriptBox[\(v\), \(ter\), \(2\)]/(\!\(
\*SubsuperscriptBox[\(v\), \(ter\), \(2\)] -
\*SubsuperscriptBox[\(v\), \(b\), \(2\)]\)))
(hope this is readable, I'm cutting and pasting from another document).

Anyway the predictions are close to that of the burnout phase for no air resistance. This model gives 43m/s and 25m for the velocity and height.

For the coast phase the equation of motion is ma=-mg-cv^2. Once again I solved for the terminal velocity to simplify the equations (different than the vter above)

Going through the process to find the height gives the expression
Subscript[m, c]/(2 c) ln(-\!
\*SubsuperscriptBox[\(v\), \(ter\), \(2\)]/(\!\(
\*SubsuperscriptBox[\(v\), \(b\), \(2\)] +
\*SubsuperscriptBox[\(v\), \(ter\), \(2\)]\)))

where m is now the mass of the rocket without fuel, c is my drag coefficient, vter is now equal to about 46m/s and vb is the velocity at burnout. Anyway, the total height predicted is now close to 270m.

 

[brainpushups]

I see that my cut and pasting didn't work.

during engine burn:

vb=vter((1-Exp[-kt])/(1+Exp[-kt]) - velocity

h = (m/(2c)) Ln(vter^2/(vter^2-vb^2))


during coast

h = m/(2c) Ln(-vter^2/(vb^2-vter^2))

 

[hikaru1221]

I see something strange in your calculations of the 2nd case where there is drag. The terminal speed in the 2nd phase is supposed to be zero, since drag and gravity should cease the upward motion. I got 69m for the additional height contributed by the 2nd phase.

 

[brainpushups]

hikaru,

thanks for the reply. Perhaps my labeling of terms was misleading. What I'm calling terminal velocity for the coasting phase comes from solving the equation ma=-mg-cv^2 for when the acceleration is zero so that vter = Sqrt(mg/c). This would be the terminal velocity for an object in free fall with quadratic drag. The only reason for this was to simplify the number of symbols for integration.

 

[hikaru1221]

How do you solve that equation, when the left side is equal to 0, and the right side is always negative? Why is there free fall here, when we are to find the height, which means we only care about the upward motion? Besides, without caring/knowing about the exact equation, in a physical sense, drag and gravity always prevent the upward motion, right?

 

[brainpushups]

I'll show you more explicitly.

During the boost stage the equation of motion is m(dv/dt)=T-mg-cv^2. To solve this differential using separation of variables you would get

(m dv)/(T-mg-cv^2) = dt

This can be cleaned up a bit by plugging in what I've termed the terminal velocity for the boost stage (set the acceleration in the equation of motion equal to zero and solve for v). During the boost stage this is vterm=Sqrt((T-mg)/c). You can see that a substitution can be made to make the equation look like this

(m/c)(dv/(vterm^2-v^2) = dt

Integrating the left hand side from v=0 to v=vboost and the right hand side from 0 to t where t is the burn time of the engine returns.

(m/(2c vterm) Ln ((vterm+vboost)/(vterm-vboost)) = t

Solve this for vboost to get the equation for the final velocity after the engine burns out

vboost = vterm ((1-Exp[-kt])/(1+Exp[-kt])) where k = (2 c vter)/m.


For the coasting phase the equation of motion is different because there is no thrust. This is m(dv/dt) = -mg-cv^2. Once again this can be solved by separation of variables. But first I make a substitution that dv/dt (dh/dh) = dv/dh (dh/dt) = (dv/dh) v to get

(m v dv)/(-mg-cv^2) = dt

But the same simplification can be applied if we set acceleration equal to zero in the equation of motion for the coast phase to see that the terminal velocity is

vtermcoast = Sqrt(mg/c)

Now the separated equation can be written as

(m/c) (v/(-vtermcoast^2-v^2) = dh

So you see that I really haven't changed anything about the equations. Solving for the terminal velocity in the equations of motion just gives me something to call a bunch of constants to simplify the look of the equations. After integration, the equation for the height during the coast phase (by the way, I used the same method to find the height during the boost phase) is

h = m/(2c) Ln(-vter^2/(vb^2-vter^2))

I hope that clarifies my approach.

 

[brainpushups]

hikaru,

I've found my error. I went through the integration of the coasting height without substituting for any new terms. After performing this integration I got the result

coastheight = (-m/2c) Ln[(mg/c)/(mg/c - burnvelocity^2)

And guess what... this model predicts a height of about 69m!

Thanks for your help.

 

[hikaru1221]

I'm not sure if you're on the right track or not. Here is the problem:

But the same simplification can be applied if we set acceleration equal to zero in the equation of motion for the coast phase to see that the terminal velocity is

vtermcoast = Sqrt(mg/c)

We have: mdv/dt = - mg - cv^2.
If dv/dt = 0, then we have: mg = - cv^2 There is a minus sign! In the 2nd phase, the acceleration cannot be equal to zero.

I think this error has nothing to do with substituting the values. Therefore your result, if it is right, may be just a coincidence. Please kindly excuse me if I misunderstood you

 

[lshcue]

Something wrong with your integration

I think your assumption is a creative one and physically meaningful.

During the integration you should consider that the mass is not a constant but a function of time. The mass can be expressed as

m = m0 + rt (r< 0).

where m0 is the initial mass and r is the flow rate of the propellant. This leads to

dm = r dt

Hence the equation "(m/c)(dv/(vterm^2-v^2) = dt" can be rewritten as follows:

(1/c)(dv/(vterm^2-v^2) = dt/m = (1/r)dm/m

If the terms c and r are constants, we can get a solution. ^.-

 

[K^2]

The zero thrust equation can be solved analytically.

\frac{dv}{dt} = -g-kv

Yields. v(t)=-v_ttan(gt/v_t)

v(t) = -v_t tan\left(\frac{g t}{v_t}\right)

Where v_t=sqrt(g/k) is terminal velocity. Notice that this equation gives v(0)=0, so you need to solve for t₀, such that v(t₀=v₀. That gives you t₀=(vt/g)tan^-1(v₀/v_t)

The height can be obtained by simple integration.

h=\frac{v_t^2}{g}ln\left(cos\left(\frac{gt_0}{v_t}\right)\right)Boost stage can be handled similarly. The differential equation above is going to become:

\frac{dv}{dt} = a-kv

Here, a=F/m-g is the acceleration of rocket in vacuum, with F being average thrust and m being average mass. This is also easy to solve.

v(t) = v_t tanh\left(\frac{g t}{v_t}\right)

Notice that it's hyperbolic tangent, and the v_t for this computation is sqrt(a/k). This time, v(0)=0 is your starting condition, so v_f=v(t_f) where t₀ is time of engine operation, will give you v₀ for coasting stage. Distance is again obtained by integration.

h=\frac{v_t^2}{g}ln\left(cosh\left(\frac{gt_f}{v_t}\right)\right)

Again, note the hyperbolic cosine function here.

Try combining these two computations and comparing them to your results. You should be able to get somewhat better estimate. I'm curious to see how close it is to your numbers. (You should have all the coefficients. What I'm calling k is your c/m. So just substitute in, and see what you get.)