ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Unit 20 -- Application: Rocket Flight

Outline of Unit 20

Initialization

> restart;

> with( DEtools ):

> with( plots ):

> with( linalg ):

Warning, the name changecoords has been redefined

Warning, the name adjoint has been redefined

Warning, the protected names norm and trace have been redefined and unprotected

>

Problem Statement

A model rocket having initial mass kg is launched vertically from the ground. The rocket expels gas at a constant rate of kg/sec and at a constant velocity m/sec relative to the rocket. Assume the magnitude of the gravitational force is proportional to the mass with proportionality constant g and that the rocket contains kg of gas.

• Formulate a model of the rocket's motion until all of the gas has been expelled.
• Assuming the initial velocity is zero, determine the velocity of the rocket and its height above ground according to this model.
• What model applies after all of the gas has been used?

>

Let denote the height above ground in meters and denote the time since lauch in seconds. Then is the velocity of the rocket. The rocket's mass at time is until ; thereafter, the mass is a constant .

Recall that Newton's second law states that force is equal to the time rate of change of momentum. Thus, for < ,

> mass := m[0] - alpha*t;

> F[gravity] := - mass * g;

> F[thrust] := alpha * beta;

> eomX1 := mass * diff( x(t), t$2 ) - F[thrust] = F[gravity];

> icX1 := x(0)=0, D(x)(0)=0;

>

After all of the gas has been released, the equation of motion simplifies to

> mass2 := eval( mass, t=m[gas]/alpha );

> F[gravity2] := - mass2 * g;

> eomX2 := mass2 * diff( x(t), t$2 ) = F[gravity2];

>

with "initial conditions" at chosen to ensure continuity with the solution at the end of stage 1.

>

Note that the second-order ODE for the height reduces to a first-order linear ODE for the velocity

> eomV1 := eval( eomX1, diff( x(t), t )=v(t) );

> icV1 := v(0)=0;

>

The solution to this IVP is

> solV1a := dsolve( { eomV1, icV1 }, v(t) );

>

which simplifies (somewhat) to

> combine( solV1a, ln, anything, symbolic );

>

Further simplifications are difficult for Maple perform, but hand manipulations quickly lead to

> solV1b := v(t) = -g*t + beta*ln( m[0]/(m[0]-alpha*t) );

>

To check that this function is a solution to the equation of motion for the velocity,

> odetest( solV1b, eomV1 );

>

Integration of the velocity, with , is most easily obtained as the solution to a trivial IVP

> solX1a := dsolve( { eval( solV1b, v(t)=diff(x(t),t) ), x(0)=0 }, x(t) );

>

which simplifies nicely to

> solX1b := collect( solX1a, ln );

>

The only additional modifications that one is likely to want to make to this solution are factoring from the coefficient of the first term and rewriting the argument of the logarithm as .

> solX1c := x(t) = beta*t - 1/2*g*t^2 - beta*(m[0]/alpha-t)*ln( m[0] / (m[0]-alpha*t) );

>

The equations of motion for the rocket after all gas has been expelled are simpler than those for Stage 1.

> eomX2;

> eomV2 := eval( eomX2, diff(x(t),t)=v(t) );

>

The only difficulty to finding the height above ground and velocity of the rocket is finding appropriate initial conditions. The initial conditions for Stage 2 are chosen to ensure continuity with the solutions at the end of Stage 1. That is, at time :

> tstar := m[gas]/alpha;

> icV2 := eval( solV1b, t=tstar );

> icX2 := eval( solX1b, t=tstar );

>

The velocity of the rocket is

> solV2 := collect( dsolve( { eomV2, icV2 }, v(t) ), alpha );

>

while the height above ground is

> solX2 := dsolve( { eval( solV2, v(t)=diff(x(t),t) ), icX2 }, x(t) );

>

The fact that the velocity of the rocket is linear and the height above ground is quadratic is consistent with the observation that the equation of motion dictates that the acceleration is constant. Moreover, note that the velocity is a decreasing function and the height is concave down. The maximum height above ground occurs when where

> Tmax := solve( eval( solV2, v(t)=0 ), t );

>

assuming that > . If this condition is met, then the maximum height of the rocket is

> Xmax := simplify( eval( solX2, t=Tmax ) );

>