home
PDF letter size
PDF legal size

Dynamics cheat sheet

September 4, 2016 compiled on — Sunday September 04, 2016 at 10:53 PM

1  Vibration
 1.1  Modal analysis for two degrees of freedom system
  1.1.1  Step one. Finding the equations of motion in normal coordinates space
  1.1.2  Step 2. Solving the eigenvalue problem, finding the natural frequencies
  1.1.3  Step 3. Finding the non-normalized eigenvectors
  1.1.4  Step 4. Mass normalization of the shape vectors (or the eigenvectors)
  1.1.5  Step 5, obtain the modal transformation matrix Φ
  1.1.6  Step 6. Applying modal transformation to decouple the original equations of motion
  1.1.7  Step 7. Converting modal solution to normal coordinates solution
  1.1.8  Numerical solution using modal analysis
 1.2  Fourier series representation of a periodic function
 1.3  Generating Transfer functions for different vibration systems
  1.3.1  Force transmissibility
  1.3.2  vibration isolation
  1.3.3  Accelerometer
  1.3.4  Seismometer
  1.3.5  Summary of vibration transfer functions
 1.4  Solution of Vibration equation of motion for different loading
  1.4.1  common definitions
  1.4.2  Harmonic loading mu ′′ + cu ′ + ku = F sin ϖt
  1.4.3  constant loading mu ′′ + cu′ + ku = F
  1.4.4  No loading (free vibration) mu ′′ + cu′ + ku = 0
  1.4.5  impulse F0δ(t)  loading
  1.4.6  Tree view look at the different cases
  1.4.7  Cycles for the peak to decay by half its original value
  1.4.8  references
2  Dynamics equations, kinematics, velocity and acceleration diagrams
 2.1  Derivation of rotation formula
 2.2  Miscellaneous hints
 2.3  Formulas
 2.4  Velocity and acceleration diagrams
  2.4.1  Spring pendulum
  2.4.2  pendulum with blob moving in slot
  2.4.3  spring pendulum with block moving in slot
  2.4.4  double pendulum
 2.5  Velocity and acceleration of rigid body 2D
 2.6  Velocity and acceleration of rigid body 3D
  2.6.1  Using Vehicle dynamics notations
  2.6.2  3D Not Using vehicle dynamics notations
  2.6.3  Acceleration terms due to rotation and acceleration
 2.7  Wheel spinning precession
 2.8  References
 2.9  Misc. items
3  Astrodynamics
 3.1  Ellipse main parameters
 3.2  Table of common equations
 3.3  Flight path angle for ellipse γ
 3.4  Parabolic trajectory
 3.5  Hyperbolic trajectory
 3.6  showing that energy is constant
 3.7  Earth satellite Transfer orbits
  3.7.1  Bi-Elliptic transfer orbit
  3.7.2  semi-tangential elliptical transfer
 3.8  Rocket engines, Hohmann transfer, plane change at equator
 3.9  Spherical coordinates
 3.10  interplanetary transfer orbits
  3.10.1  rendezvous orbits
  3.10.2  Semi-tangential transfers, elliptical, parabolic and hyperbolic
  3.10.3  Lagrange points
  3.10.4  Orbit changing by low contiuous thrust
  3.10.5  References
4  Dynamic of flights
 4.1  Wing geometry
 4.2  Summary of main equations
  4.2.1  definitions
 4.3  images and plots collected
 4.4  Some strange shaped airplanes
 4.5  links
 4.6  references

1  Vibration

1.1  Modal analysis for two degrees of freedom system

Detailed steps to perform modal analysis are given below for a standard undamped two degrees of freedom system. The main advantage of solving a multidegree system using modal analysis is that it decouples the equations of motion (assuming they are coupled) making solving them much simpler.

In addition it shows the fundamental shapes that the system can vibrate in, which gives more insight into the system. Starting with standard 2 degrees of freedom system

pict

