Parachute Problem

shock” (j(t_{d}), where j = \frac{da}{dt} is the jerk) is a “heavy but smooth tug”. The mathematical interpretation of this is that the acceleration should be (at least) C^{1} for t > 0. Explicit expressions for the acceleration and jerk can be found directly from the equations of motion:

 > A := subs( k=K[pw_const], rhs(ode(SYS,v)) );
 > J := subs( 'diff(v(t),t)'='a(t)', diff(A,t) );

A := −g− ( k1 + ( k2 − k1 ) Heaviside( t − td ) ) v( t) /m

J := − \frac{( k2 − k1 ) Dirac( t − td ) v( t)}{m}−

( k1 + ( k2 − k1 ) Heaviside( t − td ) ) a( t) /m Thus, the acceleration and jerk do not possess the necessary smoothness (unless k_{1} = k_{2}, i.e., the parachute does not open). Note that the jump in the jerk, i.e.,

[j(t_{d})] := lim_{t→t_{d}^{+}} j(t) − lim_{t→t_{d}^{−}} j(t),

is not simply the coefficient of the Heaviside function since the jump in the acceleration must also be considered. The following procedure provides a simple means to obtain the jump in any expression at a specified time.

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

The two jumps are now found to be

 > Ajump := jump(A, t=td );
 > Jjump := jump(subs(a(t)=A,J), t=td );

Ajump := \frac{v( td ) ( −k2 + k1 )}{m}

Jjump := −( −k2 + k1 )
( k1 v( td ) + m g + k2 v( td ) ) /m^{2}

Using the Nagle/Saff parameter values, the jump in the acceleration is consistent with the graphical solution:

 > subs( v(td)=subs(t=td,V), PARAM[ns], Ajump/g ):
 > 'Ajump' = "*G, `` =evalf(",6)*G;

Ajump = ( 6 − 6 e^{( −12 )} ) G, = 5.99996 G

 > subs( v(td)=subs(t=td,V), PARAM[ns], Jjump/g ):
 > 'Jjump' = " * `G/s`, `` = evalf(",6) * `G/s`;

Jjump = (\frac{48}{5} e^{( −12 )} − \frac{42}{5}) G/s, = − 8.39994 G/s

One way to overcome the deficiencies in the original model is to generalize the form of the coefficient of air resistance. The model proposed in [6] (see also [8]) adds a transition layer to the coefficient of air resistance during which the air resistance will increase from the free-fall value to the final descent value. That is,

k(t) = \{ \begin{matrix}k_{1}, & t < t_{d} \\ k_{d}, & t_{d} ≤ t < t_{d} + τ \\ k_{2}, & t ≥ t_{d} + τ \end{matrix}  (3)

where k_{d} is the coefficient of air resistance during deployment and τ is the length of time required for the chute to be fully deployed. Or, in Maple,

                                                         > Kgen := k1 + (KD-k1)*Heaviside(t-td)
                                                         >  + (k2-KD)*Heaviside(t-td-tau):

Note that choosing k_{d} to be the linear interpolant results in an acceleration that is only continuous. One choice of k_{d} that makes the acceleration continuously differentiable (on t > 0) is the interpolating cubic polynomial whose derivative vanishes at t = t_{d} and t = t_{d}+τ. Maple definitions of these functions are obtained from the following commands:.

                                                         > Kd[linear] := alpha*t+beta:
                                                         > Kd[cubic] := int( alpha*(t-td)*(td+tau-t), t )
                                                         >  + beta:
                                                         > for nam in [ linear, cubic ] do
                                                         > BC := Kd[nam]=kk;
                                                         > solve( { subs(t=td, kk=k1, BC),
                                                         >  subs(t=td+tau, kk=k2, BC) },
                                                         >  { alpha, beta } );
                                                         > DEPLOY := subs( ", Kd[nam] );
                                                         > K[nam] := subs( KD=DEPLOY, Kgen );
                                                         > od:

The form of the coefficient of air resistance is only half of the story. It is also necessary to use realistic values of the parameter. Values consistent with the Air Force Academy Training Guide [10] are

                                                         > PARAM[af] := k1=2/11*75, k2=32/15*75, td=10,
                                                         >  tau=3.2, x0=4000*0.3048, m=75,
                                                         >  g=981/100:

With these constants, the coefficients of air resistance with the linear and cubic interpolants appear as

                                                         > plot(subs( PARAM[af], { K[linear], K[cubic] } ),
                                                         >  t=9..15, k=0..170,
                                                         >  title=`Continuous resistance coefficients`);

Explicit solutions for the linear and cubic coefficients can be found by a three-step process similar to the one

\overline{4}