Trajectory and Orbits - Computational Physics - Lecture Notes | PHYS 4007, Study notes of Physics

Download Trajectory and Orbits - Computational Physics - Lecture Notes | PHYS 4007 and more Study notes Physics in PDF only on Docsity!

PHYS-4007/5007: Computational Physics

Course Lecture Notes

Section X

Dr. Donald G. Luttermoser East Tennessee State University

Version 4.

Abstract

These class notes are designed for use of the instructor and students of the course PHYS-4007/5007: Computational Physics taught by Dr. Donald Luttermoser at East Tennessee State University.

X–2 PHYS-4007/5007: Computational Physics

Figure X–1: The orthogonal (Cartesian) coordinate system.

y

z

x

P y^

z^

x^

1. In classical mechanics, the static potential field is related to a conservative force by the equation F^ ~ = −∇~U , (X-2) where U is the potential energy of the field (a scalar), and the “del” operator acting on a scalar is referred to as taking the gradient of the scalar =⇒ it converts a scalar to a vector. The “del” operator takes first derivatives on each coordinate of the vector space. As such, it has a different form depending on the coordinate system: a) Orthogonal (Cartesian) Coordinates (x, y, z):

∇^ ~ = ˆx ∂ ∂x

• ˆy ∂ ∂y
• ˆz ∂ ∂z

, (X-3)

such that ∇^ ~f(x, y, z) = ˆx ∂f ∂x

• ˆy ∂f ∂y
• ˆz ∂f ∂z

b) Spherical-Polar (Spherical) Coordinates (r, θ, φ):

∇^ ~ = ˆr ∂ ∂r

• θˆ

r

∂θ

• φˆ

r sin θ

∂φ

. (X-4)

Donald G. Luttermoser, ETSU X–

Figure X–2: The spherical-polar coordinate system.

y

z

x

P

r

^r

^ θ

^ φ

θ

φ

Figure X–3: The circular-cylindrical coordinate system.

y

z

x

P

ρ

ρ

z

z

^ ρ

z^ ^ φ

φ

Donald G. Luttermoser, ETSU X–

iii) Integrating Eq. (X-8), we get

− G m 1 m 2 r

2 1

∫ (^2) 1 dU^ =^ U^2 −^ U^1 , or U 1 − U 2 = G m 1 m 2

r 2

r 1

)

. (X-9)

b) From Eq. (X-9), we can immediately see that the potential energy for a gravitational field takes on the form

U = − G m 1 m 2 r

. (X-10)

i) From this equation we see that U → 0 as r → ∞.

ii) Also we can see that gravitational potential en- ergy is a negative energy.

1. Besides potential energy, any body in motion has a kinetic energy associated with it. a) The left-hand side of Eq. (X-7) is the definition of the work done on a particle by a force field:

W ≡

∫ (^2) 1 F^ ~ · d~r = U 1 − U 2. (X-11)

b) Using one of the forms of Newton’s Second Law of Motion in Eq. (X-1) and rewriting d~r as (d~r/dt)dt, we can write

F^ ~ · d~r =

( m d~v dt

) ·

( (^) d~r dt dt

) = m d~v dt · ~v dt.

Note that d dt (~v · ~v) = ~v · d~v dt

• d~v dt · ~v = 2d~v dt · ~v ,

so d~v dt · ~v =

d dt (~v · ~v)

X–6 PHYS-4007/5007: Computational Physics

and hence,

F^ ~ · d~r = 1 2 m d dt (~v · ~v) dt =

m d dt

( v^2

) dt

= d

2 mv

2

)

. (X-12)

c) Now using Eq. (X-12) in the first part of Eq. (X-11), we get W ≡

∫ (^2) 1 F^ ~ · d~r =

∫ (^2) 1 d

mv^2

)

mv 22 −

mv^21 = T 2 − T 1 , (X-13) where ‘T ’ represents kinetic energy.

d) Inserting this value for the work in Eq. (X-11), we see W = T 2 − T 1 = U 1 − U 2 or T 1 + U 1 = T 2 + U 2 (X-14) E 1 = E 2 , where E represents the total mechanical energy of the system.

e) Eq. (X-14) is called the conservation of mechanical energy, which is valid only for a conservative force → forces that do not depend on time nor does the work done by such a force depend upon the path taken in a trajec- tory.

1. Whereas the force is equal to the negative gradient of the poten- tial energy (see Eq. X-2), the acceleration due to a conservative force can be determined from the potential Φ of the force (‘po- tential’ is the potential energy per unit mass). For gravity ~g ≡ −∇~Φ , (X-15)

X–8 PHYS-4007/5007: Computational Physics

1. When dealing with motion in gravitational fields, there are two regimes that are typically encountered: (1) trajectories (near a surface of a large mass, e.g., Earth) and (2) orbits (where 2 masses can be considered as point-like). a) For orbits, we use the general form of the gravitational potential as described in Eqs. (X-10,16):

Ug = − GMm r

, (X-21)

where we are now using M (the larger mass) for m 1 and m (the smaller mass) for m 2 and

Φ = −

GM

r

. (X-22)

b) For trajectories, typically the maximum height (ymax = h) reached is small with respect to R⊕ and hence g ≈ constant. As such, we can write Eqs. (X-9 and X-11) as

W = U 1 − U 2 = GM⊕m

r 2

r 1

) .

i) If we take point ‘2’ to be the Earth’s surface and ‘1’ to be the position of the projectile, then

∆U = G M⊕ m

R⊕

R⊕ + h

)

= G M⊕ m

  R⊕^ +^ h R⊕ (R⊕ + h) −^

R⊕

R⊕ (R⊕ + h)

 

= G M⊕ m

  R⊕^ +^ h^ −^ R⊕ R⊕ (R⊕ + h)

 

= G M⊕ m

  h R⊕ (R⊕ + h)

  (^).

ii) If R⊕ is much greater than h (which it will be for experiments near the Earth’s surface), h  R⊕. As

Donald G. Luttermoser, ETSU X–

such, R⊕+h ≈ R⊕ and the equation above becomes

∆U = G M⊕ m

  h R⊕ (R⊕)

  (^) = G M⊕ m

  h R^2 ⊕

 

G M⊕ m h R^2 ⊕ = m

G M⊕

R^2 ⊕

h. (X-23)

iii) Using Eqs. (X-15 and X-16), we see that

~g = −∇~Φ = −∇~

( G M

⊕ r

)

= − d dr

(G M

⊕ r

) ˆr = G M⊕ r^2 rˆ (X-24)

and at the Earth’s surface,

g = G M⊕ R^2 ⊕

, (X-25)

where ‘g’ is referred to as the Earth’s surface grav- ity.

iv) Using Eq. (X-25) in Eq. (X-23), we finally get

∆U = mgh = mg∆y ,

where ∆y = y − y◦ = h is just the change in height from our initial position y◦ (typically the ground) and y is an arbitrary position in the trajectory above y◦.

v) If we arbitrarily set y◦ = 0, then the potential at that position is zero, and y represents the posi- tion above the ground (y◦). As such, the potential energy becomes

U = mgy. (X-26)

Donald G. Luttermoser, ETSU X–

ter of rotation (perpendicular to the axis of rota- tion).

iii) F~cor = − 2 m~ω × ~vr is the Coriolis force, which results from the Earth’s rotation (~vr ≡ Earth’s ro- tational velocity) =⇒ this “force” arises when an attempt is made to describe motion relative to a ro- tating body (i.e., the ground moves out from under you when the projectile is in the air).

iv) Note that the centrifugal and Coriolis forces are not forces in the usual sense of the word. They are only introduced so that the inertial (non-accelerat- ing) frame equation F^ ~ = m~af (i.e., Newton’s 2nd law) can have a non-inertial (accelerating) frame analogous equation: F^ ~eff = m~ar , so that F^ ~eff = m~af + (non-inertial terms).

B. Numerical Solutions for Trajectories.

1. Trajectories with h  R⊕ and x  R⊕. a) Combining Eqs. (X-1, 27, & 28), Newton’s 2nd law be- comes: ∑ (^) ~ F = F~g + F~r = m d^2 ~r dt^2 or m d^2 ~r dt^2 =^ −mg^ yˆ^ −^ mkv

n ~v v ,^ (X-30) where we have defined +ˆy in the upward direction.

X–12 PHYS-4007/5007: Computational Physics

b) Since the atmospheric drag force is a function of ~v, it is more convenient to write Eq. (X-30) as

m d~v dt = −mg yˆ − mkvn^ ~v v d~v dt =^ −g^ yˆ^ −^ kv

n ~v v (X-31) i) At low velocities, n ≈ 1 and the magnitude of the drag force follows Fr ≈ −B 1 v (X-32) =⇒ this is known as Stoke’s law.

ii) As v increases, n → 2 and the drag force follows Fr ≈ −B 2 v^2. (X-33)

iii) As such, we can write a general form of the drag force as Fr ≈ −B 1 v − B 2 v^2 (X-34) or Fr = −

∑N i=

Bi vi^ (X-35)

=⇒ a simple power series in v (note that Bi → 0 faster than vi^ → ∞ as i → ∞).

iv) The Bi’s are related to the drag coefficient k in Eqs. (X-30 & 31).

c) Since we are dealing with projectiles here, the drag force will simply be described by Eq. (X-33). As such, we need to solve each component of Eq. (X-31) — hence, we need to break the drag force into its component forces: Fr = −mkv^2.

X–14 PHYS-4007/5007: Computational Physics

g) Eqs. (X-43) can be then solved as an initial value problem as described in §IX with x◦, y◦, vx,◦, and vy,◦ supplied by the user. i) The ∆t steps are chosen to give the errors that follow Eq. (VII-15).

ii) Calculations are carried out until a certain τ = total time is reached or some condition of xi, yi, vx,i, or vy,i is satisfied.

h) But what is the drag coefficient k? i) Since the projectile is trying to push air of mass dmair out of the way, where dmair ≈ ρ A v dt , (X-44) where ρ is the density of the air and A is the frontal area, we can guess that k is a function of ρ [and possibly t if the object is rotating → then A = A(t)].

ii) k(y = 0) = k◦ is usually given (based on air tun- nel measurements), and the following expression is used for k(y):

k(y) = ρ(y) ρ◦ k◦ , (X-45)

where ρ◦ is the density of air when the k◦ measure- ment was made.

i) Now we need a description of ρ(y)! i) One could supply a data table of ρ as a func- tion of height (see Appendix B.1 in Fundamentals of Atmospheric Modeling by Mark Jacobson, 1999, Cambridge University Press).

Donald G. Luttermoser, ETSU X–

ii) One could solve the following set of differential equations to determine ρ(y):

P = NkB T (ideal gas law) (X-46) ρ =

∑ Ni mi (mass conservation) (X-47) dP dy = −ρ g (hydrostatic equilibrium) (X-48) dT T

  R cp

  dP P (Poisson’s equation) (X-49) dQ = cp dT − α dP (energy conserv.) (X-50)

where P = pressure, T = temperature, N = total particle density, ρ total mass density, Ni = num- ber density of species i, mi = mass of species i, kB = Boltzmann’s constant, g = surface gravity, R = universal gas constant, cp = specific heat of air at constant pressure, Q = total heat (determined from solar radiation incident on the Earth’s atmo- sphere), and α = specific volume of air.

iii) Note that the solution to these equation would still only be an approximation since we have left out condensation, evaporation, sublimation, chemical reactions, and wind from these equations.

1. Trajectories with h ∼ R⊕ and x ∼ R⊕. a) Rotating Coordinate Systems. i) Let’s consider 2 sets of coordinate axes: =⇒ one ‘fixed’ = inertial frame (the primed [′] coordinates), =⇒ one ‘rotating’ with respect to the fixed system and possibly in linear motion with respect to the fixed frame = noninertial frame.

Donald G. Luttermoser, ETSU X–

iv) The time rate of change of ~r as measured in the fixed frame is ( (^) d~r dt

) fixed

d~θ dt × ~r = ~ω × ~r , (X-53)

since the angular velocity ~ω is defined by

~ω ≡ d~θ dt

. (X-54)

v) If point P has a velocity (d~r/dt)rot with respect to the rotating system, this velocity must be added to ~ω × ~r to obtain the time rate of change of ~r in the fixed system: ( (^) d~r dt

) fixed

( (^) d~r dt

) rot

• ~ω × ~r. (X-55)

vi) This expression is not just limited to the displace- ment vector ~r, in fact, for any arbitrary vector Q~, we have   d Q~ dt

  fixed

  d Q~ dt

  rot

• ~ω × Q~. (X-56)

vii) Note that the angular acceleration ~ω˙ is the same in both the fixed and rotating systems: ( (^) d~ω dt

) fixed

( (^) d~ω dt

) rot

• ~ω × ~ω ( (^) d~ω dt

) fixed

( (^) d~ω dt

) rot

≡ ~ω .˙ (X-57)

viii) As such, the velocity of point P as measured in the fixed coordinate system is   d~r^

′ dt

  fixed

  d R~ dt

  fixed

(d~r dt

) fixed

X–18 PHYS-4007/5007: Computational Physics

  d~r^

′ dt

  fixed

  d R~ dt

  fixed

(d~r dt

) rot

• ~ω × ~r.

If we define ~vf ≡ ~r˙f ≡

  d~r^

′ dt

  fixed

(X-58)

V^ ~ ≡ ~R˙f ≡

  d R~ dt

  fixed

(X-59)

~vr ≡ ~r˙r ≡

( (^) d~r dt

) rot

(X-60)

we may write

~vf = V~ + ~vr + ~ω × ~r , (X-61)

where ~vf = velocity relative to the fixed axes V^ ~ = linear velocity of the moving origin ~vr = velocity relative to the rotating axes ~ω = angular velocity of the rotating axes ~ω × ~r = velocity due to the rotation of the moving axes.

b) The Coriolis Force. i) Newton’s 2nd law F~ = m~a is only valid in an inertial frame, therefore

F^ ~ = m~af = m

( (^) d~v f dt

) fixed

, (X-62)

where the differentiation must be carried out with respect to the fixed system.

ii) If we limit ourselves to cases of constant angular acceleration ( ˙ω = 0), using Eq. (X-61) we can write

F^ ~ = mR~¨f + m

(d~v r dt

) fixed

• m~ω ×

( (^) d~r dt

) fixed

. (X-63)