>

Differential Equations in the New Millennium: The Parachute Problem

Douglas B. Meade (meade@math.sc.edu)

Department of Mathematics

University of South Carolina

Columbia, SC 29229

Allan A. Struthers (struther@math.mtu.edu)

Department of Mathematical Sciences

Michigan Technical University

Houghton, MI 49931

April 29, 1999

>

This worksheet is a CAS supplement to a paper submitted to a special Mathematics Education of Engineers issue of the International Journal of Engineering Education (co-edited by Prof. Shirley Pomeranz, Univ. of Tulsa). A parallel Mathematica notebook is also available.

>

> restart;

> with( plots ):

>

Auxiliary Procedure Definitions (must be executed -- just press return)

As the system of ODEs and initial conditions are both sets, it will be useful to have a means to extract a specific entry from either object.

> ic := (IVP,y0) -> op( select( has, IVP, y0 ) ):
ode := (IVP,y) -> op( select( has, IVP, diff(y(t),t) ) ):

>

The position (height above ground level), velocity, acceleration, and jerk (rate of change of acceleration) provide a good description of the motion of the skydiver. Given an initial value problem for the rate of change of the velocity and the initial position, the procedure motion uses dsolve,numeric to compute a numerical approximation to the position, velocity, acceleration, and jerk of the object. The result of a call to motion is a list of procedures. plot_motion uses the information computed by motion to create a single plot containing any combination of the position, velocity, acceleration, and jerk. The plot returned by plot_motion is created with odeplot .

Note that it may be necessary to exercise special care near discontinuities of the solution and/or ODE. In order to "see" discontinuities in the solutions, a point plot is recommended; to obtain a point plot, simply add style=POINT as an optional argument to the plot_motion command. For finer resolution, I also add numpoints=200 .

>

> motion := proc( IVP::set(equation), t::name )
local ODEpos, auxIVP, soln;
options `Copyright 1999 by Douglas B. Meade`;
ODEpos:=diff( x(t), t ) = v(t);
auxIVP := { ODEpos };
soln:=dsolve(IVP union auxIVP, {x(t),v(t)},
type=numeric, output=listprocedure ):
RETURN( soln );
end:

> plot_motion := proc( soln::list,
time_int::{range(numeric),list(range(numeric))},
scale::set(equation),
IVP::set(equation) )
local acc, jerk, pos, vel, CURVES, S, TI, TIME_INT, opts;
options `Copyright 1999 by Douglas B. Meade`;
if nargs>4 then opts:=args[5..-1] else opts:=NULL fi;
if map(lhs,scale) minus {x,v,a,j} <> {} then
ERROR("invalid argument; scale factors can be specified only for x, v, a, and j, but received",scale);
fi;
soln(0); # always begin the plot by integrating from the initial point
pos := x=x(t);
vel := v=v(t);
acc := a=solve( ode(IVP,v), diff(v(t),t) );
jerk := j=rhs( diff( acc, t ) );
jerk := subs( 'diff(x(t),t)'=rhs(vel), 'diff(v(t),t)'=rhs(acc), jerk );
CURVES := subs( {pos, vel, acc, jerk}, [seq( [t,lhs(S)/rhs(S)], S=scale)] );
RETURN( plots[display]( [seq( plots[odeplot]( soln, CURVES, TI, opts ),
TI=convert(time_int,list) )] ) );
end:

>

The ``jump'' in a function at a given point is the difference between the left and right limits of the function at that point.

> jump := ( f, pt ) -> limit( f, pt, right ) - limit( f, pt, left ):

>

Equations of Motion for Linear and Quadratic Models

> IC := { x(0)=x0, v(0)=0 };

> EOM[lin] := { diff(x(t),t) = v(t),

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

> EOM[quad] := { diff(x(t),t) = v(t),

> m*diff(v(t),t) = -m*g + k*v(t)^2 };

>

Summary of Solutions and Terminal Velocity for Linear and Quadratic Models

Explicit formulae for the velocity and position can be obtained using dsolve . While the dsolve command could be used, I have been unable to force Maple to present the solution to the system as nicely as I am able to obtain by first solving for the velocity and then the position. The acceleration and jerk can also be expressed directly in terms of the velocity. The theoretical terminal velocity is also easily obtained. Note that the velocity should be negative. It might be possible to use the assume facility to impose conditions that help Maple simplify some of the results, but I prefer to avoid using assume if possible.

> i := 'i':

> for i in [lin,quad] do

>

> if i=lin then printf("\n\nLINEAR MODEL\n============\n")

> elif i=quad then printf("\n\nQUADRATIC MODEL\n===============\n")

> fi;

>

> # S := dsolve( EOM[i] union IC, {x(t),v(t)} );

> # X[i] := collect( allvalues( subs( S, x(t) ) ), {m,g,k} );

> # V[i] := normal( allvalues( subs( S, v(t) ) ) );

> V[i] := rhs( normal( allvalues( dsolve( { ode(EOM[i],v), ic(IC,v(0)) }, v(t) ) ) ) );

> X[i] := rhs( collect( dsolve( { subs( v(t)=V[i], ode(EOM[i],x) ), ic(IC,x(0)) }, x(t) ), {m,g,k} ) );

> A[i] := collect( simplify( subs( {v(t)=V[i], x(t)=X[i]}, solve( ode(EOM[i],v), diff(v(t),t) ) ) ), {m,g} );

> J[i] := subs( solve(ode(EOM[i],v),{diff(v(t),t)}), diff(subs(k=k(t),A[i]),t) );

>

> ACC[i] := rhs( op( solve( subs({diff(v(t),t)=a(t),k=k(t)}, ode(EOM[i],v)), {a(t)} ) ) );

> JERK[i] := collect( subs( diff(v(t),t)=a(t), diff( ACC[i], t ) ), {m,g,k(t),diff(k(t),t)} );

>

> Vterm[i] := rhs( ode(EOM[i],v) ) = 0;

>

> printf("\nTerminal Velocity\n" );

> print( v[T]=solve(Vterm[i],v(t)) );

>

> printf("\nSolutions: expressed in terms of position and velocity" );

> print( a=ACC[i] );

> print( j=JERK[i] );

>

> printf("\nSolutions: expressed explicitly in terms of time" );

> print( x=X[i] );

> print( v=V[i] );

> print( a=A[i] );

> print( j=J[i] );

> od:

LINEAR MODEL

============

Terminal Velocity

Solutions: expressed in terms of position and velocity

Solutions: expressed explicitly in terms of time

QUADRATIC MODEL

===============

Terminal Velocity

Solutions: expressed in terms of position and velocity

Solutions: expressed explicitly in terms of time

>

Parameters for Skydiver, T-10 Canopy, Physics, ...

This section contains values for all parameters in the model.

> param0 :=

> g=9.81, # gravitational constant (kg m/s^2)

> rho = 1, # density of air (kg/m^3)

> mu=1.5*10^(-5), # viscosity of air (kg/m/s

> h=70 * 0.0254, # height of skydiver (m)

> mb=190/2.2, # mass of skydiver (kg)

> b0=0.5, # spread-eagle cross-sectional area (m^2)

> b1=0.1, # feet-first cross-sectional area (m^2)

> d0=35 * 0.3048, # nominal diameter of canopy (m)

> mh=10/2.2, # mass of harness (m)

> mc=13.85/2.2, # mass of canopy and suspension lines (m)

> beta1=0.15, # amount of overinflation (%)

> p=1, # exponent of line extension interpolant

> t0=10, # time ripcord is pulled (s)

> Dt1=0.5, # time for lines to extend (s)

> Dt2=1.0, # time for initial fill of canopy (s)

> Dt3=1.7, # time for remainder of deployment (s)

> x0=4000 * 0.3048, # initial height (m)

> NULL ;

>

> param1 := op(subs( param0, [

> l = 0.84*d0, # length of suspension lines

> dp = 0.7*d0, # projected nominal diameter

> m = mb+mh+mc, # total mass of person and equipment

> t1 = t0+Dt1, # time when snatch force felt

> t2 = t0+Dt1+Dt2, # time when opening force felt

> t3 = t0+Dt1+Dt2+Dt3, # time when deployment is complete

> NULL ] ));

>

> param2 := op(subs( param0, param1, [

> a1 = Pi*(dp/2)^2, # projected nominal area

> alpha0 = (1.95*b0+0.35*b1*(l-h))/1.33, # initial value for Ae12

> NULL ] ));

>

> param3 := op(subs( param0, param1, param2, [

> beta0 = log(a1/alpha0), # `growth rate' in Ae12

> NULL ] ));

> param4 := op(subs( param0, param1, param2, param3, [

> Ae12 = alpha0*exp(beta0*(t-t1)/(t2-t1)), # canopy area for [t1,t2]

> Ae23 = a1 * (1+beta1*sin(Pi*(t-t2)/(t3-t2))), # canopy area for [t2,t3]

> NULL ] ));

>

> param := { seq( param.i, i=0..4 ) };

>

> paramFF := select( has, param, {m,g} ) union {

> VtermFF=-100 * 5280*1282.54/100/60/60, # terminal free-fall velocity: 100 m/hr

> VtermD =-16 * 12*2.54/100, # rate of descent for T-10: 16 ft/sec

> x0=1.52, # landing impact comparable to jump from 5 ft (1.52m) wall

> NULL };

>

Computation of Drag Coefficients

Assuming the free-fall terminal velocity is known, the free-fall drag coefficient is found by setting the RHS of the EOM to zero and solving for . The rate of descent for the T-10 parachute is given as 16 ft/sec [ebk, p. 86]. This suffices to compute the drag coefficient when the parachute is fully deployed. A second method for finding the drag coefficient when the parachute is fully deployed is to use the characterization of the impact force as being equivalent to jumping from a five-foot high wall. To find this impact velocity it is necessary to solve each model with an initial height of 5 ft (1.52m), find the time when the height is zero, then determine the velocity at that time.

It would not be unreasonable to make these computations into a function of the values in param2 .

>

> i := 'i':

> k := 'k':

> for i in [lin,quad] do

>

> if i=lin then printf("\n\nLINEAR MODEL\n============\n")

> elif i=quad then printf("\n\nQUADRATIC MODEL\n===============\n")

> fi;

>

> PARAM[FF,i] := solve(subs( v(t)=VtermFF, paramFF, Vterm[i] ), {k});

> PARAM[D ,i] := solve(subs( v(t)=VtermD , paramFF, Vterm[i] ), {k});

>

> tt := fsolve( subs( paramFF, PARAM[FF,i], X[i] )=0, t, 0..2 ); # time of impact for jump from 5-foot wall

> vv := evalf( subs( t=%, paramFF, PARAM[FF,i], V[i] ) ); # corresponding impact velocity

> PARAM[V,i,5] := { t=tt, v=vv };

> subs( v(t)=vv, paramFF, Vterm[i] );

> PARAM[D,i,5] := solve(%, {k});

>

> printf( "\nDrag Coefficients\n" );

> print( k[FF]=subs( PARAM[FF,i], k ) );

> print( k[D]=subs( PARAM[D,i], k ) );

> print( k[D,5]=subs( PARAM[D,i,5], k) );

> printf( "\nImpact Velocity (from 5 ft wall)\n" );

> print( v[impact,5]=subs( PARAM[V,i,5], v) );

> od:

LINEAR MODEL

============

Drag Coefficients

Impact Velocity (from 5 ft wall)

QUADRATIC MODEL

===============

Drag Coefficients

Impact Velocity (from 5 ft wall)

>

Analysis of the Model

Drag Coefficient

Specification of the quadratic model is complete once the drag coefficient is implemented. From this information it is possible the smoothness of the solution can be determined.

>

> kk := 1/2 * rho * piecewise( t<=t0, 1.95*b0,

> t<=t1, 1.95*b0 + 0.35*b1*l*((t-t0)/(t1-t0))^p,

> t<=t2, 0.35*b1*h + 1.33*Ae12,

> t<=t3, 0.35*b1*h + 1.33*Ae23,

> t> t3, 0.35*b1*h + 1.33*a1 );

>

With the parameter values suggested in the text, the drag coefficient becomes completely well-defined

> K := simplify( subs( param, kk ) );

> dK := diff( K, t );

>

Notice that the drag coefficient during final descent is approximately 29. This value is a little lower than eiter prediction based on the terminal velocity. This suggests that the landing velocity is likely to be a little higher than that for a jump from a 5-foot high wall and probably noticeably higher than the rate of descent published in [ebk, p. 86]. The precise amount of this over-estimate remains to be seen. The drag coefficient during free-fall is much closer to the value obtained from the literature.

>

The graphs of and its derivative provide a first test of the smoothness of this function. However, be very careful before using a picture as a proof of continuity or smoothness.

> plotK := plot( K, t=0..30, discont=true, axes=FRAMED, title="k(t) vs. t ( as given in (3.5) )" ):

> display( plotK );

> plotdK := plot( dK, t=0..30, discont=true, axes=FRAMED, title="k'(t) vs. t" ):

> display( plotdK );

>

A more reliable method of determining continuity of a function at a point is to compare its right- and left-hand limits at the point, i.e., the jump of the function at the point. The jump command defined at the top of this worksheet is designed for exactly this purpose. (This might be a nice place to use a Maple spreadsheet in Release 5.1.)

> #for T in eval([t0,t1,t2,t3],param) do

> # sprintf( "[k(%5.2f)] = %14.7f, [k'(%5.2f)] = %14.7f", T, eval( jump(K,t=T) ), T, eval( jump(dK,t=T) ) );

> #od;

>

Thus, the drag coefficient is continuous but is not differentiable. This means the acceleration will be continuous, but not the jerk.

>

Terminal Velocity

We are most interested in the free-fall and landing terminal velocities, but the terminal velocity is easily found for any time during the jump.

>

> V[term] := solve( Vterm[quad], v(t) )[2];

>

> Vt := simplify( subs( k=kk, param, V[term] ) );

The terminal velocity during free-fall and final descent are consistent with our earlier observations. The free-fall terminal velocity is extremely close to 100 mph. The final descent terminal velocity does not match reality nearlyy as well; as expected, the landing speed is higher than the speeds found in the literature.

>

Observe that the continuity of the drag coefficient says that the terminal velocity must be continuous for all time. This is not immediately apparent from the explicit definition of the terminal velocity or its graph.

> plot( Vt, t=0..20 );

>

Graphical Verification of Smoothness and Deployment

>

> IVP := subs( k=kk, param, EOM[quad] union IC );

>

> SOLN := motion( IVP, t );

>

> scale[X] := { x=1 }:

> scale[XV] := { x=1000, v=10 }:

> scale[VAJ] := { v=10, a=981/100, j=981/100 }:

> scale[XVA] := { x=1000, v=10, a=981/100 }:

> scale[VA] := { v=10, a=981/100 }:

> scale[all] := { x=1000, v=10, a=981/100, j=981/100 }:

> OPTS := style=POINT, symbol=CIRCLE, numpoints=200:

>

> interface( plotdevice=default );

> PLOT1 := display(

> {

> plot_motion( SOLN, [0..30], scale[VA], IVP, numpoints=200, xtickmarks=3 ),

> textplot( [ [ 14, -2.0, `v/10` ],

> [ 14, 2.5, `a/g` ] ] )

> },

> title="Initial Motion for T-10" ):

> display( PLOT1 );

>

> PLOT2 := display(

> {

> plot_motion( SOLN, [9.7..13.5], scale[VAJ], IVP, OPTS, xtickmarks=3 ),

> textplot( [ [ 10.6, -5.5, `v/10` ],

> [ 11.6, 5.0, `a/g` ],

> [ 10.4, 6.2, `j/g` ] ] )

> },

> title="T-10 Canopy Deployment" ):

> display( PLOT2 );

>

> PLOT3 := display(

> {

> plot_motion( SOLN, [10.25..11.75], scale[VAJ], IVP, OPTS, xtickmarks=3 ),

> textplot( [ [ 10.7, -5.9, `v/10` ],

> [ 10.8, 3.8, `a/g` ],

> [ 10.6, 7.7, `j/g` ] ] )

> },

> title="Snatch and Opening Forces" ):

> display( PLOT3 );

>

Estimation of Landing Time

> PLOT4 := display(

> {

> plot_motion( SOLN, [0..180], scale[XV], IVP, xtickmarks=3 ),

> textplot( [ [ 40, 1.1, `x/1000` ],

> [ 40, -0.9, `v/10` ] ] )

> },

> title="Position and Velocity: The First 3 Minutes" ):

> display( PLOT4 );

>

> PLOT5 := display(

> {

> plot_motion( SOLN, [150..170], scale[X], IVP, xtickmarks=3 ),

> textplot( [ [ 152, 60, `x` ] ] )

> },

> title="Estimation of Landing Time" ):

> display( PLOT5 );

>

From the last plot, it appears that landing occurs near t=162 seconds.

> SOLN( 0 ): # reset numerical integrator to initial time

> SOLN( 162 );

>

Since, by this time, the motion is essentially linear, an improved estimate of the landing time can be obtained by linear interpolation:

> t_new := 162 - subs( SOLN(162), x(t)(162)/v(t)(162) );

>

As a final cross-check:

> SOLN( 161.764 );

OK, so it's not exactly zero but this height is only 2 cm, less than one inch!

>

>

Create PostScript Files for Plots

Note: It is necessary to manually modify the Maple-generated PS files. In particular, the translation needs to be changed for each plot (fig1.ps: 126 324 -> 226 324; fig2-5.ps: 216 324 -> 316 224). Also, the bounding box in fig1.ps needs to be changed from 323 119 468 486 to 323 319 468 486.

>

> interface( plotdevice=cps, plotoptions="portrait,noborder,width=5in,height=2in,shrinkby=0" );

> interface( plotoutput="/home/meade/Parachute2/fig1.ps" );

> display( PLOT1, title="" );

> interface( plotdevice=default );

>

> interface( plotdevice=cps, plotoptions="portrait,noborder,width=2.5in,height=2in" );

> interface( plotoutput="/home/meade/Parachute2/fig2.ps" );

> display( PLOT2, title="" );

> interface( plotdevice=default );

>

> interface( plotdevice=cps, plotoptions="portrait,noborder,width=2.5in,height=2in" );

> interface( plotoutput="/home/meade/Parachute2/fig3.ps" );

> display( PLOT3, title="" );

> interface( plotdevice=default );

>

> interface( plotdevice=cps, plotoptions="portrait,noborder,width=2.5in,height=2in" );

> interface( plotoutput="/home/meade/Parachute2/fig4.ps" );

> display( PLOT4, title="" );

> interface( plotdevice=default );

>

> interface( plotdevice=cps, plotoptions="portrait,noborder,width=2.5in,height=2in" );

> interface( plotoutput="/home/meade/Parachute2/fig5.ps" );

> display( PLOT5, title="" );

> interface( plotdevice=default );

>

> interface( plotdevice=cps, plotoptions="portrait,noborder,width=2.5in,height=2in" );

> interface( plotoutput="/home/meade/Parachute2/fig6.ps" );

> display( plotK, title="" );

> interface( plotdevice=default );

>

> interface( plotdevice=cps, plotoptions="portrait,noborder,width=2.5in,height=2in" );

> interface( plotoutput="/home/meade/Parachute2/fig7.ps" );

> display( plotdK, title="" );

> interface( plotdevice=default );

>

>