ORDINARY DIFFERENTIAL EQUATIONS POWERTOOL

Unit 22 -- Application: Parachute Problem

Outline of Unit 22

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 Description

Training jumps for the Parachute Team at the United States Air Force Academy begin 4000 ft above ground level (AGL). The free-fall portion of the jump lasts about 10 seconds; free-falls longer than 13 seconds are grounds for removal from the team. The terminal velocity for free-fall is 120 mph. The landing velocity should be no worse than a free-fall from a 5' wall.

• How long after leaving the plane does the skydiver return to solid ground?

Under the assumption that the force of air resistance is proportional to velocity, the second-order equation of motion for this situation is

> eom2 := m*diff( x(t), t$2 ) = -m*g - k*diff( x(t), t );

>

with initial conditions

> ic2 := { x(0)=x0, D(x)(0)=0 };

>

The equivalent first-order system of ODEs and initial conditions is

> odeX := diff( x(t), t ) = v(t);

> odeV := m*diff( v(t), t ) = -m*g - k*v(t);

> eomXV := { odeX, odeV }:

> icX := x(0)=x0;

> icV := v(0)=0;

> icXV := { icX, icV }:

>

The terminal velocity, , is the equilibrium solution for the velocity ODE. That it, must satisfy

> Vterm := eval( odeV, v(t)=v[T] );

>

The coefficient of air resistance is piecewise-defined with general form

> kk := piecewise( t<10, kFF, kL );

>

The values for the drag coefficients during free-fall, , and landing, , with units kg/s, can be determined from the information given in the problem.

> param0 := g = 9.81,

> m = 190 / 2.2;

>

When the terminal free-fall velocity of 120 mph is converted to units of meters per second

> vFF := -120 * 5280*12*2.54/100/60/60;

>

the coefficient of air resistance during free-fall is found to be

> kFF := solve(eval( Vterm, [param0,v[T]=vFF] ), k );

>

The value of is based on the landing velocity for a jump from a 5' wall.

> x5 := 5 * 12*2.54/100;

>

The position and velocity for a free-fall jump from an initial height of m (= 5') are

> sol5 := eval( dsolve( eomXV union icXV, {x(t),v(t)} ),

> {param0, x0=x5, k=kFF} );

>

The landing time, in seconds, for this short jump is

> t5 := fsolve( eval( x(t)=0, sol5 ), t=0..1 );

>

The corresponding landing velocity

> v5 := eval( rhs( op(select(has,sol5,v(t))) ), t=t5 );

>

is the final ingredient needed to determine the coefficient of air resistance during the time the parachute is deployed

> kL := solve(eval( Vterm, [param0,v[T]=v5] ), k );

>

The complete model is

> sys := eval( eomXV, [param0, k=kk] );

> ic := eval( icXV, x0 = 4000 * 12*2.54/100 );

>

Prior to obtaining an explicit solution to this IVP, note that the velocity satisfies a linear first-order ODE. The graph of the direction field for the velocity equation and the solution curve with provides a first look at the solution to this model.

> DEplot( select(has,sys,diff(v(t),t)), {v(t)}, t=0..20,

> [[0,0]], v=-100..0, stepsize=0.01 );

>

Observe how the velocity approaches the terminal free-fall velocity until the parachute is deployed at . Thereafter, the velocity quickly approaches the new (and slower) terminal velocity.

>

The explicit solution to this model is

> sol := dsolve( sys union ic, { x(t), v(t) } );

>

This solution is quite complicated even with floating-point coefficients :

> evalf( sol );

>

A slightly simpler representation of the solution can be obtained by solving the system twice: first with and the original initial conditions and then with using the position and height from the free-fall solution at the time the parachute is deployed as initial conditions. That is,

> sys1 := eval( eomXV, [param0, k=kFF] );

>

The explicit solution during free-fall is

> sol1 := dsolve( sys1 union ic, { x(t), v(t) } );

>

The equations of motion for the skydiver after the parachute is deployed are

> sys2 := eval( eomXV, [param0, k=kL] );

>

The skydiver's height and velocity at the time the parachute is deployed are

> ic2 := eval( {x(10)=x(10.),v(10)=v(10.)}, eval(sol1,t=10.) );

>

The explicit solution after the parachute is deployed is

> sol2 := dsolve( sys2 union ic2, { x(t), v(t) } );

>

Assembling the solution for both stages of the jump into a single expression leads to

> solX := x(t) = piecewise( t<10, eval( x(t), sol1 ), eval( x(t), sol2 ) );

> solV := v(t) = piecewise( t<10, eval( v(t), sol1 ), eval( v(t), sol2 ) );

>

As expected, this solution is in a much simpler form than the original solution.

>

Note

The explicit solution can be obtained in the same manner without specifying values for the parameters in the problem. This purely symbolic solution is quite complicated, but has a number of uses, including a sensitivity analysis of the solution with respect to each of the parameters.

>

A numeric solution is much easier to obtain and can be used to answer many questions about the motion. To confirm that Maple is computing a reasonable solution, recreate a plot of the velocity function over the first 20 seconds of the jump. (Since it is so easy to do, include a scaled version of the height.)

> solN := dsolve( sys union ic, { x(t), v(t) }, type=numeric ):

> odeplot( solN, [ [t,x(t)/10], [t,v(t)] ], 0..20,

> numpoints=500, legend=["x/10","v"] );

>

Graphical Solution

To see the motion for the entire jump, extend the time interval to 200 seconds.

> odeplot( solN, [ [t,x(t)/10], [t,v(t)] ], 0..200,

> numpoints=500, legend=["x/10","v"] );

>

This suggests that the skydiver lands in just over 3 minutes (180 seconds)

> solN(180);

>

A quick search leads to a better estimate of the impact time

> solN(181.76);

>

Note that the impact velocity is essentially the same as the velocity for a free-fall from a 5' wall

> v5;

>

Analytic Solution

To determine the time of impact from the explicit solution, one might try

> q1 := solve( rhs(solX)=0, t );

> evalf( q1 );

>

but this answer is not physically realistic. Further thought leads to the realization that this "solution" is obtained from the free-fall phase of the solution but the time of impact must occur during the post-deployment phase. To convey this additional knowlege to Maple, use

> Timpact := solve( eval( x(t), sol2 )=0, t );

> evalf( Timpact );

>

and note the agreement with the time of impact obtained from the graphical and numerical solutions.

>

This version of the parachute problem is the classical model found in many ODE textbooks. It's major flaw is that the acceleration is discontinuous at the time the parachute is deployed.

The acceleration can be obtained by differentiation of the velocity

> solA := eval( a(t) = diff( v(t), t ), diff( solV, t ) );

>

The acceleration appears to have a jump discontinuity at the instant the ripcord is pulled

> t_jump := op( discont( rhs(solA), t ) );

>

The size of the jump is

> jump := evalf( limit( rhs(solA), t=10, right )

> - limit( rhs(solA), t=10, left ) );

> `` = eval( jump/g, [param0] ) * g;

>