Appendix F v6

This has been installed, to not edit here! Appendix F v4 PhL 12.18.16 Appendix F : The Dumbbell (Tethered) satellite as an example of rotating frame analysis 1 F.1 Kinematics of the satellite in rotating Frame S 2 F.2 Angular momentum of the satellite and its time derivative in Frame S 6 F.3 The torque on the Dumbbell Satellite in Frame S' 7 F.4 The fictitious torque on the satellite in Frame S 11 F.5 Equations of Motion for the satellite in Frame S (Spherical Coordinates) 15 F.6 Force analysis of the satellite in Frame S (Spherical Coordinates) 18 F.7 Numerical solutions of the equations of motion (Spherical Coordinates) 25 F.8 Force analysis of the satellite in Frame S (Cartesian Coordinates) 29 F.9 Verification of the Cartesian equations of motion and stick tension 34 F.10 Numerical solutions of the equations of motion (Cartesian Coordinates) 38 Appendix F : The Dumbbell (Tethered) satellite as an example of rotating frame analysis This is a very long appendix so we provide an overview: Section F.1 lays out kinematic details of the coordinates we use and discusses simple geometric facts of the satellite. We use the "swap notation" wherein Frame S' is the inertial frame and Frame S is the rotating frame. The masses m1 and m2 can be equal or different. Section F.2 computes the angular momentum L(0) of the satellite and its time derivative (0) in rotating Frame S where the origin of Frame S is used as a reference point. Section F.3 computes the true torque N'(b) on the satellite due to the Earth's gravitational attraction acting on the two masses. The resulting torque is computed in inertial Frame S' and has only a component as shown in (F.3.7). This computation is done in two extra ways to verify the result. In the far approximation, the expression for N'(b) simplifies to that shown in (F.3.13). Section F.4 then computes the fictitious torque which appears in rotating Frame S. This torque N(0)fict is stated in (F.4.11) and has and components. An interpretation is provided. Section F.5 then uses "Newton's Angular Law" (0) = N'(b) + N(0)fict to obtain the angular equations of motion for the satellite as stated in (F.5.7) using the far approximation. Certain special-case solutions are extracted which demonstrate the notion of in-plane and out-of-plane small-θ libration of the satellite with frequencies given in (F.5.10) and (F.5.13). Section F.6 basically starts over working only with forces, not torques. The two angular equations of motion obtained in Section F.5 are again obtained, and a third equation determines the tension T in the dumbbell stick. One can interpret this tension T as a tidal force acting on mass m1 which is distance r1 from the Frame S center of mass. An equal and opposite tidal force acts on mass m2 at the other end of the stick. When the satellite is vertically aligned, one finds T = 3m1 ω2r1 in (F.6.30). The tension T does not appear in the angular analysis of Section F.5 since T makes no contribution to torque about the Frame S origin, being aligned with the stick. Tension T is part of the radial force equation (F.6.19) of Section F.6, but there is no radial equation in Section F.5. There the stick is regarded just as a constraint, and it is typical for the forces of constraint not to be determined in the simplest analysis. Section F.7 provides some Maple numerical solutions of the far-approximation θ,φ equations of motion of the satellite stated in (F.5.7). The two libration modes are verified, and a more general case is examined. Section F.8 reworks the satellite force analysis entirely in Cartesian coordinates. The equations of motion for the variables x,y,z are shown in (F.8.15) with tension T provided by (F.8.20). Section F.9 shows that the x,y,z equations of motion of Section F.8 are entirely equivalent to the θ,φ equations of motion of Section F.5. Section F.10 examines various numerical solutions for the dumbbell satellite, working now in Cartesian coordinates. The libration modes are again examined, and a more general solution is studied. F.1 Kinematics of the satellite in rotating Frame S This section uses the "swap notation" described in the Summary for Section 1 where the rotating frame is Frame S and the non-rotating frame is Frame S'. We do this to reduce the number of primes since most activity will be in the rotating satellite Frame S. We now place the dumbbell satellite in a more general orientation than it was in Appendix D : (F.1.1) Description of the Figure Half the battle is having a clear picture of what is going on and we shall expend many words to describe the above drawing. It shows the dumbbell satellite in orbit in a completely arbitrary orientation. The two masses m1 and m2 are connected by a massless stick (not shown) of length r1+r2 = s . If we find that this stick is always in tension for some situation, we can replace it in that situation with a non-stretching massless tether. If the stick were to go into compression, the replacement tether would lose its linear shape and we don't want to deal with that problem. As shown much later in Section F.10, the stick is always in tension for normal situations. The gray-filled triangle is a part of the plane φ = constant which has normal vector . This plane is not in the plane of paper. The fill region contains two non-right triangles shown in blue and red. The red triangle is on the viewer's side of paper, while the blue one lies behind the paper. Each of these triangles contains the vector b as an edge. Frame S' is an inertial frame whose center is the center of the Earth and which is assumed fixed with respect to the stars. Frame S is a non-inertial frame whose center is located at the satellite center of mass point. We make the assumption discussed above (D.2.2) that we can ignore the tiny offset between the center of mass and center of gravity of the satellite, so we then regard the Frame S origin as travelling in a circular orbit around the Earth. Recall from (D.3.16) that this offset is about 3 microns for s = 10 meters. At time t = 0 shown in the figure, the axes of both Frames align with each other. At any time, axes y,z and y',z' are in the plane of paper. The unit vector always points from the center of the Earth to the origin of Frame S. so b = b at any time. The axes x,x' always point directly out of the plane of paper so = ' . The orbital rotation vector ω = ω also points out of the plane of paper and recall that ω = 2π/T where T is about 88 minutes for a low-Earth orbit. We have in mind an orbit at any altitude, but we do assume a circular orbit. The circled dot at the bottom is the center of the Earth (mass ME) and we have drawn the Earth surface in green. The blue circle is the orbit of the center of mass of the satellite (the Frame S origin) and it lies in the plane of paper. Note that ω is for the orbit and has nothing whatsoever to do with the rotation of the Earth. The above picture is the same whether or not the Earth rotates at its 24 hour ωE rate about some obscure E axis (not shown). One could try to treat the satellite as a reduced-mass single-particle system as is done for planetary orbits, but even if possible this would obscure details we want to be visible. We really have here a 3-body problem where the three bodies are point masses, but there is a constraint (the stick), so perhaps it is a 2 1/2-body problem. If the dumbbell were treated as a rigid object, it is then a 2-body problem where one body is not a point mass. Naming of coordinates In Frame S mass m1 has spherical coordinates (r1,θ1,φ1) while mass m2 has coordinates (r2,θ2,φ2). The reader is now forewarned about our upcoming slipshod notation. We define (θ,φ) ≡ (θ1,φ1). Thus Fig (F.1.1) shows θ,φ and not θ1,φ1. The three unit vectors of spherical coordinates for mass m1 will always be called , , and never 1, 1, 1. These spherical coordinates are in the "physics" convention: angle θ is the polar angle down from the "vertical" z axis, φ is the azimuthal angle measured from the x axis toward the y axis. See Appendix E.1 regarding conventions. Because r1 and r2 are collinear, we know that θ2 = π -θ and φ2 = π + φ. Since the Frame S origin is the center of mass, we also know that r2 = (m1/m2)r1. Our strategy is to avoid the subscript-2 coordinates whenever possible and express everything in terms of the mass m1 coordinates (r1,θ,φ). The mass m2 unit vectors are related to the m1 unit vectors by 2 = - , 2 = + , and 2 = - . Each unit vector points toward increasing parameter value for its associated mass. To see these last relations, it helps to stare at Fig (F.1.1) and think of the gray triangle as being in the plane of paper. We shall only make use of 2 = - below. See Appendix E regarding the unit vectors , and . To summarize : r1 = [x1,y1,z1] = (r1,θ,φ) = coordinates of mass m1 velocity = v1 r2 = [x2,y2,z2] = (r2,θ2,φ2) = coordinates of mass m2 velocity = v2 θ2 = π-θ φ2 = φ+π r2 = (m1/m2) r1 2 = - 1 = - . (F.1.2) In Frame S mass m1 is constrained to lie on a sphere of radius r1, while mass m2 is constrained to lie on a sphere of radius r2. The picture assumes m1 < m2 so r1 > r2. If m1 lies at a point on its sphere, m2 lies on the inverse point but on its sphere, as the spherical coordinates above show. If the masses are the same, the two spheres coincide. Due to these constraints on the vectors r1 and r2, the corresponding velocity vectors v1 and v2 of the two masses must be tangential to their respective spheres. Thus, for example, for mass m1 we can write v1 = v1θ + v1φ , whereas r1 = r1 . (F.1.3) Some Basic Kinematic Facts The internal angles of the red triangle are β1, α1 and π-θ at the origin. Thus we know from a Law of Cosines that r'12 = b2 + r12 - 2br1cos(π-θ) = b2 + r12 + 2br1cosθ . The internal angles of the blue triangle are β2, α2 and θ at the origin. Thus, r'22 = b2 + r22 - 2br2cosθ . The three Laws of Sines for the two triangles tells us that = = = // red triangle = = . // blue triangle Looking at the drawing it is clear that b + r1 = r'1 and b + r2 = r2', so one can write r'2 - r2 = r'1 - r1 = b . Recall the center of mass condition from (D.2.8) that m1r1 = - m2r2 m1r1 = m2r2 , r2/r1 = m1/m2 , 2 = – . (D.2.8) Applying the Frame S time derivative ∂S gives similar results for r → v and r → a, see below (velocity and acceleration). Summary: r'12 = b2 + r12 + 2br1cosθ (F.1.4) r'22 = b2 + r22 - 2br2cosθ (F.1.5) = = // red triangle (F.1.6) = = // blue triangle (F.1.7) r'2 - r2 = r'1 - r1 = b r'1 = b + r1 r'2 = b + r2 (F.1.8) m1r1 = - m2r2 m1r1 = m2r2 , r2/r1 = m1/m2 , 2 = – 1 (F.1.9) m1v1 = - m2v2 m1v1 = m2v2 , v2/v1 = m1/m2 , 2 = – 1 (F.1.10) m1a1 = - m2a2 m1a1 = m2a2 , a2/a1 = m1/m2 , 2 = – 1 (F.1.11) F.2 Angular momentum of the satellite and its time derivative in Frame S Our upcoming path is to obtain the equations of motion for the satellite in terms of θ, φ coordinates. As a demonstration of the method of Section 11.3 we shall do this first using Newton's rotational second law with fictitious torques. Later in Section F.6 we will do this again using Newton's linear second law with fictitious forces as a demonstration of the method of Section 8.1. For angular momentum (and later torque) we take our reference point to be c = 0 in Frame S (the origin) and thus c' = b in Frame S'. Then (we show all detail in this first calculation), L(0) = r1 x p1 + r2 x p2 = m1 r1 x v1 + m2 r2 x v2 // dim = ML2/T = m1r1 x v1 + m2[-(m1/m2)r1] x [-(m1/m2)v1] // (F.1.9) and (F.110) = m1r1 x v1 + [m1r1] x [(m1/m2)v1] = m1{ r1 x v1 + (m1/m2) r1 x v1 } = m1[ 1 + (m1/m2) ] r1 x v1 = m1[ (m2+m1)/m2 ] r1 x v1 = (m1/m2) M r1 x v1 = (m1/m2) Mr1 x v1 // M ≡ m1+m2 (F.2.1) = (m1/m2) Mr1 x (v1θ + v1φ) = (m1/m2) Mr1 [ v1θ x + v1φ x ] // (F.1.3) = (m1/m2) M r1 (v1θ – v1φ ) // (E.2.12) = (m1/m2) M r12 ( – sinθ ) . // (E.3.5) dim=ML2/T (F.2.2) The dumbbell has no angular momentum around the axis which seems very reasonable since it consists of two point masses aligned with . Taking a ∂S time derivative gives the rate of change of angular momentum in Frame S (all variables here are the "natural" ones in Frame S ), (0) = ∂S{ (m1/m2) M r1 x v1 } // (F.2.1) = (m1/m2) M (v1 x v1 + r1 x a1 ) = (m1/m2) Mr1 ( x a1 ) = (m1/m2) Mr1 x [ ar + aθ + aφ ] // (E.3.6) = (m1/m2) Mr1 [ aθ - aφ ] // (E.2.12) and then (E.3.6) for next line = (m1/m2) Mr12 [ (- 2 sinθ cosθ) – (2 cosθ + sinθ) ] . (F.2.3) Here are the conclusions so far: L(0) = (m1/m2) M r12 ( – sinθ ) (F.2.2) (0) = (m1/m2) Mr12 [ (- 2 sinθ cosθ) – (2 cosθ + sinθ) ] . (F.2.3) One obvious statement can be made looking at these equations: there is no angular momentum about the axis and this vanishing angular momentum never changes. The equation of motion for the satellite within Frame S is given by (11.3.4), but converted to swap notation, N(0)eff = (0) . (11.3.4)s (F.2.4) Our next task then is to compute the total effective torque on the satellite in Frame S which from (11.3.5) is (again converted to swap notation), N(0)eff = N'(b) + N(0)fict . (11.3.5)s (F.2.5) Here N'(b) is the Frame S' torque on the satellite relative to point b in Frame S', and N(0)fict is the fictitious torque that arises because Frame S is a rotating frame of reference. F.3 The torque on the Dumbbell Satellite in Frame S' The torque (in Frame S') of the Earth on the satellite (relative to the Frame S origin) is given by N'(b) = r1 x F1 + r2 x F2 // dim = L3 M /T (F.3.1) where all four of these vectors are shown in Fig (F.1.1). Recall from (D.2.13) that F1 = - (GMEm1/r'12) '1 = - (GMEm1/r'13) r'1 F1 ≡ |F1| = (GMEm1/r'12) F2 = - (GMEm2/r'22) '2 = - (GMEm2/r'23) r'2 F2 ≡ |F2| = (GMEm2/r'22) . (F.3.2) // The stick tension T exerts no torque on the satellite masses since T is collinear with r and r2. Therefore, r1 x F1 = r1 x [ - (GMEm1/r'13) r'1 ] = - (GMEm1/r'13) r1 x r'1 r2 x F2 = r2 x [ - (GMEm2/r'23) r'2 ] = - (GMEm2/r'23) r2 x r'2 . (F.3.3) Using (F.1.8) we evaluate the cross products making use of the first line of (E.2.15) for x , r1 x r'1 = r1 x (b + r1) = r1 x b = r1b x = r1b (-sinθ ) = -r1bsinθ r2 x r'2 = r2 x (b + r2) = r2 x b = [-r2] x [b] = -r2b x = -r2b(-sinθ ) = r2bsinθ . (F.3.4) Then inserting (F.3.4) into (F.3.3), r1 x F1 = (GMEm1/r'13) r1b sinθ (F.3.5) r2 x F2 = - (GMEm2/r'23) r2bsinθ = - (GMEm1/r'23) r1bsinθ . // (F.1.9) (F.3.6) Thus the Earth's torque on the satellite in Frame S' is, N'(b) = r1 x F1 + r2 x F2 = (GMEm1/r'13) r1b sinθ - (GMEm1/r'23) r1bsinθ = (GMEm1br1sinθ) [1/r'13 - 1/r'23] . (F.3.7) Note that this equation uses r1 and r2 and not r'1 and r'2 because the torque reference point is point b. We can confirm the two torque contributions by computing them geometrically from Fig (F.1.1) using these right-hand-rule helper drawings: (F.3.8) Then: r1 x F1 = r1F1 sin(π-β1) = r1F1sinβ1 = r1F1(bsinθ/r'1) // (F.1.6) then (F.3.2) = r1(GMEm1/r'12)(bsinθ/r'1) = (GMEm1/r'13) r1b sinθ // agrees with (F.3.5) r2 x F2 = r2F2 sin(π-β2) 2 = r2F2 sinβ2 2 = r2F2 sinβ2 [ - ] = r2F2 (bsinθ/r'2) [ - ] = r2(GMEm2/r'22) (bsinθ/r'2) [ - ] // (F.1.7) and (F.3.2) = r1(GMEm1/r'22) (bsinθ/r'2) [ - ] // (F.1.9) = - (GMEm1/r'23) r1bsinθ // agrees with (F.3.6) where we have used the fact that 2 = - 1 = - . There is a third way to compute the above torque, based on the torque theorem (D.1.3) which says that N(R) = N(0) - R x F . In our current context of working in Frame S' this reads N'(b) = N'(0) - b x F . (F.3.9) The torque N'(0) relative to the Frame S' origin is exactly 0 N'(0) = r' x F1 + r'2 x F2 = 0 + 0 = 0 since r'1 is collinear with F1 and r'2 is collinear with F2 as shown in Fig (F.1.1). We should include the stick tension/compression T, but (r'1 - r'2) x T = (r1 - r2) x T = 0 so T can be ignored. Thus we find from (F.3.9) that N'(b) = - b x (F1+ F2) = - b x[ - (GMEm1/r'13) r'1 - (GMEm2/r'23) r'2 ] = (GMEm1/r'13) b x r'1 + (GMEm2/r'23) b x r'2 = (GMEm1/r'13) (r'1-r1) x r'1 + (GMEm2/r'23) (r2' - r2) x r'2 // (F.1.8) = - (GMEm1/r'13) r1 x r'1 - (GMEm2/r'23) r2 x r'2 = - (GMEm1/r'13)[-r1bsinθ ] - (GMEm2/r'23) [ r2bsinθ ] // (F.3.4) then (F.1.9) = (GMEm1br1sinθ) [1/r'13 - 1/r'23] (F.3.10) in agreement with (F.3.7). Here are some quick checks on (F.3.10) : If the satellite is vertically aligned, θ = 0,π, then sinθ=0 and N'(b) = 0 as expected (no moment arms). If the satellite is horizontally aligned and m1= m2, then r'1 = r'2 so N'(b) = 0 (balanced moment arms). If m2> m1, then r2< r1 so for horizontal alignment one has r'2<r'1 so [ (1/r'13) - (1/r'23) ] < 0. Then if m1 is on the right, we have θ = π/2 and sinθ = 1 and then (F.3.10) has N(b) = -(positive). But for this orientation = 1 = - so N(b) = (positive) : (F.3.11) In this case the lever arms balance in the sense that m2r2 = m1r, but m2 is closer to Earth center so it feels the stronger force and so we expect N'(b) = (positive) . Far Approximation If we now assume as before that r1',r'2,b >> r1, r2 we can approximate [ (1/r1'3) - (1/r'23) ] by adding on to the Maple code shown in (D.3.5) to get where recall e1 = ε1 ≡ (r1/b). This time there is a leading linear term in ε1 so (1/r'13) - (1/r'23) ≈ 3 (r1/b) = - 3r1cosθ/ (μ2b4) . (F.3.12) Installing this result into (F.3.7) gives N'(b) = GMEm1b r1sinθ [ (1/r'13) - (1/r'23)] ≈ - GMEm1b r1sinθ ( 3r1cosθ/ (μ2b4) ) = - 3GME(m1/μ2) b-3r12sinθcosθ . (F.3.13) To this order of approximation, the torque N(b) vanishes when the satellite is horizontally aligned as well as when it is vertically aligned, and this is due to r1' ≈ r'2 in the horizontal case. F.4 The fictitious torque on the satellite in Frame S Recall the general fictitious torque expression given in (11.3.10), acting on a single particle of mass m, N'(c')fict = - (r'-c') x [ mS +mω x (ω x r') + 2m ω x v' + m x r' ] + m(' + ω x c' + S) x ( v' + ω x r' + S) – m' x v' . (11.3.10) Converted to swap notation this says, N(c)fict = - (r-c) x [ mS' +mω x (ω x r) + 2m ω x v + m x r ] + m( + ω x c + S') x ( v + ω x r + S') – m x v . (11.3.10)s (F.4.1) Here ω is the angular rotation rate of the satellite about the Earth. Our application has the torque center at c = 0 (and c' = b) so this simplifies somewhat to N(0)fict = - mr x [ S' +ω x (ω x r) + 2 ω x v + x r] + S' x ( mv + ω x [mr] ) . (F.4.2) frame centrifugal Coriolis Euler Recall from (8.1.8) that the square bracket in the above is - Ffict/m so we can trace the origin of the terms. Now write (F.4.2) separately for each of the two masses of the satellite : N(0)f,1 = - m1r1 x [ S' + ω x (ω x r1) + 2ω x v1 + x r1] + S' x ( m1v1 + ω x[m1r1]) N(0)f,2 = - m2r2 x [ S' + ω x (ω x r2) + 2ω x v2 + x r2] + S' x ( m2v2 + ω x [m2r2]) . (F.4.3) In the second line, use (F.1.9) to replace m2r2 = - m1r1 and (F.1.10) to replace m2v2 = - m1v1 : N(0)f,2 = + m1r1 x [ S' + ω x (ω x r2) + 2ω x v2 + x r2] + S' x ( -m1v1 + ω x [-m1r1]) . (F.4.4) Next, add the two torques to get the total fictitious torque on the satellite seen in Frame S, N(0)fict = N(0)f,1 + N(0)f,2 = - m1r1 x [ S' + ω x (ω x r1) + 2ω x v1 + x r1] + S' x ( m1v1 + ω x[m1r1]) + m1r1 x [ S' + ω x (ω x r2) + 2ω x v2 + x r2] + S' x ( -m1v1 + ω x [-m1r1]) = - m1r x { ω x (ω x [r1- r2]) + 2ω x [v1-v2] + x [r1- r2] } centrifugal Coriolis Euler = - m1r x { ω x (ω x [r1 + r1]) + 2ω x [v1 + v1] + x [r1 + r1] } = - m1 (1+) r1 x { ω x (ω x r1) + 2ω x v1 + x r1 } = - M r1 x { ω x (ω x r1) + 2ω x v1 + x r1 } . (F.4.5) All terms involving S' and S' have cancelled out. For the Earth orbit we have = 0 and ω = ω. Now consider this vector identity: C x [A x (A x C)] = C x [ (AC)A - A2C ] = (AC) C x A = – (AC) A x C . (F.4.6) Then r1 x [ ω x (ω x r1) ] = – (ωr1) (ω x r1) = - ω2r12 ( )( x ) // use (E.2.4) for and (E.2.13) for x = - ω2r12 ( sinθcosφ )(- cosθcosφ - sinφ ) = ω2r12 sinθcosφ (cosθcosφ + sinφ ) . (F.4.7) Next, the Coriolis term in (F.4.5) involves, r1 x (ω x v1) = (r1v1)ω - (r1ω)v1 // A x (B x C) = (AC)B - (AB)C . (F.4.8) But v1 is tangent to the radius-r1 sphere to which m1 is constrained, so (r1v1) = 0. Then r1 x (ω x v1) = - (r1ω)v1 = - r1ω ()v1 = - r1ω sinθcosφ v1 . (F.4.9) We now have this somewhat complicated expression for N(0)fict : N(0)fict = - (m1/m2)M r1 x [ ω x (ω x r1) + 2ω x v1 ] = - (m1/m2)M { r1 x [ω x (ω x r1)] + 2 r1 x (ω x v1) } = - (m1/m2)M { ω2r12 sinθcosφ [ cosθcosφ + sinφ ] - 2r1ω sinθcosφ v1 } (F.4.10) centrifugal Coriolis Now using (E.3.5), v1 = v1θ + v1φ = r1 + r1 sinθ (E.3.5) we then have, centrifugal Coriolis N(0)fict = - (m1/m2)M{ω2r12sinθcosφ [cosθcosφ + sinφ ] - 2r1ω sinθcosφ (r1 + r1 sinθ ) } = - (m1/m2)M { (ω2r12 sinθcosφsinφ - 2r12ωsinθcosφ) + (ω2r12 sinθcosφ cosθcosφ - 2r12ωsinθcosφsinθ ) } = - (m1/m2)M { ωr12 sinθ cosφ (ωsinφ -2) +ωr12sinθcosφ (ω cosθ cosφ - 2sinθ ) } = - (m1/m2)Mωr12sinθcosφ [ (ωsinφ -2) + (ω cosθ cosφ - 2sinθ ) ] . (F.4.11) How might one interpret this simple result? We examine the pieces of N(0)fict= 0 as they appear in (F.4.11). Terms involving velocities and arise from the Coriolis term, the other terms come from the centrifugal term. (a) the "frame effects" due to the motion of b (S' and S') are equal and opposite for the two masses because the origin of Frame S is at the center of mass causing m2r2 = - m1r1, as shown in (F.1.9). (b) within Frame S, the centrifugal acceleration term ω x (ω x r1) = -ω2r1 tries to push m1 to a larger radius. But r1 is constrained to lie on a sphere of radius r1 so m1 cannot go to a larger radius. This centrifugal acceleration is neutralized by part of the tension in the stick which we avoided talking about. This centrifugal term, by the way, involves "the short vector" r1 and not the long vector r'1. We discussed this situation in Section 8.2 for an Earth-based Frame S'. In our current context, the picture that corresponds to Fig (8.2.8) is the following (F.4.12) The arrow in each location during the orbit represents the position vector r1 of mass m1 where we assume that other effects are turned off so r1 stays fixed in Frame S. When these arrows are transferred to the picture on the right with common tails, we see that r1 does in fact go around in a circle of radius r1 and that is why the corresponding centrifugal force acting on m1 is -ω2r1 . (c) Within Frame S, the Coriolis force - 2m1 ω x v1 tries to deflect mass m1 "to the right" in Fig (F.1.1). But again, r1 is constrained to lie on a sphere of radius r1 so m1 cannot deflect to a different radius. This Coriolis force is neutralized by the rest of the tension in the stick. F.5 Equations of Motion for the satellite in Frame S (Spherical Coordinates) After much effort, we have arrived at this set of results for the satellite : L(0) = (m1/m2) M r12 ( – sinθ ) (F.2.2) (0) = (m1/m2) Mr12 [ (- 2 sinθ cosθ) – (2 cosθ + sinθ) ] . (F.2.3) N(b) = (GMEm1br1sinθ) (1/r'13 - 1/r'23) (F.3.7) r'12 = b2 + r12 + 2br1cosθ (F.1.4) r'22 = b2 + r22 - 2br2cosθ (F.1.5) N(0)fict = - (m1/m2)Mωr12 sinθcosφ [ (ωsinφ -2) + (ω cosθ cosφ - 2sinθ ) ] . (F.4.11) Using the orbit equation (8.6.7) applied to the satellite, GME = ω2b3 , // larger b means smaller ω (F.5.1) we can rewrite the true torque above as N'(b) = (ω2b4m1r1sinθ) (1/r'13 - 1/r'23) . (F.5.2) Writing (0) = N'(b) + N(0)fict then gives the satellite vector equation of motion, (m1/m2) Mr12 [ (- 2 sinθ cosθ) – (2 cosθ + sinθ) ] = (ω2b4m1r1sinθ) (1/r'13 - 1/r'23) - (m1/m2)Mωr12 sinθcosφ [ (ωsinφ -2) + (ω cosθ cosφ - 2sinθ ) ] . (F.5.3) Divide all three terms by the factor (m1/m2)Mr12. The coefficient of the first term on the right becomes (ω2b4m1r1sinθ)/ [ (m1/m2)Mr12] = (1/r1) (ω2b4(m2/M)sinθ) so the vector equation of motion is then, [ (- 2 sinθ cosθ) – (2 cosθ + sinθ) ] = (1/r1)(ω2b4μ2sinθ) (1/r'13 - 1/r'23) - ωsinθcosφ [ (ωsinφ -2) + (ω cosθ cosφ - 2sinθ ) ] . (F.5.4) Using now the far approximation (F.3.12) that (1/r'13) - (1/r'23) = -3r1cosθ/ (μ2b4), we get [ (- 2 sinθ cosθ) – (2 cosθ + sinθ) ] = -3ω2sinθcosθ - ωsinθcosφ [ (ωsinφ -2) + (ω cosθ cosφ - 2sinθ ) ] (F.5.5) As a reminder, the first line above is (0), the second line is true torque N'(b), and the last line is the fictitious torque N(0)fict created by the fact that Frame S is a rotating frame of reference, and we display the fictitious contributions in blue to keep track of them for a while below. Comment: Notice that the equation of motion is independent of m1 and m2 and hence of r1 and r2. If we were to vary the ratio m1/m2, we just "slide the stick" in Fig (F.1.1) so the Frame S origin remains at the center of mass point. The equation does depend on b through ω, since ω2 = GME/b3. Write out the separate component equations in the above, (- 2 sinθ cosθ) – (2 cosθ + sinθ) = -3ω2sinθcosθ - ωsinθcosφ (ωsinφ -2) - ωsinθcosφ (ω cosθ cosφ - 2sinθ ) or (F.5.6) (- 2 sinθ cosθ) – (2 cosθ + sinθ) +3ω2sinθcosθ + ωsinθcosφ (ωsinφ -2) + ωsinθcosφ (ω cosθ cosφ - 2sinθ ) = 0 . The component equations are then - 2 sinθ cosθ + 3ω2sinθcosθ + ωsinθcosφ (ω cosθ cosφ - 2sinθ ) = 0 // – (2 cosθ + sinθ) + ωsinθcosφ (ωsinφ -2) = 0 // where the fictitious torque terms are shown in blue. Changing the sign of the second equation and making a few adjustments we get, + sinθcosθ(3ω2 - 2) + ωsinθcosφ (ω cosθ cosφ - 2sinθ ) = 0 // + 2 cotθ – ωcosφ (ωsinφ -2) = 0 // (F.5.7) After much effort using the fictitious torque method we have finally arrived at the spherical equations of motion for the dumbbell satellite! These are two ordinary differential equations with time t as the variable. The equations are 2nd order in both θ(t) and φ(t) and they are non-linear due to terms like and 2. Finally, the equations are coupled, so they form a system of two 2nd order, coupled, non-linear ODE's. This system of two non-linear 2nd order ODE's can be trivially replaced with an equivalent system of four non-linear 1st order ODE's as follows, θ + sinθcosθ(3ω2 - 2) + ωsinθcosφ (ω cosθ cosφ - 2sinθ ) = 0 // φ + 2 cotθ – ωcosφ (ωsinφ -2) = 0 // = vθ = vφ (F.5.8) We mention this only because this is the first step Maple takes (internally) when it sets about solving a pair of 2nd order ODE's with its numerical dsolve command (coming soon). Right now, we make an anzats that (F.5.7) has a solution for which φ = π/2 (so cosφ = 0) and φ does not change, so that both and = 0 at all times. For such a solution, the dumbbell lies in the plane of paper of Fig (F.1.1) and the only Frame S motion is in the θ degree of freedom. In this case the two equations in (F.5.7) simplify to + 3ω2sinθcosθ = 0 // 0 = 0 . // in-plane // (F.5.9) For small θ, the first equation becomes + 3ω2θ = 0 which indicates sinusoidal oscillation in θ of the dumbbell about θ = 0 with frequency ωosc2 = 3ω2 so ωosc1 = ω Tosc1 = (1/) T ≈ 0.58 T // in-plane libration (F.5.10) where recall that ω is the orbital rotation rate of the satellite. For a low-Earth orbit with T = 88 minutes, the dumbbell initialized to φ = π/2 and a small angle θ would have an oscillation period of 88*.58 = 51 minutes. Next we look for a solution with φ = = 0 where the satellite at t = 0 is in a plane perpendicular to the plane of paper in Fig (F.1.1). In this case the equations (F.5.7) become, + 3ω2sinθcosθ + ωsinθ(ω cosθ) = 0 // sinθ – ωsinθ ( -2) = 0 // or +(2ω)2sinθcosθ = 0 // +2ω = 0 . // out-of-plane // (F.5.11) For small angles θ we have then +(2ω)2θ = 0 // = -2ω . // (F.5.12) The first equation implies θ oscillation at frequency ωosc2 = 2ω Tosc2 = 0.5T // out-of-plane libration (F.5.13) while the second equation shows that the dumbbell cannot remain very long at φ = 0 since ≠ 0, so this oscillation solution is only a temporary solution and the dumbbell is not stable in the plane φ = 0. Both libration frequencies appear on page 126 of Cosmo and Lorenzini with a reference (on C&L p 169) to the following item, 2. Beletskii, V. V. and Levin, E. M., "Dynamics of Space Tether Systems", Advances in the Astronautical Sciences , Vol. 83. (Univelt, Inc. 1993) Reader Exercise: Is the in-plane libration solution stable against small perturbations? F.6 Force analysis of the satellite in Frame S (Spherical Coordinates) Having obtained the dumbbell satellite θ,φ equations of motion in (F.5.7) using the fictitious torques method. we now set out to rederive these same equations using Newton's linear second law with fictitious forces. This time we also obtain an expression for the stick tension T which enables us to comment on the tidal forces for a static satellite positioned at θ = 0. Obtain the three equations of motion The only true forces on a dumbbell mass are gravity and stick tension T. From (F.3.2) we then write F'1 = - (GMEm1/r'13) r'1 - T 1 = - (GMEm1/r'13)( b + r1) - T 1 F'2 = - (GMEm2/r'23) r'2 - T 2 = - (GMEm2/r'23)( b + r2) - T 2 . (F.6.1) As noted earlier, the center of gravity is not quite at the center of mass in Fig (F.1.1), but the above equations are exact despite this fact. The forces are primed because they are forces in the inertial Frame S'. The reader is reminded that we are using the "swap notation" where prime↔noprime relative to the non-swap notation. In order to use Newton's Law in rotating Frame S, we must include the fictitious forces. We translate the result of (8.1.8) to swap notation to obtain Ffict,1 ≈ – m1S' – m1ω x (ω x r1) – 2m1 ω x v1 – m1 x r1 Ffict,2 ≈ – m2S' – m2ω x (ω x r2) – 2m2 ω x v2 – m2 x r2 . (F.6.2) frame centrifugal Coriolis Euler In these equations, the b acceleration is given by the swap version of (7.13) which then states S' = x b + ω x (ω x b) . // Special Case #1 (7.13)s (F.6.3) and the vectors b and ω are given as in Fig (F.1.1) by ω = ω b = b . (F.6.4) Finally we may state Newton's Law for each mass, Feff,1 = m1 a1 (F.6.5) ≈ - (GMEm1/r'13)r'1 - T 1 – m1S' – m1ω x (ω x r1) – 2m1 ω x v1 – m1 x r1 Feff,2 = m2a2 (F.6.6) ≈ - (GMEm2/r'23)r'2 - T 2 – m2S' – m2ω x (ω x r2) – 2m2 ω x v2 – m2 x r2 where we have now a set of six scalar equations. Using (F.1.9) through (F.1.11), (F.6.6) can be rewritten, Feff,2 = - m1a1 (F.6.7) ≈ - (GMEm2/r'23)r'2 + T 1 – m2S' + m1ω x (ω x r1) + 2m1 ω x v1 – m1 x r1 . Adding (F.6.5) and (F.6.7) gives 0 = - (GMEm1/r'13) r'1 - (GMEm2/r'23)r'2 - (m1+m2)S' or (m1+m2)S' = - (GMEm1/r'13) r'1 - (GMEm2/r'23)r'2 (F.6.8) This equation is just F = ma in inertial Frame S' for the total satellite where b is the center of mass. Ignoring the small offset between center of mass and center of gravity, the three equations (F.6.8) describe the circular orbit of the satellite around the Earth. We may then regard the equation (F.6.5) as a set of three scalar equations for the three unknowns θ,φ and T where recall r1 = (r1,θ,φ) in the spherical coordinates of Fig (F.1.1). Our next task is to write vector equation (F.6.5) in spherical coordinates to obtain the three equations of motion. After expanding the left side, we then consider the right side of (F.6.5) one term at a time: Feff,1 = m1 a1 (F.6.5) ≈ - (GMEm1/r'13)r'1 - T 1 – m1S' – m1ω x (ω x r1) – 2m1 ω x v1 – m1 x r1 1 2 3 4 5 6 Left side of (F.6.5): m1 a1 = m1(ar + aθ + aφ ) (E.3.6) = m1r1[( - 2 - 2 sin2θ) 1 + ( - 2 sinθ cosθ) + (2 cosθ + sinθ) ] (F.6.9) Term 1: - (GMEm1/r'13)r'1 = - (GMEm1/r'13)(b + r1) // (F.1.8) = - (GMEm1/r'13)(b + r1) // (F.6.4) = - (GMEm1/r'13)(bcosθ 1 - b sinθ + r1) // (E.2.7) = - (GMEm1/r'13)[ (bcosθ +r1) 1 - b sinθ ] (F.6.10) Term 2: - T 1 as is (F.6.11) Term 3: – m1S' = – m1 x b - m1 ω x (ω x b) // (F.6.3) = - m1ω x (ω x b) // satellite in circular orbit, = 0 = - m1(ωb)ω + m1ω2b //- A x (A x C) = -(AC)A + A2C = m1ω2b = m1ω2b // (F.6.4) = m1ω2b [ cosθ 1 - sinθ ] // (E.2.7) (F.6.12) Term 4: – m1ω x (ω x r1) = -m1(ωr1)ω + m1ω2r1 // identity shown above = -m1ω2r1(1) + m1ω2r1 1 // (F.6.4) = -m1ω2r1sinθcosφ + m1ω2r1 1 // (E.2.4)) = -m1ω2r1sinθcosφ[sinθcosφ 1 + cosθcosφ - sinφ ] + m1ω2r1 1 // (E.2.7) = - m1ω2r1 [ (sin2θcos2φ - 1) 1 + (sinθcosθcos2φ) + (- sinθcosφsinφ) (F.6.13) Term 5: -2m1 ω x v1 = -2m1 [ω] x ( vθ + vφ) // (E.3.5) = -2m1ω [vθ x + vφ x ] = -2m1ω [vθ (sinθcosφ + sinφ 1)+ vφ(-sinθ cosφ + cosθcosφ 1)] // (E.2.13) = -2m1ω [ (vθsinφ + vφcosθcosφ)1 + (-vφsinθ cosφ) + (vθsinθcosφ) ] = -2m1ωr1 [ ( sinφ + sinθ cosθcosφ)1 + (- sin2θ cosφ) + ( sinθcosφ) ] // (E.3.2) (F.6.14) Term 6: – m1 x r1 = 0 because we assume = 0 (F.6.15) Having all the bits and pieces, we now assemble the three component equations of (F.6.5). Feff,1 = m1 a1 (F.6.5) ≈ - (GMEm1/r'13)r'1 - T 1 – m1S' – m1ω x (ω x r1) – 2m1 ω x v1 – m1 x r1 1 2 3 4 5 6 1: m1r1( - 2 - 2 sin2θ) = - (GMEm1/r'13)(bcosθ +r1) - T + m1ω2b cosθ - m1ω2r1(sin2θcos2φ - 1) 1 2 3 4 - 2m1ωr1 ( sinφ + sinθ cosθcosφ) (F.6.16) 5 : m1r1 ( - 2 sinθ cosθ) = + (GMEm1/r'13) b sinθ - m1ω2bsinθ 1 3 - m1ω2r1 sinθcosθcos2φ + 2m1ωr1 sin2θcosφ (F.6.17) 4 5 : m1r1(2 cosθ + sinθ) = + m1ω2r1sinθcosφsinφ - 2m1ωr1( sinθcosφ) (F.6.18) 4 5 We now rewrite the three equations dividing by m1 and using (F.5.1) that GME = ω2b3 : 1: r1( - 2 - 2 sin2θ) = - ( ω2b3/r'13)(bcosθ +r1) - T/m1 + ω2b cosθ - ω2r1(sin2θcos2φ - 1) - 2ωr1 ( sinφ + sinθ cosθcosφ) (F.6.19) : r1 ( - 2 sinθ cosθ) = + (ω2b3/r'13) b sinθ - ω2bsinθ - ω2r1 sinθcosθcos2φ + 2ωr1 sin2θcosφ (F.6.20) : (2 cosθ + sinθ) = + ω2sinθcosφsinφ - 2ω( sinθcosφ) (F.6.21) If one uses (F.1.4) that r'12 = b2 + r12 + 2br1cosθ in (F.6.20), the pair of equations (F.6.20) and (F.6.21) can in theory be solved for θ(t) and φ(t), given appropriate initial conditions. The solutions can then be inserted into (F.6.19) to obtain a result for the stick tension T(t). Verify the angular equations of motion We can rewrite(F.6.21) as : sinθ + 2 cosθ - ωsinθcosφ (ωsinφ - 2) (F.6.22) which matches the torque equation (F.5.7), sinθ + 2 cosθ – ωsinθcosφ (ωsinφ -2) = 0 . // (F.5.7) Next, the equation (F.6.20) may be rewritten r1 ( - 2 sinθ cosθ) = + (ω2b3/r'13) b sinθ - ω2bsinθ - ω2r1 sinθcosθcos2φ + 2ωr1 sin2θcosφ or r1 ( - 2 sinθ cosθ) = +ω2bsinθ [(b/r'1)3 - 1] - ω2r1 sinθcosθcos2φ + 2ωr1 sin2θcosφ or - 2 sinθ cosθ = +ω2bsinθ [(b/r'1)3 - 1]/r1 - ω2sinθcosθcos2φ + 2ω sin2θcosφ . (F.6.23) Now assume the far approximation where r'1, b >> r. Recall that r'12 = b2 + r12 + 2br1cosθ (r'1/b)2 = 1 + (r1/b)2 + 2(r1/b)cosθ (r'1/b)3 = [ 1 + (r1/b)2 + 2(r1/b)cosθ ]3/2 (b/r'1)3 = [ 1 + (r1/b)2 + 2(r1/b)cosθ ]-3/2 ≈ 1 + (-3/2) [(r1/b)2 + 2(r1/b)cosθ ] ≈ 1 + (-3/2) 2(r1/b)cosθ = 1 - 3(r1/b)cosθ so [(b/r'1)3 - 1] ≈ - 3(r1/b)cosθ . (F.6.24) Then (F.6.23) becomes - 2 sinθ cosθ = +ω2bsinθ [- 3(r1/b)cosθ]/r1 - ω2sinθcosθcos2φ + 2ω sin2θcosφ or - 2 sinθ cosθ = -3ω2sinθcosθ - ω2sinθcosθcos2φ + 2ω sin2θcosφ or - 2 sinθ cosθ + 3ω2sinθcosθ + ω2sinθcosθcos2φ -2ω sin2θcosφ = 0 or + 3ω2sinθcosθ - 2 sinθ cosθ + ωsinθcosφ(ωcosθcosφ - 2sinθ) = 0 (F.6.25) which matches the torque equation (F.5.7), + 3ω2sinθcosθ - 2 sinθ cosθ + ωsinθcosφ (ω cosθ cosφ - 2sinθ ) = 0 // (F.5.7) At this point we have derived the same angular equations of motion (F.5.7) in two different ways. Obtaining the tension in the stick (or tether) Finally we come to the 1 equation (F.6.19), 1: r1( - 2 - 2 sin2θ) = - ( ω2b3/r'13)(bcosθ +r1) - T/m1 + ω2b cosθ - ω2r1(sin2θcos2φ - 1) - 2ωr1 ( sinφ + sinθ cosθcosφ) . (F.6.19) Inserting the far approximation (F.6.24) that (b/r'1)3 ≈ 1 - 3(r1/b)cosθ into the above gives r1( - 2 - 2 sin2θ) = -ω2[1 - 3(r1/b)cosθ] (bcosθ +r1) - T/m1 + ω2b cosθ - ω2r1(sin2θcos2φ - 1) - 2ωr1 ( sinφ + sinθ cosθcosφ) r1( - 2 - 2 sin2θ) = -ω2(bcosθ +r1) + 3ω2(r1/b)cosθ(bcosθ +r1) - T/m1 + ω2b cosθ - ω2r1(sin2θcos2φ - 1) - 2ωr1 ( sinφ + sinθ cosθcosφ) r1( - 2 - 2 sin2θ) = -ω2bcosθ - ω2r1 + 3ω2r1cosθ(cosθ +r1/b) - T/m1 + ω2b cosθ - ω2r1sin2θcos2φ + ω2r1 - 2ωr1 ( sinφ + sinθ cosθcosφ) r1( - 2 - 2 sin2θ) = 3ω2r1cosθ(cosθ +r1/b) - T/m1 - ω2r1sin2θcos2φ - 2ωr1( sinφ + sinθ cosθcosφ) ( - 2 - 2 sin2θ) ≈ 3ω2cos2θ - T/(m1r1) - ω2sin2θcos2φ - 2ω( sinφ + sinθ cosθcosφ) ( - 2 - 2 sin2θ) ≈ ω2(3cos2θ -sin2θcos2φ) - T/(m1r1) - 2ω( sinφ + sinθcosθcosφ) . The stick tension is then given by (this will later be verified in Section F.9), T = m1r1[ ω2(3cos2θ - sin2θcos2φ) + 2 + 2 sin2θ - 2ω( sinφ + sinθcosθcosφ) ] . (F.6.26) If the angular velocities are very small such that | | << ω and | | << ω, the result becomes T ≈ m1ω2r1(3cos2θ - sin2θcos2φ) . // small velocities (F.6.27) In Cartesian coordinates this becomes, T ≈ (m1ω2/r1)(3r12cos2θ- r12sin2θcos2φ) = (m1ω2/r1)(3z2-x2) . // small velocities (F.6.28) In this small-velocity limit, the tension in the stick is positive as long as |x| < |z| T > 0 // small velocity limit (F.6.29) so for small x displacements T is always positive. As noted, in general one must solve (F.6.20) and (F.6.21) for θ(t) and φ(t), then (F.6.26) gives T(t). Tidal Force If the dumbbell is static at θ = 0, we see from (F.6.26) or (F.6.27) that T = 3m1ω2r1 (F.6.30) which we associated with a "tidal force". The factor of 3 arises from (F.6.24) which in turn arises from the power 3 in the gravitational force factor in (F.6.1), (GMEm1/r'13) = (ω2b3 m1/r'13) = ω2m1 (b/r'1)3 . It happens that in Frame S' the gradient of the radial gravitational field at mass m1 is ∂r'(-GMEm1/r'12) = ∂r1'(-ω2b3m1/r'12) = -ω2b3m1 ∂r1' (r'1-2) = 2ω2b3m1/r'13 = 2ω2m1(b/r'1)3 so one can associate the factor of 3 in (F.6.30) with this gradient. However, the result (F.6.30) really comes from a sum of several terms in (F.6.5) for a static dumbbell, as we now review (the 1 equation): Feff,1 = m1 a1 (F.6.5) ≈ - (GMEm1/r'13)r'1 - T 1 – m1S' – m1ω x (ω x r1) – 2m1 ω x v1 – m1 x r1 1 2 3 4 5 6 1: m1r1( - 2 - 2 sin2θ) = - (GMEm1/r'13)(bcosθ +r1) - T + m1ω2b cosθ - m1ω2r1(sin2θcos2φ - 1) 1 2 3 4 - 2m1ωr1 ( sinφ + sinθ cosθcosφ) (F.6.16) 5 m1r1( 0 - 0) = - (GMEm1/r'13)(b +r1) - T + m1ω2b - m1ω2r1(0 - 1) - 2m1ωr1 (0 + 0) 1 2 3 4 5 0 = - m1(ω2b3/r'13)(b +r1) - T + m1ω2b + m1ω2r1 1 2 3 4 0 = {- m1ω2[1 - 3(r1/b)](b +r1)} - T + m1ω2b + m1ω2r1 1 2 3 4 0 = { - m1ω2 - m1ω2b +3m1ω2r1 } - T + m1ω2b + m1ω2r1 1 2 3 4 0 = 3m1ω2r1 - T 1,3,4 2 Thus in rotating Frame S the expression (F.6.30) for the tidal force T has contributions from the gravitational gradient (term 1) and from the "frame" term – m1S (term 3) which is the centrifugal contribution due to the acceleration of Frame S toward the earth, and finally from the "local centrifugal term" – m1ω x (ω x r1) (term 4). F.7 Numerical solutions of the equations of motion (Spherical Coordinates) The angular equations of motion for mass m1 of the dumbbell (or tether) satellite are stated in (F.5.7) which we replicate here, + sinθcosθ(3ω2 - 2) + ωsinθcosφ (ω cosθ cosφ - 2sinθ ) = 0 + 2 cotθ – ωcosφ (ωsinφ -2) = 0 . (F.5.7) (F.7.1) The first step is to enter these two equations into Maple: (F.7.2) For illustration purposes we have set the satellite orbit frequency to ω = 1 sec-1 so ω = 1 T = 2π = 6.28 sec Verification of in-plane libration The initial conditions are taken to be (see Fig (F.1.1) to see that φ = π/2 is the in-plane situation ) θ = 0.2 φ = π/2 = 0 = 0.01 . (F.7.3) Here θ = 0.2 = 11.5o is a fairly small angle. We now call Maple's ODE solver routine dsolve and plot 40 times θ(t) in red and Mod2Pi(φ) in black. The peaks of the red curve are thus at 40*.2 = 8. Without this mod routine, the black φ curve just winds up without limit. (F.7.4a) red = 40*θ black = Mod2p(φ) Rather than plot a cosine wave, dsolve keeps θ positive all the time and has φ jump by π each time the solution passes through θ = 0. In spherical coordinates the figure shows the expected θ curve! In order to avoid exactly hitting θ = 0 which is a singular point in spherical coordinates (φ is undefined there), we have added a small = .001 to cause the dumbbell to slightly miss the z axis. One can see from the second equation in (F.7.1) that the numerical integrator is faced with = 2 cotθ + stuff, and cotθ blows up at θ = 0 (and generates an error message in odeplot). We can blow up the region t = (2,4) : (F.7.4b) From (F.5.10) for in-plane libration one predicts, Tocs1 = T/ = 6.28/1.73 = 3.63 sec (F.7.5) and this value is verified by the vertical line in the above figure (each tick is .04) Verification of out-of-plane libration The initial conditions are now taken to be θ = 0.2 φ = 0 = 0 = 0 (F.7.6) We then rerun the above code with a different set of "inits" : (F.7.7a) Now the swing misses θ = 0 of its own accord. One can see that the black φ curve starts moving away from φ = 0 at about t = 0.5 and in general the black φ curve has smooth rises near small θ. We again blow up the region t = (2,4) : (F.7.7b) From (F.5.13) for out-of-plane libration one predicts, Tocs2 = T/2 = 6.28/2 = 3.14 (F.7.8) while the figure shows about 3.16, close enough. Maple can also plot (t) and (t) and here is a plot showing all four curves : (F.7.9) red = 10*θ black = φ green = 10* blue = When red θ nears θ = 0, blue has major action since φ is quickly changing by π. And the green of spherical coordinates has to make a radical change since it is in effect suddenly reversing course. A more general solution Here in addition to starting with θ(0) = 0.2 and φ(0) = 0, we provide a push in the azimuthal direction so that mass m1 of the dumbbell satellite then swings around in azimuth while it oscillates in θ : θ = 0.2 φ = 0 = 0 = 2 . (F.7.10) (F.7.11) Notice that the azimuthal velocity slows down near the peaks of red θ(t). Energy is transferred back and forth between the θ and φ degrees of freedom in this system. One can use the same odeplot routine used above to make "orbital" plots in angle space, (F.7.12) Rather than produce more plots here, we shall defer to Section F.9 where we shall plot x(t) and y(t) instead of θ(t) and φ(t). F.8 Force analysis of the satellite in Frame S (Cartesian Coordinates) Section F.6 developed equations of motion for the dumbbell satellite in spherical coordinates. Here we repeat that development but in Cartesian coordinates where things are in many ways simpler. We follow Section F.6 down Newton's Law (F.6.5) , Feff,1 = m1 a1 (F.6.5) (F.8.1) ≈ - (GMEm1/r'13)r'1 - T 1 – m1S' – m1ω x (ω x r1) – 2m1 ω x v1 – m1 x r1 . 1 2 3 4 5 6 As a reminder, this is Newton's Law (Feff = ma) for mass m1 of the satellite in Frame S where fictitious forces are included. Although this mass has coordinate r1 in Frame S, we shall refer to its components without 1 subscripts, so r1 = (x,y,z). This is similar to how we used r1 = (r1,θ,φ) in spherical coordinates where θ and φ had implied "1" subscripts. We also write v1 = v and a1 = a to reduce clutter. The quantity T is the tension in the stick (or tether). We now evaluate the six terms of (F.8.1) in Cartesian coordinates, mimicking (F.6.9) through (F.6.14): Left side of (F.8.1): m1a1 = m1a = m1(ax + ay + az ) (F.8.2) Term 1: - (GMEm1/r'13)r'1 = - (GMEm1/r'13)(b + r1) = - (GMEm1/r'13) [ b + x + y + z ] = - (GMEm1/r'13) [ x + y + (b+z) ] (F.8.3) Term 2: - T 1 = -(T/r1)r1 = -(T/r1) [ x + y + z] (F.8.4) Term 3: – m1S' = – m1 x b - m1 ω x (ω x b) // (F.6.3) = - m1ω x (ω x b) // satellite in circular orbit, = 0 = - m1(ωb)ω + m1ω2b //- A x (A x C) = -(AC)A + A2C = m1ω2b = m1ω2b (F.8.5) Term 4: – m1ω x (ω x r1) = -m1(ωr1)ω + m1ω2r1 // identity shown above = -m1ω2(r1) + m1ω2r1 = -m1ω2(x) + m1ω2 [ x + y + z] = m1ω2 [ y + z] (F.8.6) Term 5: -2m1 ω x v1 = -2m1 [ω] x ( vx + vy + vz) = -2m1ω [ vy - vz] (F.8.7) Term 6: – m1 x r1 = 0 because we assume = 0 (F.8.8) Having all the bits and pieces, we now assemble the three component equations of (F.8.1). The numbers show the Term above associated with each piece: Feff,1 = m1 a (F.8.1) ≈ - (GMEm1/r'13)r'1 - T 1 – m1S' – m1ω x (ω x r1) – 2m1 ω x v1 – m1 x r1 1 2 3 4 5 6 : m1 ax = - (GMEm1/r'13) x - (T/r1)x 1 2 : m1 ay = - (GMEm1/r'13) y - (T/r1)y + m1ω2y + 2m1ωvz 1 2 4 5 : m1 az = - (GMEm1/r'13) (b+z) - (T/r1)z + m1ω2b + m1ω2z - 2m1ω vy (F.8.9) 1 2 3 4 5 We now rewrite the three equations dividing by m1 and using (F.5.1) that GME = ω2b3 . At the same time we replace velocity and acceleration components with dot notation components like and : = - [(ω2b3/r'13) + (T/m1r1)]x = - [(ω2b3/r'13) + (T/m1r1) - ω2]y + 2ω = - [(ω2b3/r'13) - ω2] (b+z) - (T/m1r1)z - 2ω x2+y2+z2 = r12 (F.8.10) where r'12 = (r1+b)2 = r12 + b2 + 2 r1 b = r12 + b2 + 2 r1 [b] = r12 + b2 + 2bz . (F.8.11) Eq. (F.8.10) is a system of 4 equations in 4 unknowns x,y,z,T. To simplify manipulations below, define A ≡ (ω2b3/r'13) B ≡ (T/m1r1) // rescaled tension so the system of equations becomes 1 = - (A + B)x 2 = - (A + B - ω2)y + 2ω 3 = - (A - ω2) (b+z) - Bz - 2ω 4 x2+y2+z2 = r12 (F.8.12) We now wish to eliminate the rescaled tension B from the equation set. We first eliminate B between equations 1 and 2 : 1*y y = - (A + B)xy 2*x x = - (A + B - ω2)xy + 2ωx . Subtract so that the -(A + B)xy terms cancel, y - x = -ω2xy - 2ωx . (F.8.13) Next, we eliminate B between equations 1 and 3: 1*z z = - (A + B)xz 3*x x = - (A - ω2)(b+z)x - Bxz - 2ωx = -Abx -Azx +ω2(b+z)x - Bxz - 2ωx = - (A + B)xz - Abx +ω2(b+z)x - 2ωx . Subtract so that the -(A + B)xy terms cancel, z - x = Abx -ω2(b+z)x + 2ωx . (F.8.14) We now have a system of 3 equations in three unknowns x,y,z, where we now restore A = (ω2b3/r'13) 1 y - x = - ω2xy - 2ωx 2 z - x = (ω2b3/r'13)bx - ω2(b+z)x + 2ωx where r'12 = r12 + b2 + 2bz 3 x2+y2+z2 = r12 . (F.8.15) For convenience, we now reorder and rename equations 1 and 2 of this set as follows: eq3 z - x = (ω2b3/r'13)bx - ω2(b+z)x + 2ωx eq2 y - x = - ω2xy - 2ωx . (F.8.16) Now apply the far approximation in equation eq3. From (F.6.24) we know that [(b/r'1)3 - 1] ≈ - 3(r1/b)cosθ . (F.6.24) so (b/r'1)3 ≈ 1 - 3(z/b) // z = r1cosθ and A ≡ (ω2b3/r'13) ≈ ω2[1 - 3(z/b) ] . (F.8.17) Equation eq3 above then becomes, z - x = (ω2b3/r'13)bx - ω2(b+z)x + 2ωx = ω2[1 - 3(z/b) ]bx - ω2(b+z)x + 2ωx = ω2bx - 3(z/b)ω2bx - ω2bx - ω2zx + 2ωx = - 3(z/b)ω2bx - ω2zx + 2ωx = - 3ω2xz - ω2zx + 2ωx = - 4ω2xz + 2ωx . (F.8.18) This in the far approximation the equations of motion of mass m1 of the dumbbell are eq3 z - x = - 4ω2xz + 2ωx eq2 y - x = - ω2xy - 2ωx x2+y2+z2 = r12 . (F.8.19) Once these three equations are solved for x(t), y(t) and z(t), we can find the tension T(t) from (say) the first equation of (F.8.12) : /x +A+B = 0 B = - /x - A (T/m1r1) = - /x - (ω2b3/r'13) so T = - m1r1[ /x + (ω2b3/r'13) ] so T = - m1r1[ /x + ω2 ] . // far approximation (F.8.17) (F.8.20) F.9 Verification of the Cartesian equations of motion and stick tension In Section C.7 we showed that the x,y,z Foucault pendulum equations of motion were the same as the θ,φ ones. Here we repeat that task for the dumbbell satellite equations of motion. We do these verifications to strengthen our confidence in all the equations since there are not many external sources for verification. A trusting reader can just skip this section. First, here are the satellite angular equations of motion from (F.5.7), eq1 + sinθcosθ(3ω2 - 2) + ωsinθcosφ (ω cosθ cosφ - 2sinθ ) = 0 eq2 + 2cotθ – ωcosφ (ωsinφ -2) = 0 . (F.5.7) (F.9.1) These were derived in (F.5.7) using the effective torque method and were then verified using the effective force method in (F.6.22) and (F.6.25). Meanwhile, here are the Cartesian equations of motion from Section F.8, eq3 z - x = - 4ω2xz + 2ωx eq2 y - x = - ω2xy - 2ωx . (F.8.19) (F.9.2) Below we shall show that Task (a): eq2 of (F.9.2) eq2 of (F.9.1) Task (b): [(sinθ)*eq3 - (cosθsinφ)*eq2] of (F.9.2) eq1 of (F.9.1) That is to say, angular eq1 of (F.9.1) is a certain linear combination of eq3 and eq2 of (F.9.2). If we can show Task (a) and Task (b) above, then we have shown that (F.9.2) (F.9.1), and this then serves as verification of (F.9.2). Maple must replace x,y,z and derivatives with r1,θ,φ and derivatives. For coordinates and first derivatives, (F.9.3) The second derivatives are messier, but Maple is happy to do the calculations, (F.9.4) Task (a): Show that eq2 of (F.9.2) eq2 of (F.9.1) We enter eq2 of (F.9.2) and do some manipulations, suppressing the output except for the last step : (F.9.5) On the first red code line we enter eq2 of (F.9.2) and then divide the result by r12. We then replace occurrences of cos2θ by 1-sin2θ. We use lhs = "left hand side" so we end up only with the left side of an equation which says stuff = 0. Symbol % refers to the last computed quantity. The blue result can be manually transcribed as 2sinθcosθ + sin2θ - ω2sin2θ cosφsinφ + 2ωsin2θ cosφ = 0 . Now divide by sin2θ and reorder the four terms to get + 2cot(θ) - ω2cosφsinφ + 2ωcosφ = 0 or + 2cot(θ) - ωcosφ(ωsinφ - 2) = 0 (F.9.6) This is a match for eq2 of (F.9.1) so we have accomplished Task (a). Task (b): (sinθ)eq3 - (cosθsinφ)eq2 of (F.9.2) eq1 of (F.9.1) The code continues from that shown above. Equation eq2 is already entered, so we now enter eq3, form the linear combination for eq1, then process the results with a series of typical tortuous Maple steps, (F.9.7) We again manually transcribe the result -sinθcosθ 2 - 2ωsin2θcosφ +4ω2sinθcosθ - ω2sinθcosθsin2φ + or -sinθcosθ 2 - 2ωsin2θcosφ +3ω2sinθcosθ +ω2sinθcosθ - ω2sinθcosθsin2φ + or + sinθcosθ 3ω2 -sinθcosθ 2 + ω2sinθcosθ - ω2sinθcosθsin2φ - 2ωsin2θcosφ or + sinθcosθ( 3ω2 - 2) + ω2sinθcosθ(1-sin2φ) - 2ωsin2θcosφ or + sinθcosθ( 3ω2 - 2) + ω2sinθcosθcos2φ - 2ωsin2θcosφ or + sinθcosθ( 3ω2 - 2) + ωsinθcosφ (ω cosθ cosφ - 2sinθ ) (F.9.8) and after "pulling teeth" we do end up with eq1 of (F.9.1) so Task (b) is accomplished. Tension equation verification Using the angular equations of motion (F.9.1), we now show that the following two tension expressions are the same (the first is angular (F.6.26) while the second is Cartesian (F.8.20)) , T/(m1r1) = ω2(3cos2θ - sin2θcos2φ) + 2 + 2 sin2θ - 2ω( sinφ + sinθcosθcosφ) (F.6.26) T/(m1r1) = - [ /x + ω2 ] . (F.8.20) (F.9.9) Our task of showing (F.9.9) is the same as showing that x [ω2(3cos2θ - sin2θcos2φ) + 2 + 2 sin2θ - 2ω( sinφ + sinθcosθcosφ) ] = - - ω2x or x [ω2(3cos2θ - sin2θcos2φ + 1) + 2 + 2 sin2θ - 2ω( sinφ + sinθcosθcosφ) ] = - or LHS = RHS . (F.9.10) We first get the complicated left hand side LHS entered: (F.9.11) We then compute RHS = - as done earlier in this section, (F.9.12) Notice that RHS contains second derivatives and . We shall eliminate these derivatives by manually solving the angular equations of motion (F.9.1) for Tdd = and Pdd = : To show that LHS = RHS, we define d = LHS-RHS and show that d = 0: (F.9.13) (F.9.14) Thus d = 0 and LHS = RHS and the two expressions for T in (F.9.9) are the same. F.10 Numerical solutions of the equations of motion (Cartesian Coordinates) We have done a lot of "work" in this Appendix F, and now it is time to "play", making use of our hard-won Cartesian equations of motion which don't have the singularity problems had by the angular equations at θ = 0. Our task is to solve the set of equations (F.8.19) (eq1 now has a new meaning) : eq1 x2+y2+z2 = r12 eq2 y - x = - ω2xy - 2ωx eq3 z - x = - 4ω2xz + 2ωx (F.8.19) (F.10.1) These equations describe the motion of mass m1 of the dumbbell satellite in rotating Frame S as depicted in Fig (F.1.1). The position of mass m1 is (x,y,z) where x2+y2+z2 = r12. The motion of mass m2 is then determined by (D.2.8) m1r1 = - m2r2 so (x2,y2,z2) = - (m1/m2)(x,y,z). The length of the stick of the dumbbell satellite is s = r1+ r2 = r1+ (m1/m2)r1 = [1 + (m1/m2] r1. We enter eq2 and eq3 writing derivatives for example as = xdd (w = ω) , (F.10.2) At this point zd = and zdd = are unspecified. We use eq1 to compute and in terms of x and y, using eq1 above: (F.10.3) When these expressions are installed, eq2 and eq3 becomes these formidable-looking equations which contain two unknown functions x(t) and y(t) and constants r1 and ω : (F.10.4) The reader is reminded of the geometry of Fig (F.1.1) where z points up, away from Earth center, y points to the right and is in the plane of the satellite orbit, while x is perpendicular to the plane of the satellite. (F.10.5) Our plots below in the (x,y) plane are what a viewer would see looking at mass m1 "from above", that is, from a point at perhaps z = b+2s on the z axis in the above figure. In-Plane Libration As our first test, we shall look for the "in-plane libration" behavior. We examined this behavior earlier below (F.7.3) in angular coordinates, and we now look in Cartesian coordinates. The initial conditions are: x(0) = 0 y(0) = 1 (0) = 0 (0) = 0 (F.10.6) With ω = 1, we expect to get a simple swinging back in the x=0 plane with period T = 3.63 sec as shown in (F.7.5). A half period is then 1.82 seconds. The Maple code to invoke a solution is as follows (for a numerical integration from t = 0 to t = 1.82 sec): (F.10.7) We show the result below on the left, and then from t = 0 to t = 0.91 (quarter period) on the right : (F.10.8) Thus both the "orbit" and the period for in-plane libration are visually confirmed. If we run from t = 0 to t = 10, the graph is as on the left above since mass m1 just swings back an forth in the same orbit, never leaving the y-axis. Out-of-Plane Libration We examined this behavior earlier below (F.7.6) in angular coordinates, and we now look in Cartesian coordinates. The initial conditions are now, x(0) = 1 y(0) = 0 (0) = 0 (0) = 0 (F.10.9) With ω = 1, we expect to get a swinging back in the y=0 with period T = 3.14 sec as shown in (F.7.8). Here is what Maple has to say: (F.10.10) The period looks right since mass m1 swings back close to its initial position after 3.14 seconds, but one sees that the motion is not quite in the y=0 plane, so out-of-plane libration is an approximate concept as we noted earlier below (F.5.13). Note in the figure the fine scale of the vertical axis relative to horizontal. Here are orbits for a selection of final integration times: t = 4.1 t = 8.2 t = 15.9 t = 21 t = 41.6 (F.10.11) It does seem that the out-of-plane libration stays within a certain small band of deviation in the y direction which we shall leave to the reader to theoretically calculate. In each plot the time was selected to make the tail of the trace clearly visible. The author is reminded of a lecture given by Shelly Glashow on the question: What can be said about orbits on an arbitrarily-shaped-but-convex billiard table? Do they all eventually close on themselves, or might some never close? (Exercise for the reader). As regards the above plots, the fact that the ratio of the two libration frequencies is an irrational number /2 might have some bearing on the closure of the orbits. In the examples above r1 = 10 and we have used x(0) = 1 or y(0) = 1 to obtain "small oscillation". If in the in-plane libration case we use y(0) = 9, there is no change in the orbit, but the oscillation period is slightly altered. If in the out-of-plane libration case we use x(0) = 9, the orbit no longer maintains the narrow band as in the above examples. For example, going again to t = 64 seconds with x(0) = 9, (F.10.12) The Digits := 14 command tells Maple to compute the numerical integration with 14 decimal places of accuracy instead of the default 10 digits. Starting with a diagonal initial position x(0) = 1 y(0) = 1 (0) = 0 (0) = 0 t = 8 t = 64 (F.10.13) The t=64 result is reminiscent of a Lissajous pattern on an oscilloscope screen when the x and y axes are driven by different frequency sine ways. (See sine plots below.) In some sense these are the two libration frequencies. Attempting a Circular Orbit (Conical Solution) We have made many attempts to get a circular-like orbit by giving the mass m1 an initial velocity kick in some useful direction, but this system does not want to cooperate. Here is an example : (F.10.14) What starts as a rough circle is soon distorted into a narrow orbit. In the case of the spherical pendulum we had a Conical Solution in (C.5.20) where θ = θ0 and = constant. Assuming θ = θ0 in the satellite angular equations (F.5.7) gives sinθ0cosθ0(3ω2 - 2) + ωsinθ0cosφ (ω cosθ0 cosφ - 2sinθ0 ) = 0 // – ωcosφ (ωsinφ) = 0 . // Since this is two ODE's for the one function φ(t), it seems unlikely there is any general non-static solution. If = 0 the equations become sinθ0cosθ0(3ω2) + ωsinθ0cosφ (ω cosθ0 cosφ) = 0 – ωcosφ (ωsinφ) = 0 . If φ = 0 the first equation requires that θ0 = 0 or π/2 which are static vertical and horizontal positions. The same is true for φ = π/2 . Three-dimensional plots To make such plots, one must first extract the solution functions from the dsolve environment. For details on how this works and other information on dsolve (including a debugger's guide), see the author's Maple User Guide. Here we extract the functions calling them X,Y and Z , (F.10.15) The following code then creates a 3D orbit and superposes it on a contour sphere of radius r1 , (F.10.16) For a sample application, we start mass m1 at the north pole and give it a good kick in the x and y directions with (0) = 10 and (0)= 10 to get an x,y plot : ( (F.10.17) Here then is the corresponding 3D plot (F.10.18) where the sphere of radius r1 is gradually tipped down toward the viewer. Conventional plots In order to plot x(t), y(t) and so on, we first extract all functions from the dsolve system and then crudely add missing pieces like Xdd: (F.10.19) Here then is a plot of x(t),y(t),z(t) = red,black,blue for the above example: red = x(t) black = y(t) blue = z(t) (F.10.20) Here one sees red x and black y executing roughly sinusoidal motions. In a ballpark sense x and y are sinusoidal at their respective libration periods (3.14 and 3.63), and at least we see that black y has a longer period than red x (this black y period seems more like ~4). This is what creates the Lissajous pattern in our earlier figures. Meanwhile, blue z(t) is not coming down much from its maximum value of z = 10. Here is a plot of x,, = red.black,blue : red = x(t) black = (t) blue = (t) (F.10.21) Notice that red x and blue always cross the axis at the same time, which allows /x to be finite at all values of t (see below). Tension in the the stick or tether The tension in the stick (tether) is stated in (F.8.20), T = - m1r1[ /x +ω2 ], (F.8.20) which we then plot with m1 = 1 for the above example : (F.10.22) If we lower (0) = (0) = 10 to (0) = (0) = 1 to get a milder motion, the tension is less variable, (F.10.23) Here T is roughly equal to the "DC" value for the static dumbbell with θ = 0. From (F.6.30) that value is T ≈ 3m1ω2r1 = 3 * 1 * 12 * 10 = 30 N . // tether tension (F.10.24) We have been using ω = 1, but for a low-Earth orbit satellite one has Torbit ≈ 88*60 seconds so ω = 2π/Torbit ≈ .0012 . (F.10.25) Then mass m1 = 1 kg on a 20 meter tether with equal masses (r1 = 10 m) would have a tidal force of T ≈ 3m1ω2r1 = 3 * 1 * (.0012)2 * 10 = .4320e-4 N = 43 μN (F.10.26) which is a very small tension.