In the above the generalized coordinates are x1   and x2   . Hence the system requires two equations of motion (EOM's).

1.1.1  Step one. Finding the equations of motion in normal coordinates space

The two EOM's are found using any method such as Newton's method or Lagrangian method. Using Newton's method, free body diagram is made of each mass and then F  = ma  is written for each mass resulting in the equations of motion. In the following it is assumed that both masses are moving in the positive direction and that x2   is larger than x1   when these equations of equilibrium are written

pict

Hence, from the above the equations of motion are

pict

or

pict

In Matrix form

⌊        ⌋ (   )   ⌊             ⌋(   )    (     )
  m1   0   { x′′}     k1 + k2 − k2 { x1}    { f1(t)}
⌈        ⌉    1  + ⌈             ⌉       =
   0  m2   ( x′′2)      − k2    k2  ( x2)    ( f2(t))

The above two EOM are coupled in stiffness, but not mass coupled. Using short notations, the above is written as

[M ]{x′′}+  [K ]{x} = {f}

Modal analysis now starts with the goal to decouple the EOM and obtain the fundamental shape functions that the system can vibrate in. To make these derivations more general, the mass matrix and the stiffness matrix are written in general notations as follows

⌊         ⌋ (   )   ⌊        ⌋ (   )    (     )
 m11   m12  { x′′1}     k11 k12  { x1}    {f1(t)}
⌈         ⌉       + ⌈        ⌉       =
 m21   m22  ( x′′2)     k21 k22  ( x2)    (f2(t))

The mass matrix [M ]  and the stiffness matrix [K ]  must always come out to be symmetric. If they are not symmetric, then a mistake was made in obtaining them. As a general rule, the mass matrix [M ]  is PSD (positive definite matrix) and the [K ]  matrix is positive semi-definite matrix. The reason the [M ]  is PSD is that xT[M ]{x} represents the kinetic energy of the system, which is typically positive and not zero. But reading some other references 1 it is possible that [M ]  can be positive semi-definite. It depends on the application being modeled.

1.1.2  Step 2. Solving the eigenvalue problem, finding the natural frequencies

The first step in modal analysis is to solve the eigenvalue problem det([K ]− ω2 [M ]) = 0  in order to determine the natural frequencies of the system. This equations leads to a polynomial in ω  and the roots of this polynomial are the natural frequencies of the system. Since there are two degrees of freedom, there will be two natural frequencies ω1,ω2   for the system.

pict

The above is a polynomial in ω4   . Let ω2 = λ  it becomes

 2
λ (m11m22  − m12m21 )+ λ (− k11m22 + k12m21 + k21m12 − k22m11) + k11k22 − k12k21 = 0

This quadratic polynomial in λ  which is now solved using the quadratic formula. Then the positive square root of each λ  root to obtain ω
  1   and ω
 2   which are the roots of the original eigenvalue problem. Assuming from now that these roots are ω1   and ω2   the next step is to obtain the non-normalized shape vectors φ1,φ2   also called the eigenvectors associated with ω1   and ω2

1.1.3  Step 3. Finding the non-normalized eigenvectors

For each natural frequency ω1   and ω2   the corresponding shape function is found by solving the following two sets of equations for the vectors φ1,φ2

⌊        ⌋     ⌊          ⌋(    )    (  )
 k11  k12     2  m11  m12  { φ11}    { 0}
⌈        ⌉−  ω1⌈          ⌉(    )  = (  )
 k21  k22        m21  m22    φ21       0

and

⌊        ⌋     ⌊          ⌋(    )    (  )
                           {    }    {  }
⌈k11  k12⌉−  ω2⌈m11   m12 ⌉  φ12   =   0
 k21  k22     2  m21  m22  ( φ22)    ( 0)

For ω1   , let φ11 = 1  and solve for

pict

Which gives one equation now to solve for φ21   (the first row equation is only used)

(      2    )      (      2    )
 k11 − ω1m11 + φ21  k12 − ω1m12 =  0

Hence

      − (k11 − ω2 m11)
φ21 = -(-------21--)--
        k12 − ω 1m12

Therefore the first shape vector is

              (             )
     (    )   |             |
     { φ11}   {       1     }
φ1 = (    ) = |(  −((k11−ω21m11))|)
       φ21        k12− ω21m12

Similarly the second shape function is obtained. For ω2   , let φ12 = 1  and solve for

pict

Which gives one equation now to solve for φ22   (the first row equation is only used)

(      2    )      (      2    )
 k11 − ω2m11 + φ22  k12 − ω2m12 =  0

Hence

      − (k11 − ω2 m11)
φ22 = -(-------22--)--
        k12 − ω 2m12

Therefore the second shape vector is

              (             )
     (    )   |             |
     { φ12}   {       1     }
φ2 = (    ) = |(  −((k11−ω22m11))|)
       φ22        k12− ω22m12

Now that the two non-normalized shape vectors are found, the next step is to perform mass normalization

1.1.4  Step 4. Mass normalization of the shape vectors (or the eigenvectors)

Let

μ1 = φT [M ]φ
      1      1

This results in a scalar value μ
 1   , which is later used to normalize φ
 1   . Similarly

μ2 = φT2 [M ]φ2

For example, to find μ1

pict

Similarly, μ2   is found

pict

Now that μ1,μ2   are obtained, the mass normalized shape vectors are found. They are called Φ1, Φ2

             (   )
             {φ11}
                      (     )
      φ      (φ21)    { √φ11μ-}
Φ1 = √--1-=  -√-----=      1
       μ1       μ1    ( √φ21μ1)

Similarly

             (   )
             {φ12}
                      (     )
      φ      (φ22)    { √φ12-}
Φ2 = √--2-=  -√-----=     μ2
       μ2       μ2    ( √φ22μ2)

1.1.5  Step 5, obtain the modal transformation matrix Φ

The modal transformation matrix is the 2× 2  matrix made of of Φ  ,Φ
  1   2   in each of its columns

pict

Now the [Φ]  is found, the transformation from the normal coordinates {x} to modal coordinates, which is called {η} is found

pict

The transformation from modal coordinates back to normal coordinates is

pict

However,   − 1     T
[Φ ]  = [Φ]  [M ]  therefore

pict

The next step is to apply this transformation to the original equations of motion in order to decouple them

1.1.6  Step 6. Applying modal transformation to decouple the original equations of motion

The EOM in normal coordinates is

⌊         ⌋ (   )   ⌊        ⌋ (   )   (      )
 m11  m12   { x′′1}     k11  k12  { x1}   { f1 (t)}
⌈         ⌉ (  ′′) + ⌈        ⌉ (   ) = (      )
 m21  m22     x2      k21  k22    x2      f2 (t)

Applying the above modal transformation {x} = [Φ ]{η} on the above results in

⌊         ⌋    (   )   ⌊        ⌋    (   )   (      )
 m11   m12     { η′′1}     k11 k12     { η1}   { f1(t)}
⌈         ⌉ [Φ ](   ) + ⌈        ⌉ [Φ ](   ) = (      )
 m21   m22       η′′2      k21 k22       η2      f2(t)

pre-multiplying by [Φ]T  results in

     ⌊         ⌋    (   )       ⌊        ⌋    (   )        (     )
      m11  m12      { η′′1}         k11 k12     { η1}        { f1(t)}
[Φ]T ⌈         ⌉ [Φ ]      + [Φ]T⌈        ⌉ [Φ]      = [Φ ]T
      m21  m22      ( η′′2)         k21 k22     ( η2)        ( f2(t))

The result of     ⌊          ⌋
[Φ]T⌈m11   m12 ⌉ [Φ ]
      m    m
        21    22  will always be ⌊    ⌋
⌈1  0⌉
 0  1 . This is because mass normalized shape vectors are used. If the shape functions were not mass normalized, then the diagonal values will not be 1  as shown.

The result of     ⌊        ⌋
      k11  k12
[Φ]T⌈        ⌉ [Φ]
      k21  k22  will be ⌊       ⌋
 ω21  0
⌈       ⌉
  0   ω22 .

Let the result of      (     )
   T {f1 (t)}
[Φ]  (     )
      f2 (t) be (      )
{ f˜1(t)}
(      ) ,
  f˜2(t)  Therefore, in modal coordinates the original EOM becomes

⌊    ⌋ (   )   ⌊       ⌋ (   )   (      )
 1  0  { η′′}     ω2  0   { η1}   { f˜1(t)}
⌈    ⌉    1  + ⌈  1    ⌉       =
 0  1  ( η′′2)     0   ω22  ( η2)   ( f˜2(t))

The EOM are now decouples and each can be solved as follows

pict

To solve these EOM's, the initial conditions in normal coordinates must be transformed to modal coordinates using the above transformation rules

pict

Or in full form

(      )   ⌊          ⌋T ⌊         ⌋ (      )
{ η1(0)}     φ√11- √φ12-    m11   m12  { x1(0)}
         = ⌈  μ1    μ2⌉  ⌈         ⌉
( η2(0))     φ√2μ11- √φ22μ2    m21   m22  ( x2(0))

and

(      )                             (      )
{  ′   }   ⌊ φ√11- √φ12-⌋T ⌊         ⌋ {  ′   }
  η1(0)  = ⌈  μ1    μ2⌉  ⌈m11   m12⌉   x1(0)
( η′(0))     φ√21- √φ22-    m     m    ( x′(0))
   2          μ1    μ2      21   22     2

Each of these EOM are solved using any of the standard methods. This will result is solutions η1(t)  and η2 (t)

1.1.7  Step 7. Converting modal solution to normal coordinates solution

The solutions found above are in modal coordinates η (t),η (t)
 1     2  . The solution needed is x (t),x (t)
 1     2  . Therefore, the transformation {x} = [Φ ]{η} is now applied to convert the solution to normal coordinates

pict

Hence

       -φ11        φ12-
x1(t) = √ μ1-η1(t) + √ μ2η2(t)

and

        φ          φ
x2(t) = √-21-η1(t) + √-22η2(t)
         μ1          μ2

Notice that the solution in normal coordinates is a linear combination of the modal solutions. The terms φ√ijμ are just scaling factors that represent the contribution of each modal solution to the final solution. This completes modal analysis

1.1.8  Numerical solution using modal analysis

This is a numerical example that implements the above steps using a numerical values for [K]  and [M ]  . Let k1 = 1,k2 = 2,m1 = 1,m2 = 3  and let f1(t) = 0  and f2 (t) = sin (5t)  . Let initial conditions be            ′                     ′
x1 (0) = 0,x1 (0) = 1,x2(0) = 1.5,x 2(0) = 3  , hence

(     )    (  )
{     }    {  }
 x1 (0)  =   0
(x2 (0))    ( 1)

and

({  ′   )}   ({    )}
  x1(0)  =   1.5
( x′(0))   (  3 )
   2

In normal coordinates, the EOM are

pict

In this example m11 = 1,m12  = 0,m21 = 0,m22 = 3  and k11 = 3,k12 = − 2,k21 = − 2,k22 = 2  and f1 (t) = 0  and f2(t) = sin(5t)

step 2 is now applied which solves the eigenvalue problem in order to find the two natural frequencies

pict

Let ω2 = λ  hence

3 λ2 − 11λ + 2 = 0

The solution is λ1 = 3.475  and λ2 = 0.192  , therefore

ω  = √3.475-= 1.864
 1

And

     √ -----
ω2 =   0.192 = 0.438

step 3 is now applied which finds the non-normalized eigenvectors. For each natural frequency ω1   and ω2   the corresponding shape function is found by solving the following two sets of equations for the eigen vectors φ1,φ2

(⌊        ⌋      ⌊    ⌋)  (   )    (  )
    3  − 2     2  1  0    {φ11}    { 0}
(⌈        ⌉ − ω1 ⌈    ⌉)  (   )  = (  )
   − 2  2         0  3     φ21       0

For ω1 = 1.864

pict

This gives one equation to solve for φ21   (the first row equation is only used)

− 0.475− 2φ21 = 0

Hence

      0.475-
φ21 =  − 2 = − 0.237

The first eigen vector is

     (    )    (       )
     { φ11}    {    1  }
φ1 = (    )  = (       )
       φ21       − 0.237

Similarly for ω2 = 0.438

pict

This gives one equation to solve for φ22   (the first row equation is only used)

2.808− 2φ22 = 0

Hence

φ  =  − 2.808-= 1.404
 22     − 2

The second eigen vector is

     (    )    (     )
     { φ12}    {   1 }
φ2 = (    )  = (     )
       φ22       1.404

Now step 4 is applied, which is mass normalization of the shape vectors (or the eigenvectors)

      T
μ1 = φ 1 [M ]φ1

Hence

pict

Similarly, μ2   is found

μ2 = φT2 [M ]φ2

Hence

pict

Now that μ1,μ2   are found, the mass normalized eigen vectors are found. They are called Φ1, Φ2

pict

Similarly

pict

Therefore, the modal transformation matrix is

pict

This result can be verified using Matlab's eig function as follows

 
phi = 
-0.3803   -0.9249 
-0.5340    0.2196
 
0.4380 
1.8641

Matlab result agrees with the result obtained above. The sign difference is not important.

Now step 5 is applied. Matlab generates mass normalized eigenvectors by default.

Now that [Φ]  is found, the transformation from the normal coordinates {x} to modal coordinates, called {η},  is obtained

pict

The transformation from modal coordinates back to normal coordinates is

pict

However,   − 1     T
[Φ ]  = [Φ]  [M ]  therefore

pict

The next step is to apply this transformation to the original equations of motion in order to decouple them.

Applying step 6 results in

pict

The EOM are now decoupled and each EOM can be solved easily as follows

pict

To solve these EOM's, the initial conditions in normal coordinates must be transformed to modal coordinates using the above transformation rules

pict

and

pict

Each of these EOM are solved using any of the standard methods. This results in solutions η1 (t)  and η2 (t).  Hence the following EOM's are solved

pict

and also

pict

The solutions η1 (t)  , η2 (t)  are found using basic methods shown in other parts of these notes. The last step is to transform back to normal coordinates by applying step 7

pict

Hence

x1(t) = 0.925η1(t)+ 0.38η2(t)

and

x2 (t) = 0.534 η1(t) − 0.219 η2(t)

The above shows that the solution x1 (t)  and x2 (t)  has contributions from both nodal solutions.

1.2  Fourier series representation of a periodic function

Given a periodic function f (t)  with period T  then its Fourier series approximation ˜f (t)  using N  terms is

pict

Where

pict

Another way to write the above is to use the classical representation using cos  and sin  . The same coefficients (i.e. the same series) will result.

pict

Just watch out in the above, that we divide by the full period when finding a0   and divide by half the period for all the other coefficients. In the end, when we find f˜(t)  we can convert that to complex form. The complex form seems easier to use.

1.3  Generating Transfer functions for different vibration systems

pict

1.3.1  Force transmissibility

Let steady state

         {             }
           ˆF-        iϖt
xss = Re   k D (r,ζ) e

Then

pict

Hence

            | |   ∘  ---------  | |    ∘ ----------
|f  (t)|    = ||Fˆ|||D |  1+ c2ϖ2- = ||ˆF|||D |  1+ (2ζr)2
 tr   max                 k2

So TR or force transmissibility is

      |ftr (t)|         ∘ ----------
T R = ---||--||max=  |D |  1 + (2 ζr)2
         |Fˆ|

If r > √2--  then we want small ζ  to reduce force transmitted to base. For r < √2--  , it is the other way round.

1.3.2  vibration isolation

We need transfer function between y  and z.  Equation of motion

pict

Let        {  iωt}  ′     {     iωt}
z = Re  Ze    ,z =  Re  iωZe and let       {   iωt}  ′     {     iωt}  ′′     {   2   iωt}
y = Re  Ye    ,y =  Re  iωY e    ,y  = Re  − ω Y e , hence the above becomes

pict

Hence              ∘1+-(2ζr)2-
|D  (r,ζ)| = ∘-----2-----2-
            (1−r2)+(2ζr)   and                             (     )
arg (D) = tan−1(2ζr) − tan−1  21ζ−rr2 where r = ωωn  .

Hence for good vibration isolation we need |Y|max
 |Z| to be small. i.e.     ∘ ----------
|D |  1+ (2ζr)2   to be small. This is the same TR as for force isolation above.

For small |D | , we need small ζ  and small k  (the small k  is to make     √--
r >  2  ) see plot

pict

In Matlab, the above can be plotted using

1.3.3  Accelerometer

We need transfer function between u  and za  where now za  is the amplitude of the ground acceleration. This device is used to measure base acceleration by relating it linearly to relative displacement of m  to base.

Equation of motion. We use relative distance now.

pict

Let z′′ = Re {Z eiωt}.
          a  Notice we here jumped right away to the z′′ itself and wrote it as Re {Z  eiωt}
      a and we did not go through the steps as above starting from base motion. This is because we want the transfer function between relative motion u  and acceleration of base.

Now,       {      }        {        }         {         }
u = Re  Ueiωt ,u′ = Re  iωU eiωt ,u′′ = Re − ω2U eiωt , hence the above becomes

pict

Hence |D  (r,ζ)| = ∘------−12-------2
            (ω2n−ω2)+ (2ωζωn)   and                        (      )
arg(D ) = − 1800 − tan−1 2ωω2ζ−ωωn2-
                          n

When system is very stiff, which means ωn  very large compared to ω  , then D (r,ζ) ≈ −ω2n1Za  , hence by measuring u  we estimate Za  the amplitude of the ground acceleration since ω2
  n  is known. For accuracy, need ω  > 5ω
 n  at least.

1.3.4  Seismometer

Now we need to measure the base motion (not base acceleration like above). But we still use the relative displacement. Now the transfer function is between u  and z  where now z  is the base motion amplitude.

Equation of motion. We use relative distance now.

pict

Let        {     }         {       }
z = Re  Zeiωt  ,z ′ = Re iωZeiωt ,       {         }
z′′ = Re  − ω2Zei ωt ,  and let        {     }         {       }         {         }
u = Re  U eiωt  ,u′ = Re iωU eiωt  ,u′′ = Re  − ω2U eiωt , hence the above becomes

Now, u = Re{U eiωt},u′ = Re{iωU eiωt},u′′ = Re {− ω2U eiωt} , hence the above becomes

pict

Hence                 2
|D  (r,ζ)| = ∘(1−rr2)+i2ζr  and                  (    )
arg(D ) = − tan −1  21−ζrr2-

Now if r  is very large, which happens when ωn ≪ ω  , then ----1-----   -1-
(1−r2)+i2ζr ⇒ −r2   since  2
r   is the dominant factor. Therefore          2
U =  (1−-r2r)+i2ζrZa  now becomes U ≃  − Za  therefore measuring the relative displacement U  gives linear estimate of the ground motion. However, this device requires that ωn  be much smaller than ω  , which means that m  has to be massive. So this device is heavy compared to accelerometer.

1.3.5  Summary of vibration transfer functions

For good isolation of mass from ground motion, rule of thumb: Make damping low, and stiffness low (soft spring).

Isolate base from force. transmitted by machine

Equation used ftr(t) = fspring + fdamper

Transfer function                 ----------
|ftr(t|)||max      ∘         2
   ||Fˆ||   =  |D |  1+ (2ζr)

Isolate machine from motion of base

Equation used. Use absolute mass position

my ′′ + cy′ + ky = cz ′ + kz

Transfer function

            ∘  ----------
|Y-|max = |D |  1+ (2ζr)2
  |Z |

Accelerometer: Measure base acc. using relative displacement

Equation used. Use relative mass position

mu ′′ + cu′ + ku = − mz ′′

Transfer function

pict
Seismometer: Measure base motion using relative displacement

Equation used. Use relative mass position

mu ′′ + cu′ + ku = − mz ′′

Transfer function

pict

1.4  Solution of Vibration equation of motion for different loading

1.4.1  common definitions

These definitions are used throughout the derivations below.

pict

1.4.2  Harmonic loading    ′′    ′
mu  + cu + ku =  F sinϖt

Undamped Harmonic loading

pict

mu ′′ + ku = F sin ϖt

Since there is no damping in the system, then there is no steady state solution. In other words, the particular solution is not the same as the steady state solution in this case. We need to find the particular solution using method on undetermined coefficients.

Let u =  uh + up.  By guessing that up = c1sinϖt  then we find the solution to be

                          F---1---
u = A cosωnt + B sin ωnt+  k 1− r2 sin ϖt

Applying initial conditions is always done on the full solution. Applying initial conditions gives

u (0) = A

pict

Where r = ϖ-
    ωn

The complete solution is

                    (  ′             )
u(t) = u(0)cosωnt +   u-(0)− F- --r--- sinωnt + F---1---sinϖt
                       ωn     k 1− r2           k 1 − r2
(1)

Example: Given force f (t) = 3sin(5t)  then ϖ  = 5  rad/sec, and  ˆ
F = 3  . Let m =  1,k = 1  , then ωn = 1  rad/sec. Hence r = 5  , Let initial conditions be zero, then

pict

Resonance forced vibration When ϖ  ≈ ω  we obtain resonance since r →  1  in the solution given in Eq (1) above and as written the solution can not be used for analysis. To obtain a solution for resonance some calculus is needed. Eq (1) is written as

                   ( u′(0)  F    ωϖ   )         F    ω2
u(t) = u(0)cosωt +   -ω---− -kω2-−-ϖ2-  sinωt + k-ω2-−-ϖ2-sinϖt
(1A)

When ϖ  ≈ ω  but less than ω  , letting

ω − ϖ =  2Δ
(2)

where Δ  is very small positive quantity. And since ϖ  ≈ ω  let

ω + ϖ ≈  2ϖ
(3)

Multiplying Eq (2) and (3) gives

ω2 − ϖ2 =  4Δϖ
(4)

Eq (1A) can now be written in terms of Eqs (2,3) as

pict

Since ϖ ≈ ω  the above becomes

pict

Using                     (     )   (     )
sinϖt −  sinωt = 2 sin  ϖ−2ω-t cos ϖ+2ωt the above becomes

                   u ′(0)        F  ω  (   ( ϖ − ω  )    ( ϖ + ω ) )
u (t) = u(0)cosωt + ----- sin ωt+ -- --- sin  ------t  cos  ------t
                     ω           k 2Δ         2             2

From Eqs (2,3) the above can be written as

                   u′(0)        F--ω-
u(t) = u(0)cosωt +   ω  sinωt + k 2Δ  (sin(− Δt)cos (ϖt ))

Since lim Δ→0 sin(ΔΔt)= t  the above becomes

                   u′(0)        F ωt
u(t) = u (0)cosωt +-----sinωt − -----cos(ωt)
                     ω          k  2

This is the solution to use for resonance.

Underdamped harmonic loading c < cr,ξ < 1

pict

mu ′′ + cu′ + ku = F sinϖt

  ′′       ′   2    F
u  + 2ξωu  + ω u = m- sin ϖt

The solution is

u(t) = uh + up

where

         −ξωt
uh (t) = e    (A cosωdt+ B sinωdt)

and

up (t) = ∘--------F----------sin(ϖt − θ)
          (k − m ϖ)2 + (cϖ )2

where

        ---cϖ-----  -2-ξr-
tan θ = k − mϖ2  = 1 − r2

Very important note here in the calculations of tanθ  above, one should be careful on the sign of the denominator. When the forcing frequency ϖ >  ω  the denominator will become negative (the case of ϖ  = ω  is resonance and is handled separately). Therefore, one should use arctan  that takes care of which quadrant the angle is. For example, in Mathematica use

ArcTan[1 - r^2, 2 Zeta r]]

and in Matlab use

atan2(2 Zeta r,1 - r^2)

Otherwise, wrong solution will result when ϖ >  ω  The full solution is

u(t) = e− ξωt(A cosωdt + B sin ωdt)+ F-∘-------1---------sin(ϖt − θ)
                                   k  (1 − r2)2 + (2ξr)2
(1)

Applying initial conditions gives

pict

Another form of these equations is given as follows

     p0 -------1--------((     2)                 )
up =  k (1 − r2)2 + (2ζr)2 1 − r  sinϖt − 2ζr cosϖt

Hence the full solution is

                                    F        1         ((      )                 )
u (t) = e−ξωnt(Acos ωdt+ B sinωdt)+  ---------2-------2- 1 − r2 sinϖt − 2ζr cosϖt
                                    k (1− r2) + (2ζr)
(1.1)

Applying initial conditions now gives

pict

The above 2 sets of equations are equivalent. One uses the phase angle explicitly and the second ones do not. Also, the above assume the force is F sin ϖt  and not F cosϖt  . If the force is F cosϖt  then in Eq 1.1 above, the term reverse places as in

                                    F        1         ((      )                 )
u (t) = e−ξωnt(Acos ωdt+ B sinωdt)+  ---------2-------2- 1 − r2 cosϖt − 2 ζrsin ϖt
                                    k (1− r2) + (2ζr)

Applying initial conditions now gives

pict

When a system is damped, the problem with the divide by zero when r = 1  does not occur here as was the case with undamped system, since when when ϖ  ≈ ω  or r = 1  , the solution in Eq (1) becomes

            ((                )          (                                          )       )
u(t) = e− ξωt   u(0)+  F-1-sinθ  cosω t +   u′(0)-+ u-(0)ξω + F---1--(ξω sin θ − ϖ cosθ ) sinω t
                      k 2ξ           d      ωd       ωd     k 2ωd ξ                       d
                                                                         F 1
                                                                       + ----sin (ϖt − θ)
                                                                          k2
and the problem with the denominator going to zero does not show up here. The amplitude when steady state response is maximum can be found as follows. The amplitude of steady state motion is Fk-∘----212-----2-
    (1−r )+(2ξr)   . This is maximum when the magnification factor β =  ∘----212-----2
      (1−r) +(2ξr)   is maximum or when ∘ -----------------
        22        2
  (1 − r ) + (2ξr) or ∘ ---------------------
  (    (ϖ-)2)2  (   ϖ)2
   1 −  ω     +  2ξ ω   is minimum. Taking derivative w.r.t. ϖ  and equating the result to zero and solving for ϖ  gives
      ∘ ------2
ϖ  = ω  1 − 2ξ

We are looking for positive ϖ  , hence when ϖ  = ω∘1--−-2ξ2   the under-damped response is maximum.

critically damping harmonic loading     c-
ξ = cr = 1

The solution is

u(t) = uh + up

Where u =  (A + Bt )e−ωt
 h  and u  = F-∘-----1------sin (ϖt −  θ)
 p   k   (1−r2)2+ (2r)2  where tanθ = -2r-
       1−r2   (making sure to use correct arctan  definition). Hence

                      F          1
u (t) = (A + Bt) e−ωt +-∘-----------------sin(ϖt − θ)
                       k   (1 − r2)2 + (2r)2

where A,B  are found from initial conditions

pict

overdamped harmonic loading ξ = ccr-> 1

The solution is

u(t) = uh + up

where

          p1t    p2t
uh(t) = Ae   + Be

and

u (t) = F-∘--------1---------sin (ϖt −  θ)
 p      k   (1− r2)2 + (2ξr)2

hence

       p1t    p2t  F- --------1----------
u = Ae    + Be   +  k ∘ -----2-2-------2 sin(ϖt − θ)
                        (1 − r ) + (2ξr)

where tan θ = 21ξ−rr2   and

pict

Hence the solution is

pict

Solution using frequency approach to harmonic loading
pict

Let load be harmonic and represented in general as    (     )
Re  ˆF eiϖt where Fˆ  is the complex amplitude of the force.

Hence system is represented by

pict

Let        (     )
y = Re  ˆY eiϖt Hence        (        )         (          )
y′ = Re iϖ ˆYeiϖt  ,y′′ = Re − ϖ2 ˆYeiϖt , therefore the differential equation becomes

pict

Dividing numerator and denominator ω2n  gives

      ˆ
Yˆ = F-------1-------
     k (1 − r2)+ i2ζr

Where r = ϖωn,  hence the response is

       ( ˆ                  )
y = Re  F- ------1------eiϖt
         k (1 − r2)+ i2ζr

Therefore, the phase of the response is

            (  )        (        )
arg(y) = arg ˆF  − tan− 1  --2ζr--- + ϖt
                          (1− r2)

Hence at t = 0  the phase of the response will be

            (  )      −1 (  2ζr   )
arg (y) = arg  ˆF  − tan    (1-−-r2)

So when ˆ
F  is real, the phase of the response is simply      − 1(-2ζr-)
− tan    (1−r2)

Undamped case

When ζ = 0  the above becomes

pict

For real force this becomes

    F----1----
y = k (1 − r2) cos(ϖt)

The magnitude |  |
||ˆY || = F--1--
       k(1−r2)   and phase zero.

damped cases

ζ > 0

       (                    )
        Fˆ       1
y = Re  -- -----2-------eiϖt
         k (1 − r )+ i2ζr

pict

Hence for real force and at t = 0  the phase of displacement is

        (      )
− tan− 1  -2ζr--
          1− r2

lag behind the load.

When r < 1  then ϕ  goes from 0  to − 900   Therefore phase of displacement is 0  to − 900   behind force. The minus sign at the front was added since the complex number is in the denominator. Hence the response will always be lagging in phase relative for load.

For r > 1

Now      2
1 − r   is negative, hence the phase will be from     ∘
− 90 to      ∘
− 180

When r = 1

       ( ˆ        )
y = Re   F--1-eiϖt
         k i2ζ

pict

Now phase is −  90∘

pict

Examples. System has ζ = 0.1  and m =  1,k = 1  subjected for force 3cos (0.5t)  find the steady state solution.

Answer          (      )
y(t) = Re Yˆeiϖt ,      ∘ --
ωn =   k-=  1
       m  rad/sec, hence r = 0.5  under the response is

pict

pict

The equation of motion can also be written as  ′′       ′   2    F-
u  + 2ζωu  + ω u = m sinϖt  .

The following table gives the solutions for initial conditions are u(0)  and  ′
u (0)  under all damping conditions. The roots shown are the roots of the quadratic characteristic equation   2           2
λ  + 2ζω λ+ ω  λ = 0  . Special handling is needed to obtain the solution of the differential equation for the case of ζ = 0  and ϖ  = ω  as described in the detailed section below.

Summary table



ζ = 0
roots (
{− iω

(+i ω
u (t)  (|                        ′
{ϖ  = ω →  u(0)cosϖt +  u(ϖ0)sin ϖt − FK-ϖ2tcos (ϖt )
|                      (u′(0)   F  r  )        F   1
(ϖ  ⁄= ω →  u(0)cosωt +  --ω- − K-1−r2 sinωt + K-1−r2 sinϖt


ζ < 1
roots ({          ∘ -----2
 − ξω + iωn  1 − ξ
(          ∘ -----2
 − ξω − iωn  1 − ξ
u (t) = e−ξωt(A cosωdt+ B sinωdt) + FK-∘----212-----2 sin(ϖt − θ)
                                      (1−r )+(2ξr)
          F      1
A =  u0 + K-∘-(1−-r2)2+(2ξr)2 sin θ
     v0   u0ξω-  F--------1-------
B =  ωd +  ωd + K ω ∘ (1−r2)2+-(2ξr)2-(ξω sinθ − ϖ cosθ)
                   d


ζ = 1
roots (
{− ω
(
 − ω
u (t) = (A + Bt)e− ωt + F-∘----1------sin (ϖt − θ)
                      K   (1−r2)2+ (2r)2
A =  u0 + KF∘----12-----sinθ
             (1− r2) +(2r)2
B =  v0 + u0ω + ∘---F∕k2-----(ω sin θ − ϖ cosθ )
                 (1− r2) +(2r)2


ζ > 1
roots (
{− ω ξ + ω ∘ ξ2-−-1
    n     n∘ ------
(− ωnξ − ωn  ξ2 − 1
         (    ∘-2--)        (   ∘ 2--)
u (t) = Ae − ξ+  ξ −1 ωnt + Be −ξ−  ξ− 1ωnt + FKβ sin (ϖt − θ)
                        ((        )          )
     v0+u0ωξ+u0ω∘ ξ2−-1+FKβ  ξ+∘ ξ2−1 ωsinθ−ϖ cosθ
A =  -----------------2ω∘ξ2−1------------------
       v0+u0ωξ−u0ω∘ ξ2−1+ Fβ((ξ−∘ ξ2−-1)ω sinθ−ϖ cosθ)
B =  − -----------------K2ω-∘ξ2−1------------------


1.4.3  constant loading mu ′′ + cu′ + ku = F

Undamped Constant loading case ζ = 0
mu ′′ + ku = F

  ′′   2
u  + ω u = F

u(t) = uh + up

Where uh = A cosωt + B sin ωt  and up = F-
     k  , the solution is

                          F
u(t) = A cosωt + B sin ωt +--
                          k

Applying initial conditions gives

pict

And complete solution is

           (         )
       F-           F-         u-′(0)
u(t) = k +   u(0)−  k  cosωt +   ω   sin ωt

underdamped constant loading ζ < 1

The general solution is

                                   F
u (t) = e−ξωt(A cosωdt + B sin ωdt)+--
                                   k

From initial conditions

pict

Hence the solution is

            (                     (                      )       )
        −ξωt (        F )           u′(0)+ u (0)ξω − F-ξω             F
u(t) = e       u(0) − k-  cosωdt+   ----------ω------k---  sinωdt  +  k-
                                               d

Critical damping constant loading ζ = 1

The general solution is

                −ωt  F-
u(t) = (A + Bt )e  +  k

Where from initial conditions

pict

Over-damped constant loading ζ > 0

The solution is

         p1t     p2t   F-
u(t) = Ae   + Be   +  k

Where now

pict

Hence the solution is

u(t) = Aep1t + Bep2t + F
                      k

Where

pict

Summary table for constant loading solutions



ζ = 0
roots({
  − iω
( +iω
       (     F)         v         F
u (t) =  u0 − k- cos ωt+  0ω-sin ωt+  k-


ζ < 1
roots(
{ − ξω + iω ∘1-−-ξ2
          n∘ ------
( − ξω − iωn 1 − ξ2
            ( (      )         (        F-  )       )
u (t) = e−ξωt   u0 − Fk cosωdt +  v0+u0ξωω−-kξω- sin ωdt  + Fk-
                                      d


ζ = 1
roots (
{− ω

(− ω
u (t) = ((u −  F) + (v + u ω − F-ω) t) e−ωt + F
          0   k      0   0     k            k


ζ > 1
roots (          ∘ ------
{− ωnξ + ωn  ξ2 − 1
(          ∘ ------
 − ωnξ − ωn  ξ2 − 1
B =  Fkp1−u0p1+v0
       (p2−p1)
         F
A =  u0 −-k − B
            ∘ ----------               ------
p1 = − -c-+   (c--)2 − k-= − ωn ξ + ωn ∘ ξ2 − 1
       2m      2m     m
            ∘ (---)2----             ∘ ------
p2 = − 2cm-−    2cm-  − km-=  − ωn ξ − ωn ξ2 − 1
u (t) = Aep1t + Bep2t + F
                      k


1.4.4  No loading (free vibration) mu ′′ + cu ′ + ku = 0

Undamped free vibration

pict

   ′′
mu   + ku = 0

u′′ + ω2u = 0

The solution is

                   u′(0)
u(t) = u(0)cosωt + -----sinωt
                     ω

under-damped free vibration c < cr,ξ < 1

pict

mu ′′ + cu′ + ku = 0

u′′ + 2ξωu′ + ω2u = 0

The solution is

u = e−ξωt(A cosωdt+  B sin ωdt)

Applying initial conditions gives A = u (0)  and      ′
B = u-(0)+ωud(0)ξω  . Therefore the solution becomes

            (              u′(0)+ u (0)ξω       )
u (t) = e−ξωt  u(0)cosωdt + ---------------sin ωdt
                                 ωd

critically damped free vibration     c-
ξ = cr = 1

The solution is

pict

where A, B  are found from initial conditions A = u (0)  ,      ′
B  = u (0)+ u (0)ω  , hence

u(t) = (u (0)+ (u′(0)+ u (0)ω) t) e−ωt

over-damped free vibration ξ = ccr-> 1

The solution is

         λ t     λ t
u (t) = Ae 1 + Be  2

where A,B  are found from initial conditions.

pict

where λ1   and λ2   are the roots of the characteristic equation

pict

Summary table for free vibration solutions



ζ = 0
 ′′   2
u  + ω u = 0
roots(
{ − iω

( +iω
                   u′(0)
u (t) = u(0)cosωt +   ω sinωt
(
||| u(t) = A cos(ωt − ϕ)
||{     ∘ ------(-----)2
  A =   u(0)+   u′(0)
|||          (     ω)
||( ϕ = tan −1  u′(0)∕ω-
              u(0)


ζ < 1
roots(         ∘  ------
{ − ξω + iω  1− ξ2
          ∘  ------
( − ξω − iω  1− ξ2
        −ξωt(             u′(0)+u(0)ξω       )
u (t) = e     u (0) cosωdt+     ωd     sin ωdt


ζ = 1
roots(
{ − ω

( − ω
                       ′     −ωt
u (t) = (u(0)(1+ ωt) + u(0)t)e


ζ > 1
roots(              ∘ ------
{ λ1 = − ω ξ + ω ξ2 − 1
(              ∘ ------
  λ2 = − ω ξ − ω ξ2 − 1
(|          λ1ωt    λ2ωt
|||{ u(t) = Ae    + Be
      u′(0)∘−u(0)λ2
|| A =  2ω  ξ2− 1
||( B = −-u′(0)∘+u(0)λ1
        2ω  ξ2−1


Roots of characteristic equation

The roots of the characteristic equation for   ′′       ′   2
u  + 2ξωu  + ω u = 0  are given in this table




roots time constant τ



ξ < 1  {          ∘ ------          ∘ ------}
 − ξω + jωn  1 − ξ2,− ξω − iωn  1 − ξ2 ξ1ω



ξ = 1  {− ω, − ω } -1
ω



ξ > 1  {          ∘ -2----          ∘ -2----}
 − ωnξ + ωn  ξ − 1,− ωnξ − ωn  ξ  − 1 ------1----- -----1------
ωnξ−ωn∘ ξ2−1,ωnξ+ωn∘ ξ2−-1   (which to use? the bigger?)



1.4.5  impulse F0 δ(t)  loading

impulse input

Undamped system with impulse

m¨u + ku = F0δ(t)

with initial conditions u (0) = 0  and u′(0) = 0.  Assuming the impulse acts for a very short time period from 0  to t1   seconds, where t1   is small amount. Integrating the above differential equation gives

∫ t1       ∫ t1       ∫ t1
    m ¨udt+     kudt =     F0δ(t)
 0          0          0

Since t1   is very small, it can be assumed that u  changes is negligible, hence the above reduces to

pict

since we assumed u′(0) = 0  and since ∫
 t1δ(t) = 1
 0  then the above reduces to

       F
˙u(t1) =--0
        m

Therefore, the effect of the impulse is the same as if the system was a free system but with initial velocity given by F0
m  and zero initial position. Hence the system is now solved as follows

m ¨u+ ku =  0

With u (0) = 0  and u′(0) = Fm0  . The solution is

             F0
uimpulse(t) = ---sinωt
             mω

If the initial conditions were not zero, then the solution for these are added to the above. From earlier, it was found that the solution is                    u′(0)-
u(t) = u (0) cos ωt+  ω  sin ωt  , therefore, the full solution is

      ◜-----due to I◞C◟-only---◝  du◜e to i◞◟mpuls◝e
                   u′(0)          F0
u(t) = u (0) cosωt+ -----sin ωt+  ----sin ωt
                    ω           m ω

under-damped with impulse c < cr,ξ < 1

m ¨u + c˙u+  ku = δ(t)

u¨+ 2ξω ˙u+ ω2u =  δ(t)

with initial conditions u(0) = 0  and u′(0) = 0.  Integrating gives

∫ t1        ∫ t1      ∫ t1       ∫ t1
    mu¨dt +     c˙udt+     kudt =     F0δ(t)
 0          0          0         0

Since t
 1   is very small, it can be assumed that u  changes is negligible as well as the change in velocity, hence the above reduces to the same result as in the case of undamped. Therefore, the system is solved as free system, but with initial velocity  ′
u (0) =  F0∕m  and zero initial position.

Initial conditions are u (0) = 0  and u′(0) = 0  then the solution is

           −ξωt
uimpulse = e    (A cosωdt+ B  sinωdt)

applying initial conditions gives A = 0  and      (F0)
B  = --m--
      ωd  , hence

                  (           )
uimpulse(t) = e− ξωt -F0--sin ωdt
                   m ωd

If the initial conditions were not zero, then the solution for these are added to the above. From earlier, it was found that the solution is             (                             )
u (t) = e−ξωt u(0)cosωdt + u′(0)+u(0)ξωsinωdt
                              ωd , therefore, the full solution is

                      due to IC only                   due to impulse
       ◜----(-------------◞◟′-----------------)◝  ◜----(---◞◟-------◝)
u (t) = e−ξωt u (0 )cosωdt+  u(0)+--u(0)ξω-sinωdt  + e− ξωt  -F0--sin ωdt
                                ωd                       m ωd

critically damped with impulse input     c-
ξ = cr = 1  with initial conditions u (0) = 0  and  ′
u (0) = 0  then the solution is

pict

where A, B  are found from initial conditions A = u (0) = 0  and       ′              F0
B =  u (0)+ u(0)ω =  m-  , hence the solution is

            F0t  −ωt
uimpulse(t) =-m- e

If the initial conditions were not zero, then the solution for these are added to the above. From earlier, it was found that the solution is u(t) = (u (1+ ωt)+ u′(0)t)e−ωt
        0  , therefore, the full solution is

              due to IC only        due◜ to i◞m◟pu◝lse
      ◜(------------◞◟--′----)-−ω◝t    F0t − ωt
u(t) = u (0)(1+ ωt) + u (0)t e   +   -m-e

over-damped with impulse input      c
ξ = cr > 1  With initial conditions are u (0) = 0  and u ′(0) = 0  the solution is

u      (t) = Aeλ1ωt + Be λ2ωt
 impulse

where A,B  are found from initial conditions and

pict

Hence the solution is

               (    ∘----)       (   ∘ ---)
u       (t) = Ae − ξ+  ξ2−1 ωt + Be −ξ−  ξ2− 1ωt
  impulse

where

pict

Hence

                       (         )                (        )
            ----F0m----- − ξ+ ∘ξ2−1 ωt  -----Fm0----  −ξ−∘ ξ2−1-ωt
uimpulse(t) = 2ω ∘ ξ2 −-1e            − 2 ω∘ ξ2 −-1-e

If the initial conditions were not zero, then the solution for these are added to the above. From earlier, it was found that the solution is          p t     pt
u (t) = Ae 1 + Be  2  , therefore, the full solution is

                            F0               F0
u(t) = Ae λ1ωt + Beλ2ωt +-∘-m-----eλ1ωt −---∘m-----eλ2ωt
                        2ω  ξ2 − 1       2ω   ξ2 − 1

pict

Summary table



ζ = 0
 ′′   2
u  + ω u = 0
roots(
{ − iω

( +iω
               transient          steady state
       ◜---------◞◟-′-------◝   ◜--◞◟---◝
u (t) = u(0)cosωt + u-(0)sinωt + -F0-sin ωt
                     ω          m ω


ζ < 1
roots(         ∘  ------
{ − ξω + iω  1− ξ2
          ∘  ------
( − ξω − iω  1− ξ2
        ----------------transient----------------    -----steady state----
       ◜    (              ◞◟′                 )◝   ◜    (  ◞◟       ◝)
u (t) = e−ξωt  u(0)cosωdt+  u-(0)-+-u(0)ξω-sinωdt  +  e− ξωt -F0--sin ωdt
                                ωd                       m ωd


ζ = 1
roots(
{ − ω

( − ω
                       ′      − ωt  F0t −ωt
u (t) = (u(0)(1 + ωt)+ u (0)t)e    +  m e


ζ > 1
roots(              ∘ ------
{ λ1 = − ω ξ + ω ξ2 − 1
(              ∘ ------
  λ2 = − ω ξ − ω ξ2 − 1
(|          λ1ωt    λ2ωt   --F0m----λ1ωt  ---Fm0--- λ2ωt
||| u(t) = Ae    + Be    +  2ω∘ξ2−1e    − 2ω∘ ξ2−-1e
{     u′(0)−u(0)λ2
|| A =  2ω∘ ξ2−-1
||(     −-u′(0)∘+u(0)λ1
  B =   2ω  ξ2−1


The impulse response can be implemented in Mathematica as

pict

Impulse sin function

Now assume the input is as follows

pict

given by F (t) = F0sin(ϖt )  where ϖ = -2π-=  π-
    2t1   t1

undamped system with sin impulse

           (
           {  F0sin(ϖt)    0 ≤ t ≤ t1
m ¨u + ku = (
                  0          t > t1

with u (0) = u0   and u˙(0) = v0.  For 0 ≤ t ≤ t1   the solution is

                 ( v0        r  )              1     ( π  )
u (t) = u0cosωt +   --− ust-----2  sin ωt+  ust-----2 sin --t
                   ω      1 − r             1 − r      t1

where     ϖ-   π∕t1-  -T-
r =  ω =  ω  = 2t1   where T  is the natural period of the system.       F0
ust = k  , hence the above becomes

                 (           ( π∕t )   )
                 | v0   F0     -ω1-   |         F0     1         (  t )
u (t) = u0cos ωt+ ( -- − -------(---)2-) sin ωt+  ------(----)2-sin  π --
                    ω    k 1−   π∕ωt1-            k  1−   π∕tω1-        t1
(1)

When u0 = 0  and v0 = 0  then

pict

The above Eq (1) gives solution during the time 0 ≤ t ≤ t1

Now after t = t1   the force will disappear, the differential equation becomes

m ¨u + ku = 0    t > t
                    1

but with the initial conditions evaluate at t = t.
     1  From (1)

pict

since sin ϖt1 = 0  . taking derivative of Eq (1)

                     ( v0        r  )             1
˙u(t) = − ωu0 sin ωt+ ω  -- − ust----2  cosωt + ϖ -----2 cosϖt
                       ω      1 − r             1−  r

and at t = t1   the above becomes

pict

since cosϖt1 = − 1.  Now (2) and (3) are used as initial conditions to solve m¨u + ku = 0  . The solution for t > t1   is

u(t) = u (t1)cos ωt+ u˙(t1) sin ωt
                     ω

1.4.5.0.0.0 Resonance with undamped sin impulse When ϖ  ≈ ω  and t ≤ t
     1   we obtain resonance since r → 1  in the solution shown up and as written the solution can't be used for analysis in this case. To obtain a solution for resonance some calculus is needed. Eq (1) is written as

pict

Now looking at case when ϖ ≈ ω  but less than ω  , hence let

ω − ϖ =  2Δ
(2)

where Δ  is very small positive quantity. and we also have

ω + ϖ ≈  2ϖ
(3)

Multiplying Eq (2) and (3) with each others gives

ω2 − ϖ2 =  4Δϖ
(4)

Going back to Eq (1A) and rewriting it as

pict

Since ϖ ≈ ω  the above becomes

pict

now using sinϖt − sinωt = 2 sin (ϖ−ω-t) cos(ϖ+-ωt)
                       2          2 the above becomes

                                  (   (        )    (       ) )
u (t) = u0cos ωt+ v0 sin ωt+  ust-ω-  sin  ϖ-−-ω-t cos  ϖ--+-ωt
                  ω           2Δ          2             2

From Eq(2) ϖ − ω = − 2Δ  and ω + ϖ ≈ 2ϖ  hence the above becomes

                 v0           -ω-
u (t) = u0cosωt + ω  sin ωt+ ust2Δ (sin(− Δt)cos(ϖt ))

or since ϖ ≈  ω

                 v0            ω
u (t) = u0cosωt + -- sin ωt− ust---(sin(Δt )cos(ωt))
                 ω            2Δ

Now lim     sin(Δt)-= t
   Δ →0   Δ  hence the above becomes

u(t) = u0cosωt + v0 sinωt − ustωt-cos(ωt)
                 ω            2

This can also be written as

pict

since ϖ  ≈ ω  in this case. This is the solution to use for resonance and for t ≤ t1

Hence for t > t1   , the above equations is used to determine initial conditions at t = t1

                   v0            ϖt1-
u(t1) = u0 cosϖt1 + ϖ sin ϖt1 − ust 2 cos (ϖt1 )

but cos ϖt1 = cos π-t1 = − 1
            t1  and sinϖt1 =  0  and ϖt1 = π
 2    2   , hence the above becomes

                π
u(t1) = − u0 + ust
                2

Taking derivative of Eq (1) gives

                                   2
˙u(t) = − ϖu sinϖt + v cos ϖt + u  ϖ-t-sin(ϖt )− u  ϖ-cos (ϖt )
           0         0          st 2            st2

and at t = t1

pict

Now the solution for t > t1   is

pict

under-damped with sin impulse c < c ,ξ < 1
     r

                (
                { F0 sin (ϖ)    0 ≤ t ≤ t1
m ¨u + c˙u+ ku =
                (     0          t > t1

or

                 (
                 {  F  sin (ϖ)    0 ≤ t ≤ t
¨u + 2ξωu˙+ ω2u =      0                  1
                 (      0          t > t1

m ¨u + c˙u+  ku = F sin ϖt

                  F
¨u + 2ξωu˙+ ω2u = -- sin ϖt
                 m

For t ≤ t
    1   Initial conditions are u (0) = u
        0   and u˙(0) = v
        0   and u   = F-
  st   k  then the solution from above is

                                          ust
u(t) = e− ξωt(A cosωdt + B sin ωdt)+ ∘-----------------sin(ϖt − θ)
                                    (1 − r2)2 + (2ξr)2
(1)

Applying initial conditions gives

pict

For t > t1   . From (1)

u (t1) = e−ξωt1 (A cosωdt1 + B sinωdt1)+ ∘-----ust--------sin(ϖt1 − θ)
                                         (1− r2)2 + (2ξr)2
(2)

Taking derivative of (1) gives

˙u(t) = − ξωe− ξωt (A cosω t+ B sin ω t)+ e−ξωt(− Aω  sin ω t+ ω  B cos ω t)
                       d         d               d    d     d     ud
                                                      + ϖ ∘--------st--------cos(ϖt − θ)
                                                            (1− r2)2 + (2ξr)2

at t = t1

˙u(t1) = − ξωe− ξωt1 (A cosωdt1 + B sinωdt1) + e−ξωt1 (− A ωdsinωdt1 + ωdB cosωdt1)
                                                                 u
                                                     + ϖ ∘--------st--------cos(ϖt1 − θ)
                                                           (1−  r2)2 + (2ξr)2
(3)

Now for t > t1   the equation becomes

m¨u + c˙u + ku = 0

which has the solution

u = e−ξωt(A cosωdt+  B sin ωdt)

where A  = u(t1)  and B = u˙(t1)+uω(dt1)ξω

critically damped with sin impulse     -c
ξ = cr = 1  For t ≤ t1   Initial conditions are u (0) = u0   and ˙u(0) = v0   then the solution is from above

                 −ωt          ust
u (t) = (A + Bt )e  + ∘---------2-------sin(ϖt − θ)
                         (1 − r2) + (2r)2
(1)

Where         --cϖ---   2ξr-
tan θ = k−mϖ2 =  1−r2.  A, B  are found from initial conditions

pict

For t > t1   the solution is

                                  −ωt
u (t) = (u(t1)+ (˙u (t1)+ u (t1) ω)t)e
(2)

To find u(t1),  from Eq(1)

                 − ωt1   -------ust--------
u (t1) = (A + Bt)e     + ∘ -------2------2-sin (ϖt1 − θ)
                          (1− r2) +  (2r)

taking derivative of (1) gives

˙u(t) = − ω (A + Bt)e− ωt + Be −ωt + ϖ ∘---ust--------sin (ϖt − θ)
                                     (1 − r2)2 + (2r)2
(3)

at t = t1

                                              u
˙u(t1) = − ω (A + Bt1)e−ωt1 + Be −ωt1 + ϖ ∘------st-------sin(ϖt1 − θ)
                                        (1 − r2)2 + (2r)2
(4)

Hence Eq (2) can now be evaluated using Eq(3,4)

over-damped with sin impulse     c-
ξ = cr > 1  For t ≤ t1   Initial conditions are u (0) = u0   and ˙u(0) = v0   then the solution is

u = Aep1t + Bep2t + ∘------ust--------sin(ϖt − θ)
                           2 2       2
                      (1 − r ) + (2ξr)

where         2ξr
tanθ = 1−r2   (make sure you use correct quadrant, see not above on arctan  ) and

pict

and

pict

leading to the solution where         2ξr
tanθ = 1−r2   and

pict

is

pict

For t > t1   . From Eq(1) and at t = t1

          (   ∘ -2--)       (    ∘-2--)
u(t1) = Ae −ξ+  ξ− 1 ωt1 + Be −ξ−  ξ −1 ωt1 + D sin (ϖt1 − θ)
(2)

Taking derivative of Eq (1)

           (   ∘ 2--)         (   ∘-2--)
u˙(t) = ωAe  −ξ+  ξ− 1ωt + ωBe  −ξ− ξ −1 ωt + ϖD cos(ϖt − θ)

At t = t1

           (   ∘ ----)         (   ∘----)
˙u(t ) = ωAe −ξ+  ξ2− 1ωt1 + ωBe  −ξ− ξ2−1 ωt1 + ϖD cos(ϖt  − θ)
   1                                                     1
(3)

Equation of motion now is

¨u+ 2ξω ˙u+  ω2u = 0

which has solution for over-damped given by

          (   ∘----)        (   ∘ ---)
u (t) = Ae  −ξ+ ξ2−1 ωnt + Be −ξ−  ξ2− 1ωnt

where

pict

Input is given by F (t) = F0sin(ϖt)  where      2π-   π-
ϖ  = 2t1 = t1

pict

Summary table



ζ = 0
roots (
{ − iω

( +iω
u(t)  (          (|             u′(0)
||||          { u(0)cosωt + --ω-sinωt − F0k ω2tcos(ωt)  0 ≤ t ≤ t1
|||| ϖ = ω →  | (            )        −u′(0)+F0k π2t--
|{          (  − u (0)+ F0k π2( cosωt +-----ω---)1sinωtt > t1
           (|                ′     (F0∕k)( π∕t1)                     (   )
||||          |{ u(0)cosωt +   u(ϖ0)−  ---(π∕t-ω)2-- sin ωt+  --F(0π∕∕kt)2 sin πtt1 0 ≤ t ≤ t1
||| ϖ ⁄= ω →                          1−  -ω1             1−  -ω1
||(          ||( u(t )cosωt + u′(t1)-sin ωt                                  t > t
                1          ω                                              1


ζ < 1
roots (          ∘ ------
{ − ξω + iωn 1 − ξ2
           ∘ ------
( − ξω − iωn 1 − ξ2   time constant τ = ζ1ωn
      (                                                (       )
      |{ e−ξωt(A cosωdt+  B sin ωdt) + F0k∘----212-----2 sin π tt1-− θ 0 ≤ t ≤ t1
u(t) =       (                          (1− r) +)(2ξr)
      |( e−ξωt u (t1) cosωdt+ u′(t1)+u(t1)ξω sin ωdt                  t > t1
                                 ωd
           F0------1------
A = u(0) + k ∘ (1−r2)2+-(2ξr)2-sin θ
    u′(0)   u(0)ξω   F0-------1-------
B =  ωd +   ωd  +  k ω ∘(1−r2)2+(2ξr)2 (ξω sin θ − ϖ cos θ)
                      d


ζ = 1
roots (
{ − ω
(
  − ω
      (|           −ωt  F0 -----1-------  (  t-   )
      { (A + Bt) e   +  k ∘(1−r2)2+(2r)2 sin πt1 − θ 0 ≤ t ≤ t1
u(t) = |
      ( u (t) = (u(t1)+ (u′(t1) + u(t1)ω)t)e− ωt    t > t1
A = u(0) + F0∘-----1------sin θ
           k   (1−r2)2+ (2r)2
                                   (               )
B = u′(0)+ u (0)ω + F0∘-----1------ ω sin θ − π-cosθ
                    k   (1− r2)2+(2r)2          t1


ζ > 1
roots ({           ∘ -2----
  − ωn ξ + ωn ξ − 1
(           ∘ -2----
  − ωn ξ − ωn ξ − 1
      (    (   ∘----)       (   ∘ ----)          (       )
      |{ Ae  −ξ+ ξ2−1 ωt + Be −ξ−  ξ2− 1 ωt + F-βsin π t-− θ 0 ≤ t ≤ t
u(t) =            (   ∘ ----)        (    ∘k---)    t1             1
      |( u (t) = A1e −ξ+  ξ2−1 ωnt + B1e − ξ− ξ2−1 ωnt       t > t1
     ′             ∘ ---- F ((  ∘ ---)      π    )
A = u-(0)+u(0)ωξ+u(0)ω--ξ2−1+∘kβ--ξ+--ξ2−1-ω-sinθ−t1 cosθ-
                        2ω  ξ2−1
                     ∘----    ((  ∘----)           )
      u′(0)+u-(0)ωξ−-u(0)ω-ξ2−1+∘Fkβ--ξ−-ξ2−1-ωsin-θ−-πt1-cosθ--
B = −                    2ω  ξ2− 1
β = ∘-----12-----2
     (1−r2)+(2ξr)
                 (  ∘ ---)
       ˙u(t1)+u(t1)ω∘n-ξ−--ξ2−1--
A1 = −      2ωn  ξ2− 1
                (  ∘ ---)
     ˙u(t1)+u(t1)ξωn-ξ+--ξ2−-1-
B1 =       2ωn∘ ξ2−1-


1.4.6  Tree view look at the different cases

This tree illustrates the different cases that needs to be considered for the solution of single degree of freedom system with harmonic loading.

There are 12  cases to consider. Resonance needs to be handled as special case when damping is absent due to the singularity in the standard solution when the forcing frequency is the same as the natural frequency. When damping is present, there is no resonance, however, there is what is called practical response which occur when the forcing frequency is almost the same as the natural frequency.

pict

The following is another diagram made sometime ago which contains more useful information and is kept here for reference.

pict

1.4.7  Cycles for the peak to decay by half its original value

This table shows many cycles it takes for the peak to decay by half its original value as a function of the damping ζ  . For example, we see that when ζ = 2.7%  then it takes 4  cycles for the peak (i.e. displacement) to reduce to half its value.

pict

1.4.8  references

  1. Vibration analysis by Robert K. Vierck
  2. Structural dynamics theory and computation, 5th edition by Mario Paz, William Leigh
  3. Dynamic of structures, Ray W. Clough and Joseph Penzien
  4. Theory of vibration,volume 1, by A.A.Shabana
  5. Notes on Diffy Qs, Differential equations for engineers, by Jiri Lebl, online PDF book, chapter 2.6, oct 1,2012 http://www.jirka.org/diffyqs/

2  Dynamics equations, kinematics, velocity and acceleration diagrams

2.1  Derivation of rotation formula

This formula is very important. Will show its derivation now in details. It is how to express vectors in rotating frames.

Consider this diagram

pict

In the above, the small axis x,y  is a frame attached to some body which rotate around this axis with angular velocity ω  (measured by the inertial frame of course). All laws derived below are based on the following one rule

   |            |
d- ||         d- ||
dtr|absolute = dtr|relative + ω × r
(1)

Lets us see how to apply this rule. Let us express the position vector of the particle rp  . We can see by normal vector additions that the position vector of particle is

rp = ro + r
(2)

Notice that nothing special is needed here, since we have not yet looked at rate of change with time. The complexity (i.e. using rule (1)) appears only when we want to look at velocities and accelerations. This is when we need to use the above rule (1). Let us now find the velocity of the particle. From above

˙rp = r˙o + ˙r

Every time we take derivatives, we stop and look. For any vector that originates from the moving frame, we must apply rule (1) to it. That is all. In the above, only r  needs rule (1) applied to it, since that is the only vector measure from the moving frame. Replacing ˙r
 p  by V
 p  and ˙r
 o  by V
 o  , meaning the velocity of P  and o  , Hence the above becomes

Vp = Vo + ˙r

and now we apply rule (1) to expand ˙r

Vp = Vo + (Vrel + ω × r)
(3)

where Vrel  is just d-||
dtr relative

The above is the final expression for the velocity of the particle Vp  using its velocity as measured by the moving frame in order to complete the expression.

So the above says that the absolute velocity of the particle is equal to the absolute velocity of the base of the moving frame + something else and this something else was (Vrel + ω × r)

Now we will find the absolute acceleration of P  . Taking time derivatives of (3) gives

         (                   )
˙Vp = V˙o + V˙rel + ω˙× r + ω × ˙r
(4)

As we said above, each time we take time derivatives, we stop and look for vectors which are based on the moving frame, and apply rule (1) to them. In the above, V˙
  rel  and ˙r  qualify. Apply rule (1) to ˙V
 rel  gives

V˙  = a   + ω × V
 rel   rel       rel
(5)

where arel  just means the acceleration relative to moving frame. And applying rule (1) to r˙  gives

r˙= Vrel + ω × r
(6)

Replacing (5) and (6) into (4) gives

pict

Eq (7) says that the absolute acceleration ap  of P  is the sum of the acceleration of the base ao  of the moving frame plus the relative acceleration arel  of the particle to the moving frame plus 2 (ω × V   )+ (˙ω × r)+ (ω × (ω × r))
        rel

Hence, using Eq(3) and Eq(7) gives us the expressions we wanted for velocity and acceleration.

2.2  Miscellaneous hints

  1. When finding the generalized force for the user with the Lagrangian method (the hardest step), using the virtual work method, if the force (or virtual work by the force) ADDS energy to the system, then make the sign of the force positive otherwise the sign is negative.
  2. For damping force, the sign is always negative.
  3. External forces such as linear forces applied, torque applied, in general, are positive.
  4. Friction force is negative (in general) as friction takes energy from the system like damping.

2.3  Formulas

pict

pict

2.4  Velocity and acceleration diagrams

2.4.1  Spring pendulum

pict

2.4.2  pendulum with blob moving in slot

pict

2.4.3  spring pendulum with block moving in slot

pict

2.4.4  double pendulum

pict

2.5  Velocity and acceleration of rigid body 2D

pict

Finding linear acceleration of center of mass of a rigid body under pure rotation using fixed body coordinates.

In the above U  is the speed of the center of mass in the direction of the x  axis, where this axis is fixed on the body itself. Similarly, V  is the speed of the center of mass in the direction of the y  axis, where the y  axis is attached to the body itself.

Just remember that all these speeds (i.e. U  ,V  ) and accelerations (ax  , ay  ) are still being measured by an observer in the inertial frame. It is only that the directions of the velocity components of the center of mass is along an axis fixed on the body. Only the direction. But actual speed measurements are still done by a stationary observer. Since clearly if the observer was sitting on the body itself, then they will measure the speeds to be zero in that case.

2.6  Velocity and acceleration of rigid body 3D

2.6.1  Using Vehicle dynamics notations

pict

2.6.2  3D Not Using vehicle dynamics notations

pict

Derivation for F = ma  in 3D
pict

Derivation for τ = Iω  in 3D

Let A = Iω  then using the rule

pict

Then τ = Iω   can be found for the general case

pict

Derivation for τ = Iω  in 3D using principle axes

The above derivation simplifies now since we will be using principle axes. In this case, all cross products of moments of inertia vanish.

    (             )
    | Ixx   0    0 |
    ||             ||
I = |(  0  Iyy   0 |)

       0   0   Izz

Hence

pict

So, we can see how much simpler it became when using principle axes. Compare the above to

(             ) (   )    (                                                    )

| Ixx Ixy  Ixz| | αx|    | ωy(Izxωx + Iyzωy + Izzωz)− ωz (Iyxωx + Iyyωy + Iyzωz )|
|| I   I    I  || || α || +  ||ω  (I  ω +  I ω  + I  ω )− ω  (I ω  + I  ω  + I  ω )||
|(  yx   yy   yz|) |(  y|)    |(  x  zx x   yz  y   zz z     z  xx  x   xy y    xz z |)
  Izx  Iyz  Izz    αz      ωx (Iyxωx + Iyyωy + Iyzωz)− ωy (Ixxωx + Ixyωy + Ixzωz)

So, always use principle axes for the body fixed coordinates system!

2.6.3  Acceleration terms due to rotation and acceleration

pict

pict

2.7  Wheel spinning precession

pict

2.8  References

  1. Structural Dynamics 5th edition. Mario Paz, William Leigh

2.9  Misc. items

The Jacobian matrix for a system of differential equations, such as

pict

is given by

    (           )
    | dfdx  ddfy  ddfz|
    || dg  dg  dg||
J = | dx  dy  dz|
    ( dh  dh  dh)
      dx  dy  dz

For example, for the given the following 3 set of coupled differential equations in n3

pict

then the Jacobian matrix is

    (                   )

    ||  0    − 1    − 1  ||
J = |  1    a      0    |
    |(                   |)
      z(t)  0   x (t) − c

Now to find stability of this system, we evaluate this matrix at t = t0   where x(t0),y(t0),z(t0)  is a point in this space (may be stable point or initial conditions, etc...) and then J  become all numerical now. Then we can evaluate the eigenvalues of the resulting matrix and look to see if all eigenvalues are negative. If so, this tells us that the point is a stable point. I.e. the system is stable.

If X is N (0,1 )  distributed then mu + sigma  ∗X  is N (mu, sigma2)  distributed.

3  Astrodynamics

3.1  Ellipse main parameters

pict

3.2  Table of common equations

The following table contains the common relations to use for elliptic motion. Equation of ellipse is  2    2
xa2 + yb2 = 1



term to find

relation



conversion between E  and θ

   (  )   ∘ ------   (   )
tan   θ- =    1+-e-tan  E-
     2       1− e      2
           e-+-cosθ--
   cosE =  1+ ecos θ
           e − cosE
   cosθ =  ----------
           ecosE − 1



position of satellite at time t  Solve for E  , then find θ  . τ  here is time at perigee  and n  is mean satellite speed.

E − esinE = n (t− τ)



eccentricity e

                 --------
    c   ra−-rp-  ∘     2ℰh2
e = a = ra+rp =   1 +  μ2



Maor axes a

    r (1 + e)
a = -p----2--
     1−  e
 =  ra(1-−-e)-
     1 − e2
 = − -μ-
   ∘ 2ℰ------
 =    b2 + c2
      p
 =  ----2-
    1− e



Minor axes b

    ∘ -----2
b = a 1 − e



rp

    a (1 − e2)
rp =---------
      1 + e
  = a(1 − e)



ra

    a (1 − e2)
ra =---------
      1 − e
  = a (1+ e)
    h2   1
  = --------
     μ 1 − e  rrpa = 11−+ee



p

p = a (1− e2) = h2 = r (1 + e) = r (1 − e)
                μ    p          a



specific angular momentum h

h = rpvp = rava = ⃗r × ⃗v = √p-μ
    √---
h =  μr   (circular orbit)



Total Energy ℰ

    v2   μ     μ
ℰ = 2-−  r = −2a



velocity v

         ∘ --(------)-
     v =   μ  2-− 1-   (vis- viva)
              r   a
         ∘ ---
vescape =   2μ-  (escape velocity for parabola)
         ∘ -r
           μ-
 vradial =   p esin θ
         ∘ --
vnormal =   μ-(1+ e cosθ)
           p



vperigee  (closest)

    ∘ --(------)-
       μ- 1-+-e
vp =   a  1 − e
    ∘ --
  =    μ(1+ e)
       p
    ∘ --(-------)-
  =   μ   2-−  1-
          rp   a



vapogee  (furthest)

    ∘ --(------)-
       μ- 1-−-e
va =   a  1 + e
    ∘ --
  =    μ(1−  e)
       p
    ∘ --(-------)-
  =   μ   2-−  1-
          ra   a



magnitude of ⃗r

          (    2)     2
    r = a--1−-e---= h------1-----
        1 + ecosθ    μ 1 + ecosθ
rcosθ = a(cosE −  e)

    r = a(1 − ecosE )  (eq 4.2-14 Bate book)



period T

    2         ∘ a3-
T = hπab = 2π   μ-



mean satellite speed n

         ∘ ---
n = 2Tπ=    μa3



eccentric anomaly E

       ∘ ----
tan θ=    1+e-tan E-
    2     1− e    2



area sweep rate

dA-  h
dt = 2



equation of motion

¨  -μ
⃗r + r3⃗r = 0



spherical coordinates relation

cos(i) = sin(Az )cos(ϕ )  where i  is the inclination and Az  is the azimuth and ϕ  is latitude 1 pict



Notice in the above, that the period T  of satellite depends only on a  (for same μ  )

In the above, μ = GM  where M  is the mass of the body at the focus of the ellipse and G  is the gravitational constant. h  is the specific mass angular momentum (moment of linear momentum) of the satellite. Hence the units of h2
 μ  is length.

To draw the locus of the satellite (the small body moving around the ellipse, all what we need is the eccentricity e  and a  , the major axes length. Then by changing the angle θ  the path of the satellite is drawn. I have a demo on this here

See http://nssdc.gsfc.nasa.gov/planetary/factsheet/earthfact.html for earth facts

This table below is from my class EMA 550 handouts (astrodynamics, spring 2014)

pict

3.3  Flight path angle for ellipse γ

pict

pict

To find γ  , if r  is given then use

       ∘ ----------  ∘ ----------
cosγ =    ---ap----=    a2(1−-e2)-
          r(2a− r)      r(2a − r)

If θ  is given, then use

tan γ = --esin-θ---
        1+ e cosθ

3.4  Parabolic trajectory

This diagram shows the parabolic trajectory

pict

3.5  Hyperbolic trajectory

This diagram shows the hyperbolic trajectory

pict

This diagram below from Orbital mechanics for Engineering students, second edition, by Howard D. Curtis, page 109

pict

3.6  showing that energy is constant

Showing that energy     v2   μ
ℰ = 2 −  r  is constant.

Most of such relations starts from the same place. The equation of motion of satellite under the assumption that its mass is much smaller than the mass of the large body (say earth) it is rotating around. Hence we can use ν = GM  and the equation of motion reduces to

¨   μ-
⃗r + r3⃗r = 0

In the above equation, the vector ⃗r  is the relative vector from the center of the earth to the center of the satellite. The reason the center of earth is used as the origin of the inertial frame of reference is due to the assumption that M ≫  m  where M  is the mass of earth (or the body at the focal of the ellipse) and m  is the mass of the satellite. Hence the median center of mass between the earth and the satellite is taken to be the center of earth. This is an approximation, but a very good approximation.

The first step is to dot product the above equation with ⃗r˙  giving

pict

And there is the main trick. We look ahead and see that ⃗r˙⋅ ¨⃗r = r˙r¨  but        (  2)
r˙r¨= ddt  ˙r2- and we also see that ˙ -μ     -˙r
⃗r ⋅r3⃗r = μr2   but   ˙r-   d-(−μ)
μ r2 = dt  r Hence equation 1 above can be written as

  (       )
d-  v2   r-
dt  2  − μ   = 0

Hence

      2
ℰ =  v-−  r-
     2    μ

Where ℰ is a constant, which is the total energy of the satellite.

3.7  Earth satellite Transfer orbits

Hohmann transfer

This diagram shows the Hohmann transfer

pict

3.7.1  Bi-Elliptic transfer orbit

pict

3.7.2  semi-tangential elliptical transfer

pict

3.8  Rocket engines, Hohmann transfer, plane change at equator

two cases: Hohmann transfer, 2 burns, or semi-tangential. All burns at equator.

pict

Rocket equation with plane change not at equator

pict

3.9  Spherical coordinates

pict

3.10  interplanetary transfer orbits

interplanetary hohmann transfer orbit, case one

pict

The following are the steps to accomplish the above. The first stage is getting into the Hohmann orbit from planet 1, then reaching the sphere of influence of the second planet. Then we either do a fly-by or do a parking orbit around the second planet. These steps below show how to reach the second planet and do a parking orbit around it.

The input is the following.

  1. μ1   planet one standard gravitational parameter
  2. μ2   planet two standard gravitational parameter
  3. μsun  standard gravitational parameter for the sun 1.327× 108   km
  4. r1   planet one radius
  5. r2   planet two radius
  6. alt1   original satellite altitude above planet one. For example, for LEO use 300 km
  7. alt2   satellite altitude above second planet. (since goal is to send satellite for circular orbit around second planet)
  8. R1   mean distance of center of first planet from the sun. For earth use AU  = 1.495978 × 108   km
  9. R
  2   mean distance of center of second planet from the sun. For Mars use 1.524  AU
  10. SOI
    1   sphere of influence for first planet. For earth use 9.24× 108   km
  11. SOI
    2   sphere of influence for second planet.

Given the above input, there are the steps to achieve the above maneuver

  1. Find the burn out distance of the satellite r  = r + alt
 bo    1     1
  2. Find satellite speed around planet earth (relative to planet)       ∘ ---
V   =   -μ1
 sat    rbo
  3. Find Hohmann ellipse a = R1+R2-
       2
  4. Find speed of satellite at perigee relative to sun            --------------
         ∘      ( 2   1)
Vperigee =   μsun  R1-− a
  5. Find speed of earth (first planet) relative to sun      ∘ μ---
V1 =   -sRu1n
  6. Find escape velocity from first planet V ∞,out = Vperigee − V1
  7. Find burn out speed at first planet by solving the energy equation V 2   μ    V2       μ
-b2o−  rb1o-= -∞,2out − SO1I1   for Vbo
  8. Find ΔV1   needed at planet one ΔV1  = Vbo − Vsat
  9. Find e  the eccentricity of the escape hyperbola     ∘ ------------
e =   1 + V2∞Vb2or2bo
             μ21
  10. Find the angle with the path of planet one velocity vector η = arccos(− 1)
             e
  11. Find the dusk-line angle θ = 1800 − η

The above completes the first stage, now the satellite is in the Hohmann transfer orbit. Assuming it reached the orbit of the second planet ahead of it as shown in the diagram above. Now we start the second stage to land the satellite on a parking orbit around the second planet at altitude alt2   above the surface of the second planet. These are the steps needed.

  1. Find the apogee speed of the satellite             --------------
          ∘      ( 2   1)
Vapogree =   μsun  R2-− a
  2. Find speed of second planet      ∘ μ---
V2 =   -sRun2-
  3. Find V ∞ entering the second planet sphere of influence V∞,in = V2 − Vapogree
  4. Find burn in radius where the satellite will be closest to the second planet. rbo = r1 + alt2
  5. Find burn out speed at second planet by solving the energy equation             2
V2bo-− μ2-= V∞,in−  -μ2-
2    rbo     2     SOI2   for Vbo
  6. Find impact parameter b  on entry to second planet SOI b = rboVbo
    V∞,in
  7. Find the required satellite speed around the second planet        ∘ ---
Vsat =   μ2-
         rbo
  8. Find ΔV
   2   needed at planet two ΔV   = V   − V
   2    sat   bo
  9. Find e  the eccentricity of the approaching hyperbola on second planet    ∘ ------------
         V∞2V-2bor2bo
e =  1 +    μ22
  10. Find the angle with the path of planet two velocity vector           (  1)
η = arccos − e
  11. Find the dusk-line angle for second planet            0
θ = 2η − 180

3.10.1  rendezvous orbits

Two satellite, walking rendezvous using Hohmann transfer

pict


Algorithm 1: Hohmann Walking Rendezvous Orbit, case 1
1:  function hohmann_walking_rendezvous(θ0,r,altitude,μ  )
2:   θ0 := θ01π80   ⊳  convert from degrees to radian
3:   ra := r + altitude
4:          ∘ r3a- T := 2π   μ  ⊳  period of circular orbit
5:   n :=  1
6:   done:=false
7:   while not(done) do
8:   TOF := (       )   N − 2θ0π  T
9:             (          ∘ --) a :=  solve  T OF  = 2π  a3  fora                         μ
10:   rp := 2a − ra
11:   if rp < r   then
12:   N  = N + 1
13:   else
14:   done:=true
15:   end if
16:   end while
17:            ∘ -- V     :=    μ  befor      h
18:            ∘ -(2---1)- Vafter :=    μ h − a
19:   ΔV  := 2(Vafter − Vbefore)
20:   return (TOF,  ΔV )
21:  end function

An example implementation is below

And calling the above

gives

 
{7123.89, -0.0467913}
Two satellite, separate orbits, rendezvous using Hohmann transfer, coplaner

pict


Algorithm 2: Hohmann rendezvous algorithm, case 1
1:  function hohmann_rendezvous_1(θ0,ra,rb,μ  )
2:   θ0 := θ01π80   ⊳  convert from degrees to radian
3:   a :=  ra+rb-       2   ⊳  Hohmann orbit semi-major axes
4:            ∘ a3- TOF  := π   μ  ⊳  time of flight on Hohmann orbit
5:          (              )              (ra+rb)3∕2 θH := π  1−    2rb ⊳  required phase angle before starting Hohmann transfer
6:         ∘ --- ω  :=    μ-  a      r3a  ⊳  angular speed of lower rad/sec
7:        ∘ --- ωb :=    μr3-         b  ⊳  angular speed of higher satellite rad/sec
8:   if θ0 ≤ θH   then ⊳  adjust initial angle if needed
9:   θ0 := θ0 + 2π
10:   end if
11:                θ0−θ waitxtime  :=  ωa−-Hωb  ⊳  how long to wait before starting Hohmann transfer
12:   waitxtime  := waitxtime + TOF  ⊳  now ready to go, add Hohmann transfer time
13:   return wait_time
14:  end function

An example implementation is below (in Maple)

And calling the above for two different cases gives (times in hrs)

And

Two satellite, separate orbits, rendezvous using bi-elliptic transfer, coplaner

pict

In this transfer, the lower (fast satellite) does not have to wait for phase lock as in the case with Hohmann transfer. The transfer can starts immediately. There is a free parameter N  that one select depending on fuel cost requiments or any limitiation on the first transfer orbit semi-major axes distance required. One can start with N =  0  and adjust as needed.


Algorithm 3: Hohmann rendezvous algorithm, case 2
1:  function hohmann_rendezvous_2(θ0,r1,r2,N,μ  )
2:   θ0 := θ01π80   ⊳  convert from degrees to radian
3:          (              )              (r1+r2)3∕2 θH := π  1−    2r2 ⊳  Find Hohmann ideal phase angle before transfer
4:   if θ = θ  andN  = 0  0   H   then ⊳  adjust for special case
5:        r+r a :=  22-1-
6:             (∘ --) TOF  := π    aμ3
7:   else
8:         ∘ μ-- ω2 :=    r32   ⊳  angular speed of slower satellite in rad/sec
9:   t := (2π−θ0)+2πN-  2        ω2   ⊳  find time of light of the slower satellite
10:         rt+r a1 :=  -2-1-
11:   a2 :=  rt+2r2-
12:             (∘ -------)              a31   a32 TOF  := π     μ + μ ⊳  time of flight for the fast satellite
13:   r := solve(t  = TOF  )forr  t          2            t  ⊳  Solve numerically for rt
14:   end if
15:   return TOF
16:  end function

An example implementation is below in Maple

And calling the above for two different cases gives

3.10.2  Semi-tangential transfers, elliptical, parabolic and hyperbolic

pict

pict

3.10.3  Lagrange points

pict

3.10.4  Orbit changing by low contiuous thrust

pict

3.10.5  References

  1. Orbital mechanics for Engineering students, second edition, by Howard D. Curtis
  2. Orbital Mechanics, Vladimir Chobotov, second edition, AIAA
  3. Fundamentals of Astrodynamics, Bates, Muller and White. Dover1971

4  Dynamic of flights

pict

4.1  Wing geometry

pict

Cr  below is the core chord of the wing.

pict

This is a diagram to use to generate equations of longitudinal equilibrium.

This distance is called the stick-fixed static margin km = (hn − h)¯c Must be positive for static stability

pict

4.2  Summary of main equations

This table contain some definitions and equations that can be useful.




#

equation

meaning/use




1

CL= ∂CL-α
     ∂α
  =CL  αα

  =a α

CL  is lift coefficient. α  is angle of attack. a  is slope ∂CL-
 ∂α  which is the same as C
  Lα




2

CLw = CLw αα

wing lift coefficient




3

CD =  CDmin + kC2L

drag coefficient




4

Cmw  = Cmacw + (CLw + CDminαw )(h − hnw) + (CL αw − CDw ) z¯c

pitching moment coefficient due to wing only about the C.G. of the airplane assuming small αw.  This is simplified more by assuming CDw αw ≪  CLw  and (CL αw − CDw ) ≪ 1




5

Cmw  = Cmacw + CLw (h − hnw)

simplified wing Pitching moment




6

Cmwb=Cmacwb  + CLwb (h− hnw )
               ∂CLwb
     =Cmacwb +  ∂αwb αwb (h − hnw)
     =C      + a  α   (h− h   )
        macwb   wb wb       nw

simplified pitching moment coefficient  due to wing and body about the C.G. of the airplane. αwb  is the angle of attack




7

        Lt
CLt = 12ρV2St

CLt  is the lift coefficient generated by tail. St  is the tail area. V  is airplane air speed




8

L = Lwb + Lt

total lift of airplane. Lwb  is lift due to body and wing and Lt  is lift due to tail




9

CL = CLwb +  SSt CLt

coefficient of total lift of airplane. CL
   wb  is coefficient of lift due to wing and body. CLt  is lift coefficient due to tail. S  is the total wing area. St  is tail area




10

                    1  2
Mt = − ltLt = − ltCLt 2ρV St

pitching moment due to tail about C.G. of airplane




11

      ---Mt--     ltSt
Cmt =  12ρV 2St¯c = − ¯c S CLt = − VH CLt

pitching moment coefficient due to tail.       ltSt
VH =  ¯c-S  is called tail volume




12

    l S
VH =-t¯cSt
¯   ¯ltSt
VH = ¯cS

introducing ¯
VH  bar tail volume which is VH  but uses ¯lt  instead of lt  . Important note. VH  depends on location of C.G., but ¯VH  does not. ¯lt = lt + (h − hn )¯c
               wb




13

Cmt = − ¯VH CLt + CLtSt(h − hnwb)
                    S

pitching moment coefficient due to tail expressed using V¯
 H  . This is the one to use.




14

Cm
   p

pitching moment coefficient due to propulsion about airplane C.G.




15

C  =  C    + C   + C
 m     mwb    mt     mp

total airplane pitching moment coefficient about airplane C.G.




16

Cm=Cmwb   + Cmt + Cmp
    [                       ]  [  ¯           St          ]
   = Cmacwb + CLwb (h − hnw) +  − VH CLt + CLtS (h − hnwb) + Cmp
             ◜------C◞L◟------◝
             (           St )
   =Cmacwb +   CLwb + CLt-S  (h − hnw) − ¯VH CLt + Cmp

   =Cmacwb + CL (h − hnw) − ¯VHCLt + Cmp

simplified total Pitching moment coefficient about airplane C.G.




17

∂Cm-= ∂Cmacwb-+ ∂CL-(h − h  )−  ¯V ∂CLt + ∂Cmp-
 ∂α     ∂α      ∂α       nw     H ∂α     ∂α
Cm α= ∂Cmacwb-+ CL α (h − hnw) − ¯VH ∂CLt+ ∂Cmp
        ∂α                        ∂α      ∂α

derivative of total pitching moment coefficient C
  m  w.r.t airplane angle of attack α




18

                (                        )
            -1-- ∂Cmacwb   ¯  ∂CLt  ∂Cmp-
hn = hnwb − ∂C∂αL    ∂α    − VH  ∂α  +   ∂α

location of airplane neutral point of airplane found by setting Cm α = 0  in the above equation




19

∂C∂mα-= ∂∂CαL(h − hn)

Cm α=CL α (h − hn)

rewrite of Cmα  in terms of hn  . Derived using the above two equations.




20

kn = hn − h

static margin. Must be Positive for static stability




Writing the equations in linear form

The following equations are derived from the above set of equation using what is called the linear form. The main point is to bring into the equations the expression for CLt  written in term of αwb.  This is done by expressing the tail angle of attack αt  in terms of αwb  via the downwash angle and the it  angle. ∂CLwb
 ∂αwb  in the above equations are replaced by awb  and ∂CLt
 ∂αt  is replaced by at  . This replacement says that it is a linear relation between CL  and the corresponding angle of attack. The main of this rewrite is to obtain an expression for Cm  in terms of αwb  where αt  is expressed in terms of αwb  , hence αt  do not show explicitly. The linear form of the equations is what from now on.




#

equation

meaning/use




1

CLwb = ∂CLwb-αwb
        ∂αwb
     = awbαwb
 C   = a α
  Lt    t t
C    = C    + ∂Cmp--α
  mp     m0p    ∂α

awb  is constant, represents ∂CLwb
 ∂αwb  and C
  m0p  is propulsion pitching moment coeff. at zero angle of attack α




2

α =  α  − i − ϵ
 t    wb   t
 ϵ = ϵ + ∂-ϵα
     0   ∂α  wb

main relation that associates αwb  with αt  . αwb  is the wing-body angle of attack, ϵ  is downwash angle at tail, and it  is tail angle with horizontal reference (see diagram)




3

CLt = atα[t   (       )         ]
                  ∂-ϵ
    = at αwb  1 − ∂α   − it − ϵ0

Lift due to tail expressed using αwb  and ϵ  (notice that αt  do not show explicitly)




4

       [                 ]
a = awb 1 + -at-St(1 − ∂ϵ)
            awbS      ∂α

a  defined for use with overall lift coefficient




5

    a  α
     w◜b◞◟w◝b  S
CL=  CLwb + -tS CLt
             St  [    (    ∂ϵ)        ]
  =awb αwb +  S at αwb 1−  ∂α − it − ϵ0
  =a α

  = (CL )     + aαwb
        αwb=0

overall airplane lift using linear relations




6

         atSt
α= αwb − a S (it + ϵ0)

pict

overall angle of attack α  as function of the wing and body angle of attack αwb  and tail angles




7

Cm=Cm0   + ∂C∂mα-α = Cm0 + Cm αα

Cm= C¯m0  + ∂C∂mα-αwb = ¯Cm0 + Cm ααwb

overall airplane pitch moment. Two versions one uses αwb  and one uses α




8

                        (      )
Cm α=a (h− hnwb )− at¯VH  1−  ∂ϵ +  ∂Cmp--
                          (  ∂α  )  ∂α
Cm α=awb (h− hnwb )− atVH  1−  ∂∂αϵ + ∂Cm∂pα-

Two versions of ∂∂Cαm-  one for αwb  and one one uses α




9

                       ¯          [    atSt (   -∂ϵ)]
Cm0=Cmacwb  + Cmop + atVH (ϵ0 + it) 1− a S  1− ∂α
¯Cm  =Cm     + C¯m   + atVH (ϵ0 + it)
   0    acwb     op

Cm
   0   is total pitching moment coef. at zero lift (does not depend on C.G. location) but C¯m0   is total pitching moment coef. at αwb = 0  (not at zero lift). This depends on location of C.G.




10

¯Cm0p = Cm0p + (α − αwb)∂Cm∂pα-




11

hn=hn   +  atV¯ (1 − ∂ϵ) − 1∂Cmp-
      wb   a H      ∂α    a  ∂α (      )
  =hnwb +  --[---aatS-(--∂ϵ)] ¯VH 1 − ∂∂ϵα  − ---[--a-1S(---∂ϵ)]∂C∂mαp-
           awb 1+awtb-tS 1−∂α                awb1+ awtb-St1− ∂α

Used to determine  hn




4.2.1  definitions

  1. Remember that for symmetric airfoil, when the chord is parallel to velocity vector, then the angle of attack is zero, and also the left coefficient is zero. But this is only for symmetric airfoil. For the common campbell airfoil shape, when the chord is parallel to the velocity vector, which means the angle of attack is zero, there will still be lift (small lift, but it is there). What this means, is that the chord line has to tilt down more to get zero lift. This extra tilting down makes the angle of attack negative. If we now draw a line from the right edge of the airfoil parallel to the velocity vector, this line is called the zero lift line (ZLL) see diagram below.
    Just remember, that angle of attack (which is always the angle between the chord and the velocity vector, the book below calls it the geometrical angle of attack) is negative for zero lift. This is when the airfoil is not symmetric. For symmetric airfoil, ZLL and the chord line are the same. This angle is small,    0
− 3   or so. Depending on shape. See Foundations of Aerodynamics, 5th ed, by Chow and Kuethe, here is the diagram.

    pict

  2. stall from http://en.wikipedia.org/wiki/Stall_(flight)
    In fluid dynamics, a stall is a reduction in the lift coefficient
    generated by a foil as angle of attack increases.[1] This occurs when
    the critical angle of attack of the foil is exceeded. The critical
    angle of attack is typically about 15 degrees, but it
    may vary significantly depending on the fluid, foil, and Reynolds number.
  3. Aerodynamics in road vehicle wiki page
  4. some demos relating to airplane control http://demonstrations.wolfram.com/ControllingAirplaneFlight/

    http://demonstrations.wolfram.com/ThePhysicsOfFlight/

  5. http://www.americanflyers.net/aviationlibrary/pilots_handbook/chapter_3.htm
  6. Lectures Helicopter Aerodynamics and Dynamics by Prof. C. Venkatesan, Department of Aerospace Engineering, IIT Kanpur http://www.youtube.com/watch?v=DKWj2WzYXtQ&list=PLAE677E56C97A7C7D
  7. http://avstop.com/ac/apgeneral/terminology.html has easy to understand definitions airplane geometry. ”The MAC is the mean average chord of the wing”
  8. http://www.tdmsoftware.com/afd/afd.html airfoil design software

4.3  images and plots collected

These are diagrams and images collected from different places. References is given next to each image.

pict

pict

pict

pict

This below from http://www.grc.nasa.gov/WWW/k-12/UEET/StudentSite/dynamicsofflight.html

pict

http://www.grc.nasa.gov/WWW/k-12/airplane/alr.html

pict

From http://en.wikipedia.org/wiki/Lift_coefficient and http://en.wikipedia.org/wiki/File:Aeroforces.svg pict

from http://adg.stanford.edu/aa241/drag/sweepncdc.html

pict

Images from http://adamone.rchomepage.com/cg_calc.htm and Flight dynamics principles by Cook, 1997.

pict

From http://chrusion.com/BJ7/SuperCalc7.html

pict

From http://www.willingtons.com/aircraft_center_of_gravity_calcu.html

pict

From http://www.solar-city.net/2010/06/airplane-control-surfaces.html nice diagram that shows clearly how the elevator causes the pitching motion (nose up/down). From same page, it says ”The purpose of the flaps is to generate more lift at slower airspeed, which enables the airplane to fly at a greatly reduced speed with a lower risk of stalling.”

pict

Images from flight dynamics principles, by Cook, 1997.

pict

Images from Performance, stability, dynamics and control of Airplanes. By Pamadi, AIAA press. Page 169. and http://www.americanflyers.net/aviationlibrary/pilots_handbook/chapter_3.htm

pict

Image from http://www.americanflyers.net/aviationlibrary/pilots_handbook/chapter_3.htm

pict

Image from http://www.americanflyers.net/aviationlibrary/pilots_handbook/chapter_3.htm

pict

Image from FAA pilot handbook and http://www.youtube.com/watch?v=8uT55aei1NI

pict

Image http://www.youtube.com/watch?v=8uT55aei1NI and http://www.youtube.com/user/DAMSQAZ?feature=watch

pict

pict

4.4  Some strange shaped airplanes

Image http://edition.cnn.com/2014/01/16/travel/inside-airbus-beluga/index.html?hpt=ibu_c2

pict

pict

Image from http://edition.cnn.com/2014/01/16/travel/inside-airbus-beluga/index.html?hpt=ibu_c2

pict

Image from http://www.nasa.gov/centers/dryden/Features/super_guppy.html

pict

Image from http://www.aerospaceweb.org/question/aerodynamics/q0130.shtml ”Boeing Pelican ground effect vehicle”

pict

4.5  links

  1. https://3dwarehouse.sketchup.com/search.html?redirect=1&tags=airplane

4.6  references

  1. Etkin and Reid, Dynamics of flight, 3rd edition.
  2. Cook, Flight Dynamics principles, third edition.
  3. Lecture notes, EMA 523 flight dynamics and control, University of Wisconsin, Madison by Professor Riccardo Bonazza
  4. Kuethe and Chow, Foundations of Aerodynamics, 4th edition