General Disclaimer

One or more of the Following Statements may affect this Document

This document has been reproduced from the best copy furnished by the 

organizational source. It is being released in the interest of making available as 

much information as possible.

This document may contain data, which exceeds the sheet parameters. It was 

furnished in this condition by the organizational source and is the best copy 

available.

This document may contain tone-on-tone or color graphs, charts and/or pictures, 

which have been reproduced in black and white.

This document is paginated as submitted by the original source.

Portions of this document are not fully legible due to the historical nature of some 

of the material. However, it is the best reproduction available from the original 

submission.

Produced by the NASA Center for Aerospace Information (CASI)
G ,

The '1nvo*Ugotioe o>t Tethered

s ca11tsi Rate. DVnattws ` it

Contract *A88-36160

Quarterly Report A2

For the period 1 December 1984 through 14 . February 1985

Principal .Investigator

Dr. Enrico torenaini

(la6^S-C.-171^_^) .i1. :t Vi_" l:.C., G8^-2 3

^'^altELl} r^^clt, 1 IEC. '^f^ - 14 l.L.

!`. (`_ul!t.:^CL13I. t=IIC 1,y.1 (701 1 E t 4dt6 `111CIdE

117 F !^^ l:b^hl A 1 CSC:. «L 6.3/18 1<<4''

!March 1985

Prod fqr.

National Aeronaut lde`r Spacer Administration

Marshall Space Flight Center, Alabama 35612

Smithsonian Institution

Astr0^fiy8ical:Obaervatory

Cambridge, Massachusetts 02138

tAsiro # tlil Giervetrory

g

7"

s of the

Sarva.rd-Smithoonian diker for stre h sic^a

f
2.1.3

2.1.4

2.1.5

2.1.6

2.1.7

2.1.8

2.2

2.2.1

2.2.2

2.2.3

2.2.3.1

2.2.3.2

a

13

15

20

21

24

28

39

39

40

44

45

55

S

r.

a

gG

1

CONTENTS

.RODUCTION . . . . . . . . . . . . . . . . . . . .

;HNICAL ACTIVITY DURING REPORTING PERIOD AND PROLAM STATUS . . . . . . . . . . . . . . . . . . .

Css Rotational Dynamics Analysis . . . . . . . . .

:ntroductory Remarks . . . . . . . . . . . . . . .

4atrix Transformation From Inertial Reference Frame

To Orbiting Reverence Frame . . . . . . . . . . .

Relations Between The Euler Angular Velocities And

The Body-Axis Angular Velocities . . . . . . . .

Torque And Force Package Update . . . . . . . . .

Current Status Of Rotational Dynamics Code . . . .

Implementation Of The Transformation To The Rotating Orbital Reference Frame . . . . . . . . . . .

Rotational Dynamics Code First Level Update . . .

Rotational Dynamics Computer Simulations . . . . .

Slack Tether Studies . . . . . . . . . . . . . . .

Introductory Remarks . . . . . . . . . . . . . . .

Treatment Of Damping . . . . . . . . . . . . . . .

Full Three Dimensional Treatment . . . . . . . . .

Preliminary Mathematics . . . . . . . . . . . . .

Slack3 Implementation . . . . . . . . . . . . . .

Page

4

4

4

4

5

i
Y /^

CONTENTS ( Cont.)

Page

SECTION 2.2.4 High Resolution Luss-Of - Tension Model . . . . . . 57

2.2.5 Analytical Studies Of The Slack Tether Problem 63

2.2.5.1 Introduction . . . . . . . . . . . . . . . . . . . 63

2.2.5.2 Motion Of A Medium Under Its Own Elastic Forces 66

2.2.5.3 Longitudinal Stress Waves For Simple Viscoelastic

Constitutive Relations . . . . . . . . . . . . . 68 i

- • 2.2.5.4 Analytical Studies Status . . . . . . . . . . . . 72

r

2.2.6 Concluding Remarks . . . . . . . . . . . . . . . 74

2.3 Development Of A Computer Code To Model The Electric Potential Around A Severed Tether Immersed In

A Plasma . . . . . . . . . . . . . . . . . . . . 76

C

2.3.1 Introductory Remarks 76

2.3.2 Generating The Coordinate System By Conformal

Transformation . . . . . . . . . . . . . . . . . 77 -.

2.3.3 Elliptical Coordinates . . . . . . . . . . . . . . 80

2.3.4 The Gradient, Divergence, And Laplacian In Prolate ( j

Spheroidal Coordinatas . . . . . 86 i

r ;

2.3.4.1 The Metric Coefficients . . . . . . . . . . . . . 66

t

2.3.4.2 The Gradient, Divergence And Laplacian In Generalized Coordinates . . . . . . . . . . . . . . 90

2.3.4.3 The Gradient And Laplacian In Prolate Spheroidal

t

Coordinates . . . . . . . . . . . . . . . . . . . 93

r

2.3.5 Solution Of Laplace'S Equation For The Electric j

Scalar Potential In Prolate Spheroidal Coordinates 95

e

2.3.6 Electric Scalar +'otential Around A Solid. Permeable

Prolate Sphero}.d Immersed In A Uniform Field 101

2.3.6.1 Introduction . . . . . . . . . . . . . . . . . . . 101

2.3.6.2 Matching Boundary Conditions At The Surface Of The

Spheroid . . . . . . . . . . . . . . . . . . . . 101

2.4 Potential Build-Up Due To Plasma Contactor Failure 113

ii i

i
CONTENTS (Cont.)

Page

SECTION 3.0 PROBLEMS ENCOUNTERED DURING REPORTING PERIOD . . . 114

4.0 ACTIVITY PLANNED FOR THE NEXT REPORTING PERIOD . . 114

f t

e 1

iii
The Investigation of Tethered

Satellite System Dynamics

Page 2

Contract NASS-36160

Qwarterly Report q2

For the period 1 December 1984 through 3.4 February 198S

Principal Investigator_

Dr. Enrico Lorenzini

Co-Investigators

Mr. David A. Arnold

Dr. Mario D. Grossi

Dr. Gordon E. Gullahorn

Mr. William Harrold

Prof. Robert G. Hohlfeld

March 1985

Prepared for

National Aeronautics and Space Administration

Marshall Space Flight Center, Alabama 35812

Smithsonian Institution

Astrophysical Observatory

Cambridge, Massachusetts 02138

The Smithsonian Astrophysical Observatory

is a member of the

Harvard-Smithsonian Center for Astrophysics

PRECEDING PAGE BLANIC NOT FILMED'

t
Page 3 1

Abstract

i

i

1

i

This reports deals with three major topics related to TSS dynamics:

a) The update of the SAO rotational dynamics computer code. The prole;

gram is now suitable to deal with inclined orbits. The output has been

also modified in order to show the satellite Euler angles referred to the

4

rotating orbital frame.

b) The development of the three-dimensional high resolution computer

program SLACKS. The code simulates the three-dimensional dynamics of a

tether going slack taking into account the effect produced by boom rotations. Preliminary simulations on the three-dimensional dynamics of a recoiling slack tether are shown in this report.

c) The development of a program to evaluate the electric potential

around a severed tether immersed in a plasma. The potential is computed on

a three-dimensional grid axially symmetric with respect to the tether longitudinal axis. The electric potential variations due to the plasma are

presently under investigation.
Page 4

1.0 INTRODUCTION

This is the second quarterly report submitted by SAO under contract

NAS8-36160, "The Investigation of Tethered Satellite System Dynamics," Dr.

Enrico Lorenzlni, PI, and covers the period from 15 November 1984 through

16 February 1985.

2.0 TECHNICAL ACTIVITY DURING REPORTING PERIOD AND PROGRAM STATUS

2.1 Tss Rotational Dynamics Analysis

2.1.1 Introductory Remarks -

The SKYHOOK computer code was updated with a preliminary rotational

dynamics package in 1982. This study is documented in the report "Study of

Tethered Satellite Active Attitude Control," October 1982, NASA Contract

NAS8-33691. The satellite was modeled as a rigid body and the tension

transmitted by the tether was the only external perturbation acting upon

the satellite. The three Euler rigid body equations were added to SKYHOOK

along with the nine equations relating the body-axis angular velocities to

the derivatives of the direction cosines. Once the direction cosines are

obtained, the Euler angles are computed. This method has the major drawback of integrating nine additional equations that can be reduced to three.

The body-axis angular 3elocities, wl , 4/3, W3 of the Euler rigid body

equations can be converted to Euler angular velocities and then integrated

a second time to get the Euler angles. The direct integration of the Euler

angular velocities will therefore increase the computer code efficiency and

'i
Page S

reduce the CPU time. The only drawback related to this approach is that

the equations become singular when the roll angle approach 0 0 , that is to

say (with the SKYHOOK reference frame conventions) when the vertical body

axis of the satellite is perpendicular to the equatorial plane. This condition however is a limit condition that can be encountered only in special

situations. Both integration options will be kept available in the SKYHOOK

computer code to have the greatest flexibility.

The other major modification is toe implementation of a new force

package and of a new torque package that models the effect of thruster

activation. In-plane, out-of-plane and in-line thrusters provide, in an

ideal situation, forces applied to the satellite center-of-mass (c.m.). In

the actual situation they also produce spurious torques because of

misalignments. These torques must appear on the right hand side of the

rigid body Euler equations. However, this doesn't mean that a thruster

activation produces a satellite rotation only through misalignment. Even 4

if a lateral force acting upon the satellite is directed exactly toward the

5

satellite c.m., rotations are still produced by the different moments of

inertia of the satellite and the tether. This effect is modeled in SKYHOOK

so that the tether bends differently depending upon the point of application of the external force and the mass distribution of the system.

The update of the output is also necessary in SKYHOOK. Presently the

y

computer cod. plots the satellite attitude with respect to the inertia]

reference frame. Plots cannot be interpreted immediately especially if the

orbit is not equatorial. Reference frame transformations are necessary to

convert the inertial Euler angles to the rotating reference frame Euler

ti

angles. In this way the satellite attitude is referred to the local vertical and the local horizon and the output has an immediate interpretation.

t

V e
Page 6

This way of dealing with inertial Euler angle has made it difficult to set

up the initial conditions for SKYHOOK. For this reason, the rotational

dynamics SKYHOOK has been used, in the past, for equatorial orbit only (in

this case the problem is greatly simplified). The matrix transformation

implemented to modify the output can also be used for setting up initial

conditions in order to generalize the program to any type of orbit. All

I

these modifications are in progress. Some of them will be fully operative

in our computer code in the near future.

2.1.2 Matrix Transformation From Inertial Reference Frame To Orbiting Reference Frame -

The attitude dynamic equations integrated by SKYHOOK make use of the

direction cosine matrix. From that matrix we derive the Euler angles with

respect to the inertial reference frame ( x,y,z) (see Figure 2.1.1). If we

want to refer the satellite attitude to the orbiting reference frame (xi,

yl, z t ) we must derive the rotation matrix that relates the orbiting reference frame to the inertial reference frame. If we call [Ro I ] the rotation

matrix that transform orbiting axis into inertial axis and IRBIJ the rotation matrix that transforms satellite body axis into inertial axis we have:

xIxe

YI = [Re,^ Ya

z Izo

a

t

I

I'
Page 7

and

xi x,

yi = [Rai '2r (2.1.2)

z I

I

so that

xf xn

yl_ [Roil' [Rat] Yo (2.1.3)

zI zo

where

[Roil' = [Rio] (2.1.4)

is the transposed of the matrix [Rol] while (xa, Ys, zo) is the satellite

body-axis reference frame as shown in Figure 2.1.1. The matrix [Rol ] is

derived through three sequential rotations of the inertial reference frame

(xi , yl , z I ) in order to transform it into the orbiting reference frame

(x i , y l , z ? ). The first rotation is around the axis z I and defines the

line of nodes of the orbit so that the angle P is the argument of nodes,

the second one is around the transformed axis y I l and defines the orbit

Inclination i and the third one is around the transformed axis z=l and

defines the argument of latitude y. Because of the orientation of the

orbiting reference frame adopted in SKYHOOK and depicted in Figure 2.1.1

the inertial reference frame needs to be further rotated 3.80 0 around the

=1

7

4
Page 8

r .

1 0 0 cosy siny 0 1 0 0 cosy sing 0

[RIO ] = 0 -1 0 -siny cosy 0 0 cost sini -sins cos^3 0

0 0 -1 0 0 1 0 -sini cosi 0 0 1

. ,

axis x i . This extra rotation implies a change in sign of the second and

third row of the rotation matrix so that:

and finally:

cosjcosfl- cos'ysinp ♦ sinlsini

sin-ysinflcosi sin'ycospcosi

[Rio] = sin'ycosO+

sinfisin7- -coslsini

cosrycosisinfi cosrycospcosi

L-sinisinQ sinicosp -cosi (2.1.6)

The rotation matrix between the satellite body reference frame and the

inertial reference frame was derived in the Interim Report "Study of Tethered Satellite Active Attitude Control," G. Colombo, October 1982, Contract

NASB-33691, The reverse transformation matrix, from body-axis to inertial

't.

axis is as follows:

i'
Page 9

and

cosocoso- -sinocoso- sin0sino

coa0sinosino cos0sinmcoso

[Rat] = ososinm+ -sinosino+ -sin0coso

cosOcosmsino cosOcosOcoSO

sinosinO cosOsinO cos0 (2.1.7)

In equation (2.1.7) 0 is the pitch angle, 0 is the roll and m is the yaw

angle of the satellite. All of them are referred to the inertial reference

frame so that they give the inertial attitude of the satellite.

The matrix to transform the body reference frame into the rotating

reference frame is therefore given by:

[ Roo] = [Rio] [Rot] (2.1.8)

where the matrix [RIO] and [RD.] are given by expression (2.1.6) and

(2.1.7) respectively. The matrix multiplication is not performed analytically since it is more convenient to implement it on the computer. Once

the matrix [Roo] is obtained ws can derive the Euler angles with respect to

the orbiting reference frame as follows:

al as as

[ Rao = fia #a Na

71 72 73 (2.1.9)

wx

i

x

r

r
p age 10

tan m, _ _

8y

#9

tan 01 = 71

77

tan Of 7l' +7i'/ 7^ (2.1.10)

The last relation has an unavoidable sign ambiguity. The sign is chosen

positive for convenience.

The matrix (Roo] can also be written starting from the Euler angles

relative to the orbiting reference frame (01, 01, 00). Referring to equation ( 2.1.9) the elements of the matrix can be written as follows:

'k

al = coso l coso l - cosOlsinmisinol

a2 = -siny. l cosy. l - cosOlsinolcosO1

as = sinOlsinol

#1 = coso t sindl l + cosOlcosmisinot

,6 7 = -sinol sinm l + cosOlcosolcoso'

Y3 = -sinOlcosof

-y1 = sinOtsinOl

12 = cos#)# sinOl

ryy = cos0 l (2.1.11)

Because of equation ( 2.1.8) the matrix of the body - axis direction cosines

In inertial axis is given by:

[RB11 = [R1o] T LRool (2.1.12)
Page 11

The matrix [Rio] can be built in a similar way as shown in equation

(2.1.11) by just substituting (ry, i, P) to (0', O f , 01) in the same order.

Once the matrix (RD,] is obtained, the inertial Euler angles can be derived

with expression formally similar to equations (2.1.10). This reverse process is necessary in setting up inertial initial conditions in the SKYHOOK

preprocessor DUMBEL, in the case of a generic orbit, starting from rotating

reference frame initial conditions. Before this time, the conversion had

been done in the simplified case of circular and equatorial orbit.

The sa:..,.lite angular velocities can also be derived fairly easily by

using the •, elocity composition and the matrix transformation, as follows:

W1^ fJi t)x

02 t —[Rat]= s.y

r

031 W3 lla (2.1.13) !

In equation (2.1.13) (wl l, W2 1 1 met) are the body-axis angular velocity

components as measured from an observer moving with the rotating reference

frame; (WI, w2 , W3) are the inertial angular velocity components in bodyaxis; (ll# , Ck, Cla) are the inertial angular rates of the rotating reference frame in inertial axis while the matrix [RBt] T = [Rte] transforms the

components from inertial axis to body-axis. The vector ( W1, W2, W3) has a

special significance because its components are the quantities measured by

a gyro package strapped down to the satellite. (w 2 f, w2 l , W3 1 ) are similar

components but they are obtained from the inertial vector (WI, m2 , W3)

after filtering out the orbital motion of the system.

All these matrix transformations allows the simulation program to

switch from the inertial reference frame representation to the orbiting
12

nt

of

Figure 2.1.1 Satellite Rotational Dynamics.

System Reference Frames.

r

w

G

rr

g!
^r

i

t

Page 13 :I

2.1.3 Relations Between The Euler Angular Velocities And The Body-Axis

Angular. Velocities -

In the present version of SKYROOK the body-axis angular velocities are

derived from the integration of Euler equations of the satellite. The

body-axis angular velocities are then related to the derivatives of the

direction cosines. After a second integration the direction cosine matrix

is obtained and the inertial Eisler angles are derived by means of equations

like (2.1.10). This process can be strongly simplified as follows: a)

after the first integration of the Euler equations, the body-axis angular

velocities can be analytically transformed to Euler angular velocities, b)

Inertial Euler angles can then be directly obtained with a second integration.

The body-axis angular velocities are related to Euler angular velocities as follows: E

N

W11 sinesino cosO 0

ml = sinBcos^i -sin^i 0 A

W3 cosB 0 1 ^+ (2.1.14)

c

)

4

If the transformation matrix is reversed (note that this is not an orthogonal matrix) we get:
Page 14

ilno/sinO coso/sing 0 1

cosh+ -si.no O W2

sln+ycotg0 -cosocotg0 1 ma (2.1.15)

The problem associated with equation (2.1.15) is that it becomes singular

when the inertial roll angle 0 = 0 0 that is to say when the satellite

longitudinal axis is perpendicular to the equatorial plane. This means,

for example, that in the case of an equatorial orbit, equation (2.1.15)

becomes singular when the satallite is at 90 degree with respect to the

tether (because of the reference frame conventions in SKYHOOK 0 = 0° in

this case). As the orbital inclination increases the margin on the inertial roll angle that makes equation (2.1.15) singular is reduced. Ear from

the polar orbit, however, there are no problems in running SKYHOOK with

this more efficient procedure.

The line that we are following is to have both options available in

the computer code in order to use the novel and more efficient approach

when possible. The standard option will be used for the special cases when

the satellite attitude is near a singularity.

3

` 1

I

i

i
Page 15

2.1.4 Torque And Eorce Package Update -

The external torques acting upon the satellite are as follows:

a. The yaw torque provided by the yaw control thrusters

L. The spurious torques due to misalignments of the out - of-plane,

in-plane and in-line thrusters

C. The torques due to the passive dampers of the satellite

d. The torque due to the tether tension

Since, in our computer code, the satellite attitude dynamics is modeled, by

Euler equations, all the external torques must be expressed in satellite

body-axis. The Euler equations of the satellite are given by:

Iswl - w2w3 (I3-I3)

I3(V3 - w3w1 ( I 3 - I 3 ) = NT.. + NM + NMIS + ND,m (2.1.16)

H

I 3(^3 - w1w2 ( I 1- I 3) @

In (2.1 . 15) , the vector Nza„ is the torque due to the tether tension (it is

already modeled in the computer code), the vector NyT is the torque due to

the yaw thrusters, the vector NM38 is the torque due to thruster misalignments and ND.m its the one due to the dampers on board the satellite.

I 3 , I 3 , I 3 are the satellite moments of inertia and w l , w3 , w3 are the

components of the satellite rotational velocity in body-axis. The torque

Nn has one component only, along the vertical axis of the satellite, so

that:

1

,t
Page 16 r

r 0

Nrr = 0

^N37T (2.1.17)

N 3Yr must be modeled according to the latest version of the AERITALIA control logic. In general it is a pulse modulated signal, presently with an

y

• amplitude of 1 Nm. Pulses are modulated according to the yaw angle q5 and

i

the yaw angular rate ^ in order to stabilize the yaw attitude of the satellite or to maintain a constant spin rate. )

The torque NaIa can be itemized in terms of in-plane, out - of-plane and

1

e

p

in-line thrusters as follows:y

FIT 0 0

r

NMIS 0 X 1b, + FOT X b2 + 0 X b 3 (2.1.18)

S

0 ' 0

FIL

In equation ( 2.1.18) £ IP , Fop and £ iL are the thrust level of the in-plane,

4 i

out-of-plane and in-line thruster respectively. The vectors bl, b' and b3

y

J

are the arms of the thrust vectors with respect to the satellite c.m.f

Since the satellite is spherical and the thrusters are located on the outer

r

JI

shell, the arm vector can be expressed in terms of angular misalignments as

follows: E

,

,

s Snag,, y

i

t

bk= R sinark,3 (2.1.19)

sinnk, 3

u

•rte" 4C'4017

$$jj
Page 17

where:

k = (1,2,3) and Ofk ,j = 0 if k = j (2.1.20)

ak , l0 if k j

ak , j are the cartesian components of the angular misalignment between the

actual thruster longitudinal axis and its theoretical position. In the i

case of an in-plane thruster, for example, the overall angular misalignment

is given by:

a= I = sin-' [sin 3 (a=, 3 ) + sin' (a=,3) J 1/2 (2 .1.21)

4

where xB = 1, yB = 2 and zB = 3; xB, yB and zB, being the satellite body

axis.

The thrust levels F=p, Fop and F II, must be modeled according to the ) .j

R

selected control logic. In particular the in-line thruster has a two- h

staged thrust level and it is switched on, in a continuous operation mode,

every time that the tension falls below 2N as follows:

-p

F IB = 1N if IN < T < 2N

FIC = 2N if T < IN (2.1.22)

where T is the tether tension.

The most efficient way to activate the out-of-plane thrusters is when ^'•

the satellite crosses the local vertical so that its velocity is maximum.

To further improve the efficiency it is convenient to let the out-of-plane

^i u

J
^Ly 1

Page 18

oscillation build up to 20 0 or 300 and then fire the thruster at the next

crossing. The firing time is given by:

Atop = ( m.. + 3 mt) Vop..3Foe = ( m.. * 3 mt )

(20 B... t)/Foe (2.1.23)

In equation (2.1.23) m., is the satellite mass, mt is the tether mass (that

is negligible for short tether length), Vop... is the out-of-plane velocity

of the satellite at the crossing of the local. vertical, B... is the amplitude of the out-of-plane angle of the tether before thruster activation.

All the other terms are self explanatory. It is interesting to notice that

the satellite velocity Vop... is fairly constant with the tether length.

In fact from the out-of-plane rotational energy conservation (it holds true

In first approximation because the damping is very low) we get:

Vop... = 20B.e</ [ ( m.. + 3 μe ) / (ma. + 3 Al.) ] 1/2 (2 .1.24)

In equation (2.1.24) the index "o" stands for values computed at the beginning of retrieval and A is the linear density of the tether. For tethers

few kilometers long, equation (2.1.24) is essentially independent from the

tether length 1.

The in-plane thruster control logic is theoretically similar to the

out-of-plane one. However, the in-plane control is more complicate because

the tether oscillation has a variable bias dependent upon the Coriolis

force and the air drag. The tiring time will have an expression similar to

equation (2.1.23), as follows:

AtIp = (m.. + 3 e) vtpm../F :p = (man * 7t) ( Vf3 n A... t ) /Flp (2.1.25)
Page 19

This time, however, all the angular quantities must be computed from the

in-plane bias and the in-plane thrusters must be ac'^. i vated when the satellite crosses it. The evaluation of the in-plane oscillation bias could be

performed by software processing of the radar observations (if this is at

all feasible). A simplified analytical expression for it, valid for an

exponential retrieval control law, is as follows:

1 _1 2 _ i S (R,+H 2

Ab = 3sin-1 t.ncos (Ab ) 2 m., CDPA e

] (2.1.26)

Equation ( 2.1.26) can bei further simplified, for small angle, by assuming

cos(Ab) = 1 so that it becomes very easy to compute the angular bias. In

equation (2.1.26) t. is the retrieval time constant, C D is the drag coefficient of the satelli *: e, S its cross section (the tether cross section is

neglected because the bias angle builds up when the tether is short), P A is

the atmospheric density, R. the earth radius and H the orbital altitude.

The software implementation of the torque and force package in our

computer code is presently underway. The torque PIT,,, due to the tether

tension is already available. The torque NDaop due to the satellite attitude dampers has been examined only in a simplified way. The implementation and design of roll and pitch passive dampers has been recently studied

by the Italian prime contractor. We will contact AERITALIA in the near

future to understand what the present status is in order to model the

damping torques accordingly.
w

P

i

Page 20

2.1.5 Current Status Of Rotational Dynamite Code -

Under NASA Contract NASB-32199, a version of the DUMBEL program was

created which includes the rotational dynamics of the subsatellite. This

code has a model for a damping torque proportional to the rate of change of

the angle between the wire and the vector from the center of mass of the

subsatellite to the attachment point of the wire. The program is docuy, mented in Appendix III of the Final Report "Study of the Dynamics of a

(

Tethered Satellite System (SKYHOOK)," Kalaghan et al., March 1978.

Under NASA Contract NASB-33691, the code developed in DUMBEL, but

without the damping model, was implemented in the SKYHOOK computer program.

This work is described in the Interim Report, "Study of Tethered Satellite

Active Attitude Control," O. Colombo, October 1982. The code was tested

with a rudimentary thruster model that applied a constant torque about any

of the three body axes. A special test case was also devised in which the

subsatellite processes under the torque of the wire at the same rate as the

orbital angular velocity. In this case, the subsatellite maintains a fixed

orientation between the slmmetry axis and the rotating orbital reference

frame. In order to show the orientation in the orbital frame, an equatorial orbit was run and the direction cosine matrix was rotated about the zaxis by the orbital angle in order to derive Euler angles in the rotating

reference frame. The code was not developed to the point of being able to

give Euler angles with respect to an arbitrary orbital reference frame.

N
Page 21

2.1.6 Implementation Of The Transformation To The Rotating Orbital Reference Frame -

In order to perform simulations for arbitrary orbital elements and

obtain meaningful displays of the subsatellite rotation it is necessary to

be able to set up initial conditions in a rotating orbital reference frame

and be able to display the output orientation angles in that coordinate

system. In the method of direction cosines, the coordinates of a point 'x1

with respect to the body axes are given as a function of the inertial

coordinates z by the equation

R1 = A z (2.1.27)

where

al

^1 7i

A = a3 Y3 72

a3 Y3 73 (2.1.28)

The matrix A can be expressed as the product of three matrices representing

rotations by the Euler angles. That is,

A = B(0)C(9)D(0) (2.1.29)

where

cos+ sino 0

B (0) _ -sino coso 0

O 0 1

,e -- --- - - ----
Page 22

1 0 0

• cos0 sinO

• -sin g cos0

coso sinO 0

D (d) = -sinm cosm 0

0 0 1

In the cage of an equatorial orbit with the Shuttle initially on the

x-axis, the orbital reference frame initially coincides with the inertial

reference frame so that equation (2.1.29) gives the initial conditions. In

general the transformation from the inertial frame to the rotating orbital

reference frame can be given as:

E = B(ry)C(i)D(fl) (2.1.30)

where 9 is the argument of the ascending mode, i is the orbital inclination, and y is the argument of latitude (the angle from the node to the

Shuttle in the plane of the orbit). Equation (2.1.30) has the same form as

equation (2.1.29). The transformation L gives the coordinates of a point

in the orbital reference frame as a function of the coordinates in the

inertial frame. For arbitrary orbital elements the transformation from

y

F = GE (2.1.31) 

R
Page 23

and

C = B (O 1 ) c ( 01 ) D (m # ) (2.1.32)

The .angles 0 0 , 01, and 0 0 give the relationship between the orbital frame

and the body axis frame. Equation (2.1.31) giving F as a function of six

rotations can be used to set up the initial direction cosine matrix for

arbitral.ry orbital elements and arbitrary orientation of the body axes with

respect to the orbital reference frame.

The tension in the wire is computed using the distance from the Shuttle to the attachment point of the wire. In order to set up equilibrium

initial conditions, the positions of the attachment point and the center of

mass of the subsatellite must be known in inertial coordinates. If the

coordinates of the attachment point in the body axis frame are given by the

vector F1, the coordinates in inertial space are given by

"r = F= 7 j (2.1.33)

where £T is the transpose of the matrix F. In the DUMBEL program the

vector is subtracted from the position Gf the attachment point to compute

the initial position of the center of mass of the subsatellite.

The rotational motion of the subsatellite is obtained by integrating

vs time the nine elements of the matrix F and the three components of the

angular velocity along the body axes. The angular quantities are not integrated directly. At any instant of time the nine direction cosines can be

Interpreted either as the results of three rotations as in the case of the

matrix A, or as six rotations as in the case of matrix F. A convenient

. a

t
Page 24

form of output would be to give the three angles, 0 1 , 0 0 , and O f which

define the transformation G from the orbital frame to Lody area. Since E

GE we can obtain the matrix G by post multiplying by E T to give

EEr = GEE T = G (2.1.34)

We can then derive the angles 00, 01, and O f from the matrix G. A method

of obtaining these angles is derived in Section 2.1.2 of this report, The

results of the derivation are

tano f _ -a'/j63

tan0 l _ 'y l /rye (2.1.35)

tan0 f = 112 _ + -Y2 2 /T3

where a and ry are direction cosines as defined for the matrix A. The

explicit representation of the direction cosines in terms of Euler angles

is given in equations (2.1.11) of the above cited Section.

2.1.7 Rotational Dynamics Co_j First Level Update -

As a first step in the present study, a few of the cases used to test

the program have been rerun and the output checked against the original

'test runs. The original version of the program uses explicit expressions

for the elements of the matrix in equation (2.1.28). The matrix E involves

six rotations, which would be very cumbersome to do analytically. Therefore, a set of general purpose subroutines have been written to manipulate

3 x 3 matrices numerically. Subroutine ROTMAT (IAXIS, ANGLE, A) sets up a

1,
Page 25

matrix A to perform a rotation about axis number "IAXIS" through an angle

"ANGLE." The explicit form of the matrices for each axis with ANGLE = D are

1 0 0

1 = 0 cos# sin#

0 -sins Cosa

Cosa 0 -sin#

2 = 0 1 0

sinO 0 cos0

Cosa sine 0

3 = -sin# Cosa 0

O O 1

Subroutine MATMUL (A,B,C) performs the matrix multiplication C = A x

B. Subroutine STATVEC (THETA, PUI, PSI, A) creates the direction cosine

matrix A = B(0)C(0)D(m). This has been tested against the original subroutine which uses analytic expressions for the nine direction cosines and

gives the same results. Subroutine TRNSPOS (A, AT) performs the operation

AT = AT . Subroutine VECMUL(A. V, VP) performs the operation VP = A x V

where V and VP are vectors and A is a 3x3 matrix. The original version of

subroutine ROTSTAT in program DUMBEL has been replaced by a more general

version ROTSTAT (AN,RI,AL,THETA,PHI,PSI,DTHETA,DPHI,DPSI,Y,A,DIRCOS), where

AN is the argument of the node, RI is the orbital inclination, AL is the

argument of latitude, THETA, PHI and PSI are the Euler angles with respect

to the rotating orbital coordinate frame. DTHETA, DPHI, and DPSI are the

E

j

^i(

t(C

7^

. - '. rte?'.' '^.::>>^•k;..:..
Page 26

rates of change of the angles in inertial space but referred to the axes of

the rotating orbital reference frame. The angular velocity in inertial

space is the angular velocity in the rotating frame plus the orbital angular velocity. The way the Euler angles are defined, ^f is parallel to the

orbital angular velocity, so that the angular velocity in inertial space

can be obtained by adding the orbital angular velocity n to the P which is

the rate of change of O f with respect to 'the rotating orbital frame. The

components of the angular velocity along the body axes can be obtained

using the equations

w, = ^slnO f sinO f + AcosOl

w3 = ^sinO f cosO f - Bsin0l

W3 = yCOSB' + 0

where J, ^, and 0 are the rates in inertial space.

Subroutine ROTSTAT creates the matrices

E = B(AL)C(RI)D(AN)

and

G = B(Of)C(B1)D(O1)

F r

using subroutine STATVEC and then performs the operation F = GxE. The

matrix F Is returned as DIACOS, and placed in A, the rotational part of the

state vector. ( The positional part, Y, of the state vector is not used in

this version). The position of the attachment point of the wire is comI 
Page 27

puted by the equation

R. = ETR'a

and returned via common.

• In order to display the rotation of the subsatellite in the rotating

orbital reference frame, we must have the instantaneous orbital parameters.

These can be obtained Leom the position P and velocity V. The unit vector

normal to the orbit n is in the direction 90. The orbital inclination is

cos-' n.. The argument of the node is tan-' n./n,,. If we call R a unit

vector along the node and V a unit vector in the direction n x x, we can

get the sine and cosine of the argument of latitude from the components of

P along V and z. Subroutine ELEM (P,V,AN,RI,AL) computes the argument of

the node AN, the inclination RI, and argument of latitude AL from the

position P and velocity V.

A routine called ANGROT has been written previously to rotate the

direction cosine matrix about the z-axis and compute Euler angles from the

rotated matrix. This routine has been used for equatorial orbits to give

the orientation of the subsatellite in a rotating orbital reference frame.

A new version ANGROT (AN,RI,AL,A,ANG) has been written where AN is the

argument of the node, RI is the inclination, and AL is the argument of

latitude, A is the instantaneous direction cosine matrix called E in the

general case, and PXG is a three element array containing the Euler angles

with respect to the orbital reference frame. Subroutine ANGROT performs

the operation

G = EET
Page 28

using subroutines STATVEC, TRNSPOR, and MATMUL. The vector ANG is then

obtained from G. The calls to ELEM and ANGROT are coordinated by subroutine ROTCDORD (Z,A,ANG), where Z is the positional part of the state vector, A is the matrix of direction cosines and ANG is the Euler angles in

the rotating orbital reference frame.

' 2.1.8 Rotational Dynamics Computer Simulations -

Various test cases have been run with the new software described in

Section 2.1.6. In the first case the system is in an equatorial orbit with

no initial rotation or angular velocity in inertial space. Angular damping

is included in the model. In the original version of the program, the

rotation of the wire at the orbital angular velocity causes the subsatellite to acquire the orbital angular velocity. Measured in the orbital

reference frame, the subsatellite initially rotates away from equilibrium

and then returns exponentially to equilibrium. The new software produces

the same results and agrees numerically with the previous results.

The second test case is the same as the first but with an orbital

inc'ination of 28 0 . This run appears to have integrated properly but the

output angles were disrupted by a singularity in the angular representation. When the angle 01 is zero, one cannot distinguish between O f and 01.

In this case the program arbitrary sets 01 = 0 and attributes all the

rotation to o f . In an equatorial orbit, the z coordinates are always

Identically zero in the integration. However, in an inclined orbit, numerical rc.tndoff %,'_lowed Of to assume a slightly non-zero value, causing an

erratic distribution of the angle between O f and O f after the first four

output points.
Page 29 I

In order to avoid the singularity in the angular representation, a

pair of runs was done with non-zero initial values of .1, .2, and .3 degrees for 01, 01, and O f in the rotating orbital reference frame. Runs at

00 and 280 orbital inclination gave identical time histories of the angles

measured in the orbital frame... Figures 2 . 1.2, 2.1.3, and 2.1 . 4 give plots

of the Euler angles 01, g5f, and 0 1 for the 00 inclination case, and Figures

2.1.5, 2.1.6, and 2.1.7 give the angles at 28 0 inclination. Small difference appear in the least significant digits.

The last test run is the special case where the subsatellite precesses

at the same rate as the orbital angular velocity so that the symmetry axis

maintains a fixed orientation with respect to the orbital reference frame.

The initial conditions are with ^ equal to the orbital angular velocity,

and a 280 orbital inclination. The angles O f and O f remained constant

during the run just as in the equatorial case run with the original version

of the program. Figures 2 . 1.8, 2.1.9, and 2.1.10 show plots of the angles

VS. time.

Mow
.c

No

• N M

o

N II

.n -

0 O

O 3

-I D

roc

H ro

H

0

p o

ro N 7

O II

Cro o

6

c

•N

0

U .

U

O

h I I v

C/ O

E m

H

N r

c

O

u •j

a,

.N

^ r

^ G

C O

ro U

TJ .-I

ro ro

H •N

u N

.N

O •N

N

*-I

N

O

Ha

b0

•N

u.

Page 30

N N N N N N N N N N N N N N N N N N N N N N N :1 N N N N N N N N N N N N N N !V N N N N N N N N N N N N

C C C C C C C C C O C C C C C O C C O O C C O G G C O O O C O O C C C C o O C C C C C G C G O C C G C I I/ I A/ I/ I/ I/ I A I B/ 0 1 0// I I A I I B B O A B B// I B B B 1 1 0 1 1 1 1 1 1 1 1 1 i u; of LL.0 +: :a: W W a ai :^W i li W r r'V WW'r W ti 91W ai Wi a1 W'i W 4: iu 'J ^.: '^. +: W:• N S T:D 11NT r,Tl1NF O r W L 7 T Jf II N •DT@UI@FmU1P Lh UI r•p 1 r!1 @nla O JIrv@^ Dh s ^• 3^FF T P:LNTh.DlNNG N •^G d10 a¢fP fffTO Oal rFFF@FTS @i w.lf d@h•OL G O-0C n•y

ry T C T S r 0• N T .O w OT S T N C r r O w T@ O P O r T M N •D •O M T T C T T D J O N L O w a D T J T y O r@ t 0 P.Qwr C PNGO CFNIl:NC d¢1M¢T y P MF bFJIT J; OM ONNYIG T STNC dT^^'T V T+• T'•I@h.O ThNJIT ON+1 .'l. J ^ T @TF.Dh•llaaa •O @ d +1T NP T TFFF wF r.. :vn a J P 9 T TTTTMTlF dG dr¢¢TrHld TdF¢ Jb LN a L r m r !@ Ta @Fwwi•• O.• v •CiFwL ary D ,f^ JI A elf JI 11 JI J! UI .O 9 .0 D@ D a a L JI ^ P S J1 N F 0• ¢ J1 JA D T a N JI D P a Jf L O^ O h M L (T O¢ JA O C O @ v d d d b d b d d y v d d d •I y y @@ d b d@@@@ b b@@ y@ d@@@ d O A J • @ b +^ ^r d y d Y' J1 F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F F h F F F F F F h F F F F F

C C C C G O O O C C C G O C C C O O O O C C C O O C C O O G O C C O O G C C C C O O G O O O O O G O O

O O O O O J O C O O O C O O O O O O J J O O O O O O O G O O O O O O O J O 0 0 0 0 0 0 0 0 0 0 0 0 0 9

O-0 44 @Jf OFy T O wNT@Jf LF¢O.O -•rcnoui ar ®a o wNa @JI .O ^ a a ^-•:v-rein .a^mau r w r w w r w w w w N N N ry N N h N N N T T T l^l T T T T T T@@ b!@@ d 0@@ Mf

I!.

R

Y
_ Ige 31

t

NO

• N M

oN 11

•N

A

F+ O C ^

Rf G

• N W

H

O e

F+ o

tC N

7

O'

N I I

G •

N O

6

G

•r n

O

n '-1

Q!N 11

N •O

E m

.r+

H N

F

. p

N •N

> N

yr C

G O

tC V

b ti

Itl N

f•+ •rl

v y^

.'I _ G

•t3-,r-1

N1

'-I

N

B)

7

b0 •N(L

M

N N T T N N N w r r r •'• r N r r r w N r r r r r r r r r r r r r w w r r w w r r r r r r r r r r r

C C G -COG C C C C C C C C O C--- 0 0 0 C C C C C C C O C O C C C G C C O G^ G C^

``^

/ 1 1 1 9 1/ I/ B I B B B B B B O^ O O I I I I/ O e 0 0 1 0 0 1 8 1 1 0 1// 1 1 1/ I O B

u: ^: ILii1 Y: r •: 'i 'V +: '^:W v ii i:4 V: O:•1s W ia:i i i V4i i:d +i +I i V`a:i i ai+Ik i'r r'ai r:V ii

O i d J N P T d O M1 T N d y d N 11 m T r d r O O N m L O O T 'O w d r M1 d F¢ M1 T L T T +IF m r CT r O

C O F P S p w P F T P N r O G N r M1 O F F M1 T m N 0 0 L P T m¢ F T fI P w F a N T P O d !D ¢ J: C T .•ff

^ prT Troy 7rD NF1f 7L MJ NmLaNFflf .•!SP NON y rC PNNOPTadM1POTrNN

S a N rTry ANNM•G P Yf P 6 d rli s' F frG DN d T P LT rO ON r f TC TF P FG L r Nd 0 ' D• ^. M1

f d P l I F G 1: N 0 0. a O S F a P T. o a T 7 O N r D a ro a G N O d N d O d J f P G P w 0 r 0 O b O P P m

r0 r<I 4 O V V di N 1 4 am am•C¢P amof r0

00 9 

CNS

'^¢ 

QNJIM1 N O m GNrG TO GFr a• dM1r r0 r w O rD aP m TNd ff am am O PC : ^wy a TJ1 wF O TTrFam C C? —0-0 0. L O r C J'O a¢R S OF PF Tr/1C •4O Gr y C <¢ r¢ O STP S ¢P S O m LT C GFP b G •'S tr La bP

<:: NTAM1.^. TN DFNFro a C a L N Srfl NNTrr V ® Ortl PrNN r TP PaL O LT.(l rmll rT OJ

T l v r f. P N •/ 2 w N T V 1 •O T f d F O T F O F r G T r f F O N r r fl F• f. N N r G F N r N ^ F J 1 r +. F F Y• r T F F

• e e• e• e s e a e e e s s e e• e e s s s e s e s s e e e s e e e e s

C C C C rp C C O O G O C O G O C O O O O C O G C G ^'. C O C

e

 G

•

Cs O C G O

e

 C

e Ce 

C' G

e ^s 

C C

e

G G O

• Ce 

C G

s

C• O Ce

0 1 1 0 e 0 B 0 0 1 I B 1 0 1 0 1 1 1 1 1 i 1 e 1 1 1/

^ 0 0 0> 9 0 0 0 0 0 0 O O O G O p 0 0 9 0 0 0 g D O COG O O O O O O O O C O G O O G .^ 7 0 0 0 G O ^• s• e• e s s e s e e e e e e•• e e e• s e e•• s e e e e• e•• s•• e e•• a s e•• s• •^ wNTdJILM1b P OrNTdJf LM1m P7 w:V •f d'J1 OF 9P0 rNMd fl aF DP Ow.'V Td lI LM1 D P O

w w w w w w w w w r N N N ry N ry N N N N T T T T T T T T T T O O d y V d d d d d S

y
N N N N N N N r r r r r w N r r r r N r s r w r r r r r .1 w w r r r w w w r r ^• r M w w• r w r w w r r

- C O O C O C C C C C O C O C C G O C G O O C C O C OC C C C C O C O G• C C C C O C G O O O C O C C C

1 9 1 1 1 e 0 1 0/ 0/ e/ B B// O I B O B I B/ 0 1 6 0 1/ 1/ 1 0 1 0 1/ 1/// 0 1 1// e l

:J i •• i i' "4.:i 4.: i "" "W a:V O O1-Zam V is V-1i++1a 4 •J .f1 ^J?rGl. d^ ^ 4 M ^1 r b .JP D D m •w a c JIPNN T L N'D a+l ^1Nrr y wrD +ld r TN 4f.w °I JIwrF N ^ ^ ^ Ph'•1 q ShD q /1w y1rR G q F hN q ^1 S ]^r¢ Nd r DM^r ^^f. '^ ]'11 ^ ^rN F J7 if. rONrn ••1S rNmrN T 7 .77r^ll q 77d +1 aJf JI wO f NYI 4 b y a 4 ^^ q aTy! 77 N ^ a ] ^r D

n D d s••1f TPN q %IS JN q P C¢•o LC¢ DMIf.OhO q p nnrDr ^hC ¢'f. p L ^ L xM LT^xx

C• JI Nlf axrmPrr•f ^• r0 Q q F ] Nr DMN y^ O ^ Drp C ^ D i 4nq¢n•^ 00 ^OP.OrN'nq b a.p

f1 a OP rIl:NN^IrF GdbDrn r +147 dPM a +Ir TODrr g T rn arr jf1 o e^asrm^i LC p

^ r P r T P N P N r n P Y D D G F S N M1 p G D T N¢ Jf q if C P. n¢ O p¢ F V n C r 7 b J g Q r L D Y N TrN•• S • prrP N f4 NwDP 4 •'4/•1 D r Nr'^1077 D D D Dwi(IN N.N wrDt4r NP NNw+I T T S D •4 O. •' •l if! r x lfl r T bP y' p O d p d• D ^h^ N Dm,/1 dr D^r p r j. wlw^rLhr ylF 1'N NDF Nr N.an 1!4 wlnF r ^f. .• nrn

. e s • s s e • s e s e o e s a• s e s s s o •. e e••• a s e • o • s s e•. s C C G Oe O C C O Ce C Oa

Oe C O G• O O G O

e

Ge C C Ge O G C G G C G G O C C C C O G C C C G C C C G O G O G G level epee e/e/1 e1111

000000•]eeeeuGe^e000e000000eoeaooauooc000ccoo>>0000a e s o s

^9 r,r •n q ^ll Dw 9 P JrNN0J1 .Or D PO wN^1 q NDr a T 7 wh !•1 OJ1 Dn D T O.•:aM d' ,n e•^ mP ^ r wwwwrw •.r rNNNNNNNhNNM M••1 ••1 /•1w n ww4 n4gq q db• 00 d dJl

e 32

r ,

•N M

N I I

.N

.0 -

0 O

O 3

, •O

m C

•ti ro

o

0

.

e-' o

t0 N

7

Q'O II

GIs O

C

0

u H

U

N

P I I v

O O

E o .'I

N r

C

• O

V •J

N

.ti

m r

N C

C O

tC U .ti

7 m

O O

H •1

v F7

.N

- C

^ •N

H

N

UH

a

DO .N

L:.
0

M

^ 11

• N -

k o

0 3

0

v C F

o m .N

N

r o C N

•Nti

U 11

F

• N -

O

O 6

N w

0

Id .^

C

• N II

r,-

U O

d m

N

v y

C

N O

E •.-I

• .I N

N •N

ro

• F

^ C

7 U

V co

C "i

n1 iN

•rl •M

O F O •N

F

N

- •N

m 3

N

N

N

Q

F+

7

W

.nl

LI•

I

Pag© 33

N N N N N N N N N ry N N N N N ry ry N N n N N N N ry N N N N ry N N N N N N N N N N N N N N N N ry N N N N - C C C O C C C C C C^ C C O G G O G C O O C G G C C C C C C G O C C C O G C C C C O O G C C C C C C 1 1 1 !. B 6/ I I I// 1 1// 1 1 1 0 1/ 1 1 1 1 1 1 1 1 1 8/ I B 1/ 1/ 1 1 1 1 1 1 1/ / 1 1 1 i.W 4: I^iW W W W W:i.W li 4: 'i i4 Wii Wli 'iiW4: 'iW W W W4: W W W WWfi W:r W4.W 4i{• W W i n DM SNN¢roNm2wP.MMG IIIPMMwhr ^¢ C •p C Jf D •Dh w 2 r D9rYlM hwDN •p Td O MMw

r h!• ln 0¢N ►.rr rwwhMryf 1llor MiO Tr sro r MillD ^ rr D)NMNIC rM V whlf.Mt e C PO N T OMQ+ rO NMDwTe D ^ T¢rr0)¢ndNre OMd O V ¢ yINC •p Md.p ® •! NrGmrin NhNr O r d t0 O^Gr w lCO N OOGlCr Y1N Ol'l ¢Yf¢ "l dr If. '-mr Vi "e YIC a CNNOP'•1¢NwC m'^rwN n wMMMm r p MrNJ1MSNM ^C O Nd nh DnJIT O^O•^Dd 9MMry P O+T w r p rrrrNN 00'1•d MMMMMwIMervor..¢¢M..maM,craGV^u:oa^-o:—asMe.a r.. .. l.. G r^mr.-woN¢ .n^u+v+u+v+^nu+^n ^c¢¢DeDma p ^nhrnmminra min vlm o•S In in¢emnwaorin p o.e sin ^occ < v m v d o m d e e e e d m m a e < m e m m v e a d o e v d m e e m o d d m m e In m d d e v e m e in u+ rrrrrrrrrhrhrrrrrrrrrhhrrnrnrrnhhrrrrrrrrrrrhnhrnrr

-.. r r ^ rw^ ^ r .. w

ar r

r

c G O G C C 0 0 cp 0 G C G O C'G O C G O G C G O C C G O C C O O O C O C O O O C O C O C O C C C G G O

0 0 0 c O 0 0 0 0 0 0 0 0 G O 0 0 0 0 0 0 0 G O 0 0 0 C O O 0 0 0 0 0 0 G O 0 0 0 0 0 0 0 C O O C? O

• e a• e..•.. e e e e e.. e s e. e. e s e e s• e. e s e s. e.• s .. e s •.. e•.. G rNMeNDh 9 T O r NMeN^p h y PO r NM din p r ao• G rN M V UI Oh 2 P OrNMAN.Ohm O^ ^• r r w w w w^^ r N ry N ry N N N N n N M M M M M M M M M M O m 1 e m m e V d e N
Y

Y

.L M

H

O

I

F

O

•r O

Y 3

N

F C7

'r C

N

U

C

• N O ^N

o

a0

N II

t0 -

O

F 9 .N

—0

U

CY

v.

.^ I

m n

E m

•rr

Y F

O

Y

.N

n •O

r F

F O id U

•b .1

to N

H •N v Y

F

N

O

H

a

00 .N

rL

Page 34

NNn +^ NNNwM.+w rNrr pyrNwr ----------- ... .wwr r r

C C C^ O L O G G C C C C C O C O C C C C C C C C G C G O C C C C C C C O C C C C C C C C C C C C G C

1/// 1 1 1// I I I I Y I I I I I I I I B I B 1 1 1 1 1 1 1 1// 1 1 8 1 1 1 1 1 1 1 1 1 1 1 1

:: i..14:4i itii ^ W'r'li IiLL li :al4: fi ^i i4 i i i44'4J fia Y; ILL Y4 ^r 4W 14 V 4: ili '^4:i ^4: 1LL ti 4i ti V: r;r 0@ N S n N h h V n O r G C Ill N P C. P n O. r 0 P h .p N m n ^IOC d .p w N 9 r •q N 1/! O¢ P h J L N Vf d o CO r rO+.0 P'O w¢ 3+V: d .RNi P @ pwr J r OOP+000 NrY1wC C^OPC •OC O KbN¢ •OrC r vtb R N w AmNb w ^ .. GYP li if: PN.+nr V^V j• rb rrm.fldn f.p .00 dill m nN O dV: X10 m a SOT ¢ILNrrrvOr¢C ¢ r V CN V n r Nr@ ^n rnr¢p 0 6 CO IIIr OP d C iIC .Ch N n e06r d n V dP V ^e N NT P KS.•• G OTP UTA JYP of nP ^^'tlr O •!. ^ nNn ¢N(Irr n ry nn mhh@ T!f.N ir0 rr e C Onh6l + ¢fnr 2 Ch01 w l(.r mr.Oni r C v:6C m JCL f O.p n p n ONPNnnwn rnO¢ •GPmnN@NnmCro@rvO.CPC OmmOOm •Y d J+1 Yf erCPNOrPmR¢mrC^pPO Oh O V¢.C¢ 00 LrSrnNJbOC..d Obm..4 Q' O r. O O a0¢G¢m V KOrObCnV V ¢.G •CO @GJYn Jth enrvmF Nh NN 'Jd dNmrNiV N ..I^wdR d NTroNr nTPmd X010 nu;rmJf r0 @T nNV•o NVmrNnvrnrera+rfrtn—vr.cNrinrinNN^c rlr. ..Nmrmrnrrul •-nrr

L C C C C L C C O 1:7 C C C O C C C C L G C C C C G C O C C C O C O C O C G G C C O C G C C C G C G C C

1/ 8 1 1 1/ 1 1 1 / B B B B 1/ 1/ 1 1/ 1 0 1 1 1 1

Geee000000c0000eooeeoeooeooccoeooe00000ee0000cJOOee

e

O wN n.@1/1 Or sP O r ronON Cr 2 AO r Nn Olf).O hm T e rN+l @ NOh mP a rN nd Jl •o r.,AO ... r r w w w r .. e. r ry N N N N N N N h N n n M n n M M n M n 0@ f f d 0@@ d f Jf

' -R
e 35

41

u

.N

u a

•N rn,

p

H

O II

F -

O O

n

ro O

F F

•N ro

U

U

Fo

•.+ Cl

n

.A II

N

ro O

F

•N

0

U

0)H I I

Q O Em

V r

F

• O ^ •N

au

0rh F

F o

ro U

.ti

c ti

ro ro

^, u

.N

_ F

n

N

U

fy

7

M

.N

Lt

^•. N Y NNNN r r r r r r N r r r ry r r r w r r r r r w w.^ wr r r r r r w r r w r r r w r w w r w

C C O C C 0 0 C C C C C O C C C G C G C C C C G G O C G C C G G C C C O O C t C G C C C C G C G C C

t I 1 1 1 1 1 1 1 e 1 1 1 1 0 1 1 1 1 1 e 0 0 1 1 1 1 1 1 0 1 1 1 1 0 1 1 1 o I I 1 o 0 1/ 1 1 B/ e

.. f. w W 1I 'a" n ;i 'J WW..; IV to w fi t.; 4: w w w w w w w fi 'a •O :[ fi lJ w i w w W w w w W: w w w w w IC :J •w w O .> L NNh O •JlP Y10 A wl OC P J1 m h •.IPYdm N •+1 O NO O O70 P OP Om ahNNNll di'J er+.m O

1• r d^^• Pdh +^M a T hObHd a 111 P•Lm wC wPGmNO ^• P V admen 1l ahba P Cb drax '^. 11 r J I V 1 J I O T A 0 0 0^ I h D^ l O N O •• 1 r P N m 0 0 O h h 0 r w 1 a P e O P N h .7 N M O G¢ 11 0 1(: 0 h 1; a d N P 06 ^• If O r rh S h rOh 6 P h d^mNa?.HOhP lV r bICC PY• rNr1.1FlN Nh PC MrwNd 2 D T •• wA.. rN NP NN.f1^ a•^ ••MC 1ti ^ N hhd 11 G r.^! d D'V PUt n S fl rh N-1 .O 11 •w N^1 d P•] J J1ry J1? S h 6 Ph h W 111'• •• Iw •wen r OP11100 •^OOhONV wh dNNmNNm A nHan r wdS 9N

l J 0T w;fl YN ••1h h O OwmNh A1h Om+lP+1hP dmrrTP rN Shr OJ11-N 01fl O N CwO 'O 11 •^

O J w J P N P N^ b p b a a e h s N N• e a a N¢ J^ O N C J h D P^¢ P O h G h T h ry 0 N 0 d h h P N mN r 11 A D^1 +• mh Nh^1 a J O • O a Dr1/iry N il •ewbSh NP NNr g O P ¢ a O .O •wN a a 3.'r P 9P 1( 00 d !T m ^• ry T^ a T111 Oh a 1•• h 0h •C °1wNh bHr IC N 11ry Nbh N wNbh Nw •.1hh I(1 r1°1hh

• s • •

C C C C C O O C G G O G O O O O C C C G O O O C C O C C C 0 0 0 4 C C t •: O C C C O C O C O O O C O co

01110 e/ee1 eoole Iloe1

n^oo•>>oo.^e000eo^oce^eoo^•^eoo•a000000000ae00000000000

O rN ^• d if ah 9 A 0 wNi•1011 ah m POrN+ld 11 aN m PO rro+1 d11 O^ mP0 rN+1d A Oh mP 0 w. ..r r rw ww.. wNN Nry NNNNNN •^1 ^1 +Ir ^I /•I I.11RM/.Id d 0000 OO dd11

V
c c^ o e-- c o c. o: G e C G e o d o e o o c e G c o 0

- 0 0 0 G S C C O C C O d G C O C O C d O C O e C C COO d ♦> 0 6 0 '9 O O

fc W W. 1

0• 0W 0 0i 0 A 0 . 0 0 0 0 0 0 0 0 0 0 0 

li i^'v i W'. a5 i`Wti Sr{•: :w 'a:W W W W..:W W'v'W1•'iW a....; -^amrnNa -+^nroronrin v.ort^roenNne;n -.

_' ¢ d w +` ; Yf 0¢ d r N m d V N D -^. G n r N¢ .f h C 1f1 G h N .f n D OV L ro f. d^ d why •! • ?OrS N^rO^ D n^i 11Q;ro0 20h mM Y ¢ ¢¢P /1;;hn! OCr NmC O.P Of•I+Rr?^+ T :^a^•1 d OO M O q •! O'l ./O Q •fV O']n1w M!•1 w+1w rwgq^

d;; C: d t L;; d d;;;;;; i; O .Dd;; o; OOi •O G O ,^ D O S O ^^ O•D C O O O O G O d O 6 C C^ J e O G A G O

S 3 0 J S O D C O T 0 0 0 S O D S O O C O O C 0 O^ C O DOT O. T •S9 V OOd lROV?? A'f '! V 409C •90 } •f e00^0d '1 nl /I1 w1 I•^I.1 IA Al IA AIT 1M^MI.1 N1 ^1 T •II r+! I.1 Al PI T ,w rn el•• I•^ I•I

CO

• O C* oeeoeoeeeeocoeocoeeoeecooc co CO

C 9 O O D O p O O D O O D 0 0 0 O O O O d O O O J O O O O O

e

O O J C O? O O O C O e e d O D d O 0 0 d O 0 A v 0 0 0 9 e 0 ^h w 01f^;h ¢O.O r nom• O A;.SPG rry w 0 W;hm? e

.......... N ry N ro n N N ry ry N w

Page 16

d

s o

N ti o ^

w a

v

N d

A H

N •N

O

O

C

O m0

• .i N

N

r C

C O .H .N

U h

u ^n

C U

•e•^ U

d

o F+

N a •

N ?,

N H

C L•u

F u

c o

•H u ti

N O

o v >

E s F 3 N N co

o ..

co 7

co

U r> C

e

n N

r 

O

V 0

C E

n1 V N

b .-1 A

ro H F

N •U O

V O

11 y

W L

m 4J W

00

N

v

F

.N

44
Page 37

d

,c o N N

H m o

w o

oN u

A N

H •M O O

G N

O to

•^ H i

r 

C

G O

.N .N

U 

N

U N

G O • N U

Y

o N

0o a •

N ?,

d N

N u

c o

N a>

H >

d o > N H

y

N ro

^ ^

e

N

> 

u 

bOpf

U C

Y

n 0

N N w

G d tC

N N N

co Cd F

H •^ o

0 c

d N N

OlNOtM

.Nw
_.aJ

Page 30

d

t o

u a^

O ^tl

W O

O'

H O •rl

F. •N

O

N

C N O {7

iJ

^ O

U n

J in

G O

• ri J

O

o F+

00N P N

of 1 u

N J

C O

N d

E C ?

E t 1 3

V

v Cd w 7

> J r

n

^ " f (t

4

v .^

o CdH N N • .• I O v J

3 in N

O

ti

N

O

H ^

.N

W

O C G w w w w r r- : J r r r ... r r w ^^^ G r M r

O C G G C C G- G C C^ - O C G C C C G G C G C O G C C C I P P P P 1 P P P P P O O P P P P P P 6 0 P O P P P P P P P Y 4 W i Y Y` ^ V .•^ {L • L •L 'i ^•^ •1 id 1Y J^ Vi .^ .l {V i V: 'Y Vi {V v ^ 6

oY s • n n a r O C n a-0 a N O A N e m 9 n n d r a^ a r r o O ro

irtr.s d rv .n D^i •C- a—o.a ieeoa•T M ©•r¢ N —roeG

n wa s wr .n ON O NwNllr yrTn Nr 0 +1C V OP V.D n a

G cnco: e•mnnTrry •stnmm:^nGry m++e%nrGd d.s Ta>>+^c srdwnarrvTrinrrornnryv Or ry n rr d

«• ^• nanaoGrrarrins¢r^s o^n ^on¢x^nNrdrr•o '.G s n •nr T r9n J1 '^T a` a a J1n T r •.^nNwrT 9T(Or ^' •+ ^: — nr tr f. '•^m.f. cw S wn?¢rn sN • -dJ10 i+•^OPP e o z d': d J.¢ ^e a ••• rrrp oal OfI G n ry rNJ: •a ry rrva r

d a O - r NNNN w w^p wl ..ywwry ry i•! Nry r rrry ryPwrry

s . e o e e s e 0

G C C

• 

C

o Ce G 

C

e 

C

e C. O 

Ca G C CO CO O• C C O 06 1:7 O C G• G C O G

o 

G

e . O

1 1 1 0 1 0 1 I I I B 6 B/ 1 1 1 1

D: J 0 0 0 O O O O O J O D O O O C C O 0 0 0 0 0 0 0 0 0 0 0

C J 0 9 0 0 0 0 0 0 0 0 0 0 G 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 r N r1 O n i r T P. O w ry w d if1 O n m T O r ry M y n i r x P 0

r r r r r r r r r r N ry ry ry ry ry ry ry N M .•1
Page 39

2.2 Slack Tether Studies

2.2.1 Introductory Remarks -

Several advances have been made in the consideration of the slack

tether ( tether break) problem during the reporting period.

I Perhaps most significantly, the SLACK2 program described in the last

it

quarterly report has been extended to SLACK3, a fully three - dimensional

treatment:

- The motion of the masses is three dimensional.

- The deployment boom may point in any direction, out of the orbital plane as well as in -plane.

- The tether may be deployed in any direction ( though it is still

restricted to a straight line deployment).

- The Shuttle is allowed to rotate about an arbitrary axis at a

specified rotation rate and/or rotational acceleration.

- The vibration of the boom on release takes account of the various

directions involved, e.g. the direction of initial tension due

to the tether.

SLACK2 previously treated damping in the tether ( when the masses

4f caused a section to come into tension and "bounce") in an ad hDr. fashion.

A more correct treatment was developed and applied in SLACK3.

^ RECOIL, a high resolution program to model the initial loss of tension

process following a tether break, was written. The results ( for a perfectly elastic tether) verify the initial conditions we have been using in

SLACK2 and SLACK3: the tether is recoiling at a constant veloci ty along

si

its length, with no deformation. These initial conditions had been used on

'i

i

v^

i

V
Page 40

the basis of some low resolution SKYHOOK runs and plausibility arguments.

SAO, with Prof. Robert Hohlfeld of Boston University, has also begun

consideration of slack tether dynamics, and the tether loss-of-tension

problem, from an analytical vitrwpoint. These investigations are expected

to provide valuable results not readily obtainable from computer modeling

(e.g., the typical scale size of loop formation in the cut tether), to

provide a check on computer results, and to guide the development of modeling algorithms for partially slack situations. Prof. Hohlfeld has contributed a section, 2.2.5, detailing the results so far.

2.2.2 Treatment Of Damping -

The algorithm for treating damping in SLACK2 was simple and ad h=:

When two of the masses in the model separated by the natural length of the

(assumed massless) tether segment between them, they undergo a "bounce,"

reversing their momentum components along the segment joining them. Allowance was made for damping by reducing the velocity of the rebound by a

constant percentage. This much is reasonable, simply from consideration of

the linearity of the processes (see below). What is ad bass about our

treatment is the way in which we determine what the percentage velocity

reduction is. In the absence of any better knowledge, we simply input a

percent value, which applied to all segments regardless of length.

We have now put the damping considerations on a more physical .Level,

so that we not only have the percentage velocity reductions for the various

segments to proper scale, but compute them from an input constant characteristic of the tether material.
Page 41

Although our results are certainly not new, we have not found them in

a search of readily available sources in any form easily applicable to the

situation at hand. Hence, we shall give a simple derivation from first

principles. We wish to derive the dependence of tie damping constant for a

piece of tether material on the length of the piece, under some very general assumptions. Specifically, we assume that the material is yj,acous but

not pl.Astic. That is, there will be resistance to motion, but no permanent

• change in the equilibrium state. We assume a simple Hooke ' s law form for

the elastic portion of the stress, and the viscous damping will add a term

proportional to velocity.

Consider a simple physical system: a length of tether fastened to a

wall at one and and a mass M at the other. Neglect the tether mass. Let x

be the extension of the segment past natural length. Then the equation of

motion is

Mx = — Kx — Bx (1)

where K is the Hooke's law spring constant and B is some damping constant.

Now suppose we cut the tether in its center and place a small mass m

between the pieces; denote the extension of the segment between the wall

and mass m by y. Then the equations of motion are twofold:

MR = —k(x —y)

mY = — ky — by (Z)

+k(x—y) +b(z—y)

If we let E = (x-y), a bit of manipulation arrives at the equation

MR = — ;kx — ^-.bk -- ;m (z — £) (3)
:imz = - kx - bk

force) becomes

(5)

Page 42

Let m - 0 and compare the result to equation (1). If the accelerations in the third term of (3) remain finite, which we shall simply assume

since this is not meant to be a rigorous proof, agrees it between these two

results for the motion of mass M requires that k = 2K, b = 2B. Generalizing the argument, we have that

k o, b I (4)

where £ is the (natural) length of the tether segment. For the spring

constant k, the proportionality constant is simply EA, Young's -modulus

times the cross sectional area. For the damping constant b, we shall

follow a remark in Bodley and Park (1983; p. 54) and denote the proportionality constant by C,,. It seems reasonable that the energy dissipation,

hence the resistive force, will be proportional also to the cross sectional

area of the tether; hence, we define c„ = C„/A. The two material properties E and c„ are thus input to a program such as SLACK3, which can then

compute the spring motion parameters.

We now need to determine how the velocity after a "bounce" in SLACKS

depends on the relative velocity of the pair of masses before the "bounce."

Consider only the velocity component along the tether; the orthogonal

component cannot be changed by tether tension forces. Suppose we have

equal masses m at either end of the tether (this won't be strictly true

since SLACK3 allows the segment length to vary, and the masses depend on

segment length), and a segment of length 1. Let the separation of the

masses be x+B, so that x is the stretch of the segment. The displacement

of each mass will be half x, so the equation of motion (acceleration =
Page 43

Make the usual trial substitution x(t) = e°, and write P = b/m, w2 =

2k/m. Then the exponential will be a solution if

02 t 2fla t w2 = O (6)

There are two linearly independant solutions corresponding to the roots a =

-^ ty 2 - w2 of this characteristic equation. We consider here only a

small damoino. approximation, i.e. << W. The complications become more

complex for larger damping, and the implications of extreme cases for the

model are still being investigated; the results should be good enough for

most purposes when 0 < w. The two roots then form a complex conjugate

pair and the solutions are more familiarly written as exponentially damped

sine and cosine terms. If the bounce begins at t=0, we have initial conditions x(0) = 0, *(0) = v in , say. The cosine term will drop out and the

solution becomes

x (t) = vi" a ft cos (w t)

W

(7)

where we have replaced P2-0 by W. The tether segment comes out of the

bounc3 (i.e., loses tension again) when x(t) is next zero, i.e. at t =

or/w. Taking the derivative of (7) at this time we have

V tut = Vi. enP/w (8)

The computation of segment damping parameter, and this small damping

approximation for the velocity loss on bounce, have been implemented in the

version of SLACK3 described in Section 2.2.3.2 below.

References to Section 2.2.2

Bodley, C. S., and Park, A. C., Analysts gf Tethered Satellite System gMbitai Qynamig,s fQn Selected Mission Profiles, Report TSS-83-ACP-065,

Martin Marietta Denver Aerospace, June 1983.
Page 44

2.2.3 Full Threa Dimensional Treatment -

SLACK2 restricted the motion of the masses to the orbital plane. The

same physical model has been generalized to a three dimensional program,

SLACKS. This involved the following modifications:

- The equations of motion of the masses now allow for out-of-plane

motion.

- The deployment boom may point in any direction, out of the orbital plane as well as in-plane.

- The tether may be deployed in any direction (though it is still

restricted to a straight line deployment).

- The Shuttle is allowed to rotate about an arbitrary axis at a

specified rotation rate and/or rotational acceleration.

- The vibration of the boom on release takes account of the various

directions involved, e.g. the direction of initial tension due

to the tether.

- Randomization of tether a.,gment direction is in a cone about the

nominal tether direction, rather than just in the orbital plane.

Some additional modifications have been made (or are being coded), such as

the improved damping model discussed in Section 2.2.2, and streamlining of

the input.

There are still some restrictions from a fully adequate three dimensional treatment. The two most significant seem to be

i

- The air drag (wind) is still in the orbital plane.

- The deployment boom is assumed to be attached to the Shuttle at

the Shuttle center of mass.

d `

E

r i

i!a

1

6

f
Page 45

2.2.3.1 Preliminary Mathematics -

The equations of motion of a free mass in thL orbiting coordinate

system (i.e., relative to the orbit center a the center of mass) are

x" = 3x + 2y'

y - 2x' - D (1)

z" _ -z

where the time has been scaled by the orbital angular velocity 11 to r = nt;

()' means d()/dr; D is a scaled drag parameter, assumed constant; and the

coordinate system is as shown in Figure 2.2.1. The third equation (for z)

has solution z(r) = z(0)sinr + z'(0)cost, and addition of these solutions

to the program is trivial.

Much the more complicated part, both conceptually and in the actual

programming, has been dealing with the three dimensional geometry and various angles involved. Although mathematically elementary, there is much

room for confusion and failure to coordinate conventions, so we will spell

out in some detail what we are doing.

Figure 2.2.1 shows the orientation of our reference coordinate frame,

the standard frame centered at the center of mass (Shuttle) and rotating so

as to keep the x-axis vertical. Also shown are the orientations of the

pitch, roll and yaw rotations discussed below.

Figure 2.2.2 points up a distinction (in a two dimensional case for

clarity) that it is extremely important to be clear about: that between

physically rotating a physical object or vector by prescribed angles, and

then asking for the coordinates of the new object in the original 4QXreference) coordinate system; and rotating the coordinate system in pre-
5

n

n

PITCH

e

ROLL^

^ S

n

n'

F I p '

Y'n

40000-

"C

ORBIT

EARTH

x,C Page 46

YAW

n n— ---

n,^

S

S'

Figure 2.2.1. Orientation of the reference coordinate system and the

pitch, roll and yaw rotations defined in the text.
Page 47

Y

ROTATION

Y Y'^

COORDINATE ROTATION

R

®

0000p

^

/

i

x

i

0

x^

Figure 2.2.2. Physical vs. coordinate rotation. Physically rotating a

vector or object about the coordinate rotation, and seeking the new object's coordinates in the original system, is the inverse of rotating the

coordinate system and seeking an unmoved object's coordinates in the new

system.
Page 48

scribed fashion and then seeking the coordinates in the nsx system of an

unmoving object. In both cases, the "input data" are the coordinates of

the original object in the original system. We shall call these two situations "physical rotation" and "coordinate rotation" for brevity. Brief

reflection, perhaps while pondering Figure 2.2.2, reveals that these are

inverse situations: i.e., a physical rotation by a given angle produces

the same numbers as a coordinate rotation by the opposite (negative) angle.

We shall consider the physical rotations and their generating matrices to

be the fundamental concepts.

The three fundamental units from which we shall construct all our

rotations and angles are pitch, roll and yaw. These operations are always

taken relative to a coordinate system fixed in the Shuttle (body coordinate

system). They are somewhat easier to describe if we imagine the Shuttle to

be oriented as in Figure 2.2.1, with nose pointing along-orbit and the

landing gear pointed toward the Earth. Pitch by a Dositive angle then

brings the nose down; roll causes a roll about the Shuttle's longitudinal

axis clockwise as viewed by the pilot; and yaw by a positive angle swings

the Shuttle's nose to the left. Unless otherwise specified, B will refer

to pitch, (p to roll, and 0 to yaw.

Each of these elementary maneuvers can be represented by a 3x3 matrix.

These matrices are used to transform column vectors, e.g, for pitch

s

I

i
Page 49

x cos 0 - sin 0 U x

y sin 0 cos 0 0 y (2)

z 1 0 0 1 z 0

or in more compact notation

f = po 'ro (3)

Here, the subscripts O and 1 refer to before and after the physical rotation, and the r = [x y z] T are the coordinates in a fixed reference frame,

coincident with the body frame before the pitch maneuver, of some physical

object or vector which rotates with the Shuttle. Note that r is not a

physical vector, but a column matrix of components. If the system in which

these components are measured is not obvious, we may use notation such as

ri l wm

The roll and yaw matrices are similarly:

cos 0 sin rp

Rn = 0 1 0 (4)

- sin rp 0 cos ip

1 0 0

YO = 0 cos+( -sin (5)

0 sin + cos h+

These rotations and matrices are defined so that the body coordinate system

coincides with the reference coordinate system when the Shuttle is flying

ttitude, i.e. nose forward and head up. There may be some

sign from the conventions used in, e.g., aircraft opera-
Imo. ..

_ ;r,•c- ^.;

Page 50

tions.

These matrices are examples of orthogonal matrices, i.e. those for

which the transpose is the inverse. For these particular matrices, replacing the rotation angle (e.g. D) by its negative results in the inverse

matrix. From the comment made above about the inverse nature of physical

rotation and coordinate rotation we see that the matrices for transferring

the representation of a constant vector into a system rotated with the

Shuttle will be just the transpose of those written down above; a general

result for orthogonal rotation matrices.

It now becomes relatively simple to compute the rotation matrix corresponding to successive rotations, e.g. a sequence of pitch/roll/yaw maneuvers. In each case we suppose the rotation is described in the body coordinate system. First, perform a rotation described by matrix A, then one

described by matrix B. These will transform a physical vector t=o to F, to

r ' . Initiallv, the body system is coincident with the reference xyz system; after the first rotation, it will be coincident with a system, say

x'y'z'; and finally, with x"y"z". Now, almost by definition, the coordinates of the final vector in the final system are numerically identical to

those of the original vector in the original system. I.e.,

F2 Ix• y'.• = ro Ixyz

and we may then write immediately r21xy. = A {r21x-y,.-} = A {B ralx•yr•}

= AB rolxy.• This final expression is just in the form required for she

total rotation matrix; i.e., the total rotation is described by the matrix

AB. In general, a series 21 rotations. first A` then IL then S,- will

produce a rotation described by ,tha matrix ABC.... The order is important,

R

.

i

r
(Ali'

Page 51

since matrices do not commute.

With this background, the required computations are almost obvious.

Any direction, e.g. boom or tether, will be specified by two angles, first

a pitch and then a roll which are pictured as rotating the vertical (x)

axis into the desired direction. The full rotation matrix is easily computed as PR, and the unit vector in the appropriate direction is given by

PR (1 0 O] T = [cosOcostp , sinBcosp , -sinp ] T . These unit vectors allow

us to create suitable initial conditions.

Dealing with the boom vibration and Shuttle rotation is slightly more

complicated. First, suppose we have computed the boom's position and velocity Bp and Bo relative to an unrotated Shuttle. To describe the Shuttle

rotation we have the angles O and (p of the rotation axis, which allow us

to compute a transformation matrix PR; and the rotation angle 0(t) as a

function of time. Our prescription for performing the desired rotation can

be stated: First, physically rotate the Shuttle with PR so that the body

E -axis (vertical in the body frame) lies along the axis of rotation; second, rotate the Shuttle by 0 about the body vertical axis in a yaw maneuver, with matrix Y; third, rotate the Shuttle in the inverse manner (PR]T

so that the body vertical axis is coincident again with the reference frame

vertical. From the theorem given above on combining rotations, this is

equivalent to using a single total rotation matrix

A( t) _ (PRl Yo(t) [PR]T (6)

that is, the Shuttle position at time t is given by

91 = A( t ) go (t) (7)

4
i

Page 52

The velocity is found simply as

B3 = A(t) Bo(t) + A(t) Bolt) (8)

where

A(t) = [PP] 2' (PR' T 0(t) (4)

We still need to know the unrotated boom vibration, i.e. relative to

the Shuttle. We already in SLACK2 compute the frequency and amplitude of

i

the boom vibration, as detailed in a previous report. There, to get the

amplitude we used the tether tension before the break based on vertical

deployment. If the tether pulled away from vertical by air drag, and is in

equilibrium, the tension is the same as if it had been vertically deployed

with only the gravity gradient force; if the tether was in motion or held

off vertical by t-hrustero, the force depends on factors and scenarios not

input to SLACK3. Thus, we simply use the vertical deployment tension, with

an optional "fudge factor" to adjust to specific cases, To treat the three

dimensional case we use the two unit vectors along the boom and tether, An

and 6T, which we have computed as above based on the input angles. The

,

displacement of the boom tip due to the deployed tether will be along a

direction perpendicular to the boom, in the plane defined by the boom and

I

i

tether direction. First note that the vector (613x45=) is orthogonal to tY.is

i

plane. Then a suitable vector is

XP = (6B X 6T) X 613 = 613 — ( 613 • 60) 613 (10)

which we then normalize to form a unit vector

eP = XP I I9P I (11)
Page 53

(In the exceptional case that the tether and boom are collinear, zp is zero

and cannot be normalized, but there is no displacement and we simply use 0

in place of 6p.) The displacing force, used to compute the boom deflection,

kill then be simply

Ep = {Tension} (6= • 60 . (12)

The initial, pre-break, displacement of the boom tip is then

— Ep 3(EI) 6

1? (13)

where b is the boom length and EI its stiffness. The motion of the boom

tip, prior to overall rotation of the Shuttle, is then

Bo (t) = B. t (AB) cos [2orft• ] (14)

and its velocity

BO(t) = —2af (AB) sin[2irft) (15)

where B. = bee is the boom equilibrium position and f = c l

W

EI is the

frequency of the first mode of a clamped beam; we take values specified by

NASA, EI = 1.3 x 108 lb in 8 , b = 849 in, weight W = 106.2 lb; the value of

gravity is g = 386.06 in/sec t , and cl is 0.56 cycle/sec. Rotating (14) and

Lth (7) and (8), we specify completely the motion of the boom tip.

jll

V ]
Page 54

One other application of these concepts is to a problem that at first

sight appears trivial, the generation of vectors randomized about a given

direction. We do this in randomizing the directions of the tether segments

when we generate the initial conditions. We first compute (random) angles

B with respect to the reference direction and rp azimuthally about the reference direction (typically B might be uniform with 5 0 mean and 1p uniform

on [0,27r) ). Then we create a vector with these angles about, say, the xaxis in the reference coordinate frame:

cos B

s = sin B cos (16)

sin B sin p

Finally, we physically rotate this vector using any rotation which brings

the x-axis into the desired reference direction. If the reference direction is eT , we need a i4t , ltary matrix satisfying:

1

U 0 =eT UUT=I (17)

0
Page 55

If we partition V = L u3 I u3 I u3 J, then we see immediately that u l must

be eT, and (m l , u" u3 ) must form an orthonormal set. This may be easily

done by choosing an_v arbitrary vector t not parallel to u l (we use unit

vectors along the coordinate axes for simplicity; we will never have to

try more than two) . Then form u ixt and normalize it to get u 3 i and then

take u 3 = uyxu3 . We may now perform our rotation, which will bring the

segment into an angle B with respect to the tether direction; the effect

on the azimuthal direction is not obvious, but it will leave the uniform

distribution unchanged.

2.2.3.2 Slack3 Implementation -

The above considerations have been implemented in SLACK3. The improved damping of Section 2.2.2 is also included, as well as improvements

in the randomization processes. We are in the process of coding streamlined input procedures, in which there are default values of many

parameters and a menu system allows changing only those of interest. Only

a few test cases have been run (which verify correct deployment of boom and

tether, and boom rotation); the familiar side and rear view configuration

diagrams are shown in Figure 2.2.3 for a sample run which includes boom

rotation. SLACK3, with minimal further effort, is now suitable for running
Page 56

I^

b 1 ^^vEN-6 ^^rry^p

Hfoddr^refi^nr.

no

Figure 2.:t.3. Sample results from the fully three dimensional SLACK3. A

sequence of configurations of the tether are shown in the familiar format.

Successive configurations are displayed from left to right. The "in-plane"

figure is a side view, with the orbital direction to the left. In the

"out-of-plane" figure, the system is viewed from behind.
Page 57

some complex cases in the next reporting interval to determine the effect

of deployment directions, rotation maneuvers, tether remnant length and

other quantities of interest.

2.2.4 High Resolution Loss-Of-Tension Model -

The initial conditions used by SLACK2, that the tether was recoiling

with constant velocity and no compression initially, were based on some old

SKYHOOK simulations with only a few masses, and on some plausibility arguments. It was decided to check these assumptions by means of a lumped mass

model similar to SKYHOOK, but simplified to the extreme to allow a large

number of masses to be included.

The model, embodied in a program RECOIL, is capable of handling up to

100 masses with reasonable efficiency. The restrictions made are primarily:

- The on _lv forces included are the internal tether forces. Specifically, no gravity gradient, drag, or Coriolis forces are included.

- All motion is one dimensional.

- All segments are the same length.

- Once a segment goes slack we ignore any tension that may develop

if it again becomes stretched.

The equations of motions are solved with a simple fourth order Runge-Kutta

integrator, starting from the moment the tether is severed. For simplicity, the problem is scaled so that the only free parameter is the damping

constant.

r•

i

^ 1
kage 58

VEL05 DRAG=0.0 50M

Figure 2.2.4. High resolution loss of tension: development in time. This

Is a view of a three dimensional surface whose height represents the (negative of the) recoil velocity. This is a function of the time and the

distance along the tether, in the "horizontal" plane: the tether/boom

attachment is at the right, the cut end at the left, and successive configurations recede into the distance. Fifty mass model with RECOIL; zero

damping.
Page 59

Sample results are shown in Figure "2.2.4, which shows a three dimensional plot of the velocity vs. time and distance along the tether. Note

the propagation of a sharp wave front, behind which the velocity is approximately constant. The state of the tether at the moment it goes completely

slack is shown in Figure 2.2.5. It is seen that these detailed calculations confirm the initial conditions we have been using.

Attempts to run RECOIL with non-zero damping lead to sharply oscillating results. There is probably a coding error in the program or a numerical instability. Ile shall investigate this in the next reporting period.

The efficiency of the approach used suggests the possibility of adding y

extra forces (e.g. gravity gradient) an well as damping.

It is interesting to see how this lumped mass model relates to the

tether considered as an elastic continuum. By letting the number of tether

segments become infir:ite, and scaling the ball masses and segment spring

constants appropriately, we can derive a partial differential equation,

which is (as expected) a simple wave equation. Using a standard low order

finite difference approximation we arrive at exactly the same equations as

in the ball-and-spring model. E,x eot that the cut end boundary condition

Is not treated appropriately. There are several ways to treat such free

end (zero stress) boundary conditions which maintain the order of approximation; simply considering it as another ball with the same mass as the

interior balls in the lumped mass model, however, is equivalent to treating

the boundary condition with a '*.nwer order 9.f approximation than the interior points. It it; a truism of numerical analysis that the order of the

total solution is strongly influenced by the weakest link. (This sort of

consideration is hidden by the lumped mass approach of first discretizing

. .,

i

f

t

n°4
Page 60

u

O

a

w

rn

a

K

O[

4U

aW

. \ M Q1\ ^TaM^

a^

0 .l .2 .3 4 6 6 .7 ,6 .0 I

DISTANCE FR0,9 FIXED END

Figure 2.2.5. High resolution loss of tension: Final state. Here we plot

the (negative of the ) velocity, and the (negative of the) fractional residual strain immediately after the tether goes completely slack. By the

fractional residual strain we mean the compression in the slack state divided by the original elongation. Note that strain values greater than 1

indicate that the tether segment is now extended in the opposite direction

by more than its natural length; we do not model the resultant stress.

One hundred mass model with RECOIL; zero damping.
Page 61

the physical system conceptually before writing down any equations, which

Is why lumped mass models should always be examined carefully.) By properly

treating the free end we might substantially improve the accuracy of the

model (which can be easily checked in the zero damping case). The possibility of improving accuracy by, for instance, using a finite element approach, is also appealing.

It is also possible to derive the velocity and strain upon loss-oftension by theoretl.,cal arguments. This solution is implicit in the material of Section 2.2.5, but not actually derived there. Thus, it is of

value to write down the argument explicitly. We assume there is no damping

or other energy loss. The loss of tension will then occur via the propogation of a constant amplitude tension wave with velocity c = N /E—/P . Since

the wave is propagating along a characteristic (in the space-time plane)

with no dissipation, it will remain constant in form and amplitude. Thus,

the velocity of recoil after it passes will be the same at all points. To

be concrete, suppose the tether is fixed at x = 0, is of natural length L,

and has been stretched by a small fraction 6. Then, where tho coordinate s

describes the natural length of tether from the fixed end, and x its actual

position, what we have described may be written:

x(s, 0) = (1 + 6) s (evenly stretched initially) (1)

8x (a, t) = 0 t < L - s (2)

-V t > L - s

c

where V is the recoil velocity which we shall determine below. Then we

immediately have:

x(s,t) _ (1+6)s t < L -

c s (3)

(1 ^t - .+6)-v L - s 1 t > L - s

C J c
Page 62

Upon evaluation at t = L/c, when the wave reaches the point of attachment,

this is simply

x (a, L/c) _ ^1 + 6 — Y l (4)

The final step is to compute the recoil velocity. Our assumption of no

dissipation allows us to simply equate the elastic energy stored in the

tether segment with the kinetic energy after recall. The spring constant

of the nether as a whole is k = EA/L, and the force upon extension to L+u

is just £ = -k u. Integrating Fdu from 0 to u = 6L, we get the total

elastic energy stored when the tether is stretched to length ( 1+6)L,

E•ta•eic = hk (6L) l = hEAL6 2 . Equating this with the total recoil kinetic

energy, Eki"tic = ( mass)(veloeity ) 2 = ';(pAL) V3 , gives

V = 6 = c 6 (5)

P

and then evaluation of (4) gives simply x(s,L/c) = a. I.e., there la M2

residual strain. Thus, the initial conditions we have been using in

SLACK2/ 3 are &xes5 under the assumption of zero dissipation.
Page 63

2.2.5 Analytical Studies Of The Slack Tether Problem -

2.2.5.1 Introduction -

A program of analytic studies of tether dynamics pertinent to treatment of the slack tether problem has been undertaken, While it is not

realistic to expect that a complete understanding of the slack tether problem may be obtained solely by analytic techniques, compelling reasons

motivate analytic studies as an extremely useful adjunct to numerical modeling of the slack tether. Among these reasons are:

- Analytic calculations may give insights into the physical processes important for the dynamics of the slack tether. Some of

these processes may be obscured by the complexity of the numerical simulations.

- Numerical modeling of the slack tether using the present balland- spring techniques is computationally intensive. Analytic

studies will help to direct the numerical modeling, and perhaps

to suggest new, efficient modeling techniques.

- Analytic calculations will aid in the validation of the results

of the computer-generated numerical models.

- Information generated by the analytic work will frequently complement the results of numerical simulations.

The analytic calculations will be initially focussed on processes cccuring just as the tether is going slack, i.e. going from a tensioned to an

untensioned state. In this context, it is appropriate to treat the tether

_k -s i-a.

i

I

k

i
tether suggests we apply the method of characteristics to the dynamics of a

F

E

j

Page 64

as a one-dimensional system, as long as the "loops" one would expect to

develop in the slack state do not grow to large amplitude in the transverse

directions.

The desired results of the initial theoretical effort is a calculation

of the characteristic size for loops or wrinkles in the tether after it has

gone slack. The characteristic size is of immediate importance for the

safety of Shuttle operations in a slack tether situation. We expect that

loops of small scale size (say on the order of 1 meter) will be much. less

apt to entangle the Shuttle than larger-sized loops. We shall outline

below the progress being made toward this and other objectives.

x

r

The first step in this program of calculations is to determine conditions pertaining to the development of patterns of tensioned and slack

(untensioned) regions along the length of the tether. For siml to cases,

such as a tether break near the tethered satellite, it is obvious that the

tether will go slack through the mechanism of propagating a longitudinal

(negative amplitude) stress wave along its length. More complex situations

may be envisioned as a superposition of (positive and negative) finitek!

amplitude longitudinal stress waves propagating along the tether. An immediate result (which will not be given in detail here) is that when the

tether fails, it is possible in principle, to have a nearly arbitrary distribution of tensioned and slack regions along its length. In particular,

it is possible to construct initial conditions for the subsequent dynamics

after tether failure consisting of a tensioned region bounded by two slack

regions.

Such a discussion of longitudinal stress waves propagating through the
Page 65

the tether as it is going slack and for times shortly thereafter (Morse and

Eeshback 1953, Courant and Hilbert 1962). Recall that a characteristic of

a tPave equation (hyperbolic partial differential equation) is the locus of

points in time for the propagation of an infinitesimal amplitude wavelet.

A cable or tether going slack can be modeled (for times before the

formation of well-developed loops) as a column capable of supporting a

great deal of force in tension, and of supporting a comparatively small

(but nonzero) force in compression (Landau and Lifshitz 1959, §§1'1, 20).

The formation of transverse bends in the tether which will subsequently

develop into loops is a generalization of the well-studied problem of column buckling.

Since the compressive force the tether can support is very small, a

number of effects will have to be considered which are usually neglected

for models of the tether in tension. While the tether in tension .an be

modeled as a string with a Hooke's law constitutive (stress-strain) relation, that slack tether will require a nonlinear viscoelastic or viscoplastic constitutive relation describing its resistance to transverse bending.

Also, mathematical treatments of the slack tether will require a careful

treatment of the tether inertia.

Recall that the dominant "mode" of column buckling for a column supporting a large compressive force and with negligible inertia is an approximately half-sinusoidal displacement over the entire column. Qualitatively, in this instance of generalized buckling, the tether inertia will

favor wrinkling or buckling on smaller scales and the resistance to bending

will tend to favor larger length scales. Consequently, we can hope to

identify a fastest-growing scale size for buckling which will characterize

NNo
il"7

Page 66

loop formation in the slack tether.

2.2.5.2 Motion Of A Medium Under Its Own Ela.atic Forces -

We shall develop a general equation of motion for a medium under the

influence of its own elastic. forces. A special application will be made to

longitudinal stress waves in z cable, such as the tether of the TSS, in the

next section. While the general formalism is not necessary for this first

simple application, it is valuable to introduce this notation for later

calculations pertaining to the slack tether.

Following Morse and Feshbach (1953), define the stress dyadic

6 = F,i+Fyj +F,k (1)

where 1, j, and k are the usual unit vectors in rectangular coordinates

The force across surface element dA is then S•dA.

The stress dyadic is relai;ed to the strain dyadic defined in (3) by

8 = AIGII + 2μG (2)

where I is the unit dyadic. Equation (2) is immediately apparent by transformation to the principal axes. The quantity, u, is the shear modulus of

the solid, 3 (3a+2μ) is the bulk modulus of an elastic isotropic medium,

and μ(3a+2μ)/(a+,u) is the Young's modulus. The strain dyadic may be written in terms of the displacement, s, of the solid from its equilibrium

position,
Page 67

d = (Os + sV) (3)

The equation of motion for the solid then, under the influence of its own

elastic forces is

a2s

p at, = 6•[ai6•B + pVs + μsV]

= ( a+Iy)V(V•s) + A. (Ve)

= ( a+2p)V(V•s) - pVx(Vxs) (4)

The form of the equation of motion suggests that at least part of the

vector displacement, s, may be expressed in terms of the gradient of a

potential, 0. Let

s=VO

and so

V30 = 1a2O C. 2 at=

where

C. z =

a+2μ

P

which is the usual wave equation, a 'ryperbalic partial differential equation.

M

(s)

.e
Page 68

If the shear modulus is zero, the wave speed squared is just the bulk

modulus divided by the density, as axpected. The longitudinal character of

these waves is ivmediately apparent, since if the solution of equation (4)

Is taken to be a scalar, the dyadic operator

4

Y

Is symmetric, and the rotation dyadic

$ _ - ?(Ox6)XI = 0

and so D = O, the strain dyadic. For such a solution, there is no twisting

of the medium, only stretching and squeezing, and so the wave is purely

longitudinal ( compressional). Therefore, for purposes of this discussion,

it is possible to treat longitudinal waves in a cable as a special case of

longitudinal waves in general solid media without loss of generality.

2.2.5.3 Longitudinal Strer,s Waves For Simple Viscoelastic Constitutive Relati. -

The general case of a time - varying tension in a cable or tether can be

treated using the method of characteristics as the superposition of a number o- finite-amplitude longitudinal stress waves ( Norse and Feshbach,

1953; Courant and Hilbert, 1962). For the case of a completely linear

elastic constitutive relation ( Hooke's law), wave propagation is non-dispersive. However, 'for the case of loss of tension in the tether, viscoelasticity must be taken into account.

^i
W

is

c,

Page 69

We shall examine wave propagation in the context of two models of

visco-elastic body strain, the Maxwell model and the Kelvin model (Drucker

1967). More complex models may perhaps b^ -equired ultimately to account

for details of slack tether dynamics, but toy the present we shall use

these models to illustrate changes in elastic wave propagation by small

perturbations of viscosity.

In the Maxwell model, the strain is taken to be the sum of elastic and

viscous components, i.e.,

E = CO + ev

and the elastic and viscous stresses are constant at all times,

ov=o°=a

The elastic stress and strain are related by

0 = a/Ep

and for the viscous component,

dEv

dt — a/CM

independent of all prior history. The equivalent mechanical system for the

Maxwell model is a spring and dashpot connected in series.

;T .
r4r

^i

Page 70

We can write a differential equation for the total strain, e,

dede— •de° 1 do a

+ _ — +

at

dt dt EM dt CM

(Note that treatment of the strain, e, in this section is equivalent to the

treatment of the displacement, s, of the previous section.) Since there is

no restoring force associated with the displacement of the "dashpot," i.e.

the viscous strain, the perturbed wave equation for the Maxwell model is

XZ E° =

>>

F.

2 (E° + Ems)

but de"/dt = a/CM, and so the e4uatiort can be rewritten as a wave equation

for the stress in the cable,

a2a 82 EM8a

8x z - c,a at z " = c,=cM at

We shall consider the viscous term as a small perturbation on the

hyperbolic partial differential equation and attempt to find a solution of

the form

o= f (xtc^t) exp [-al x-x, 1]

a wave spatially damped from some initial point of propagation, x,. Substitution of this form into the wave equation gives, to lowest order in

viscous terms,

_ EM

a

2C,CM
Page 71

which relates the model parameters to the damping parameter of the wave.

Therefore, treating a damped wave in the context of the Maxwell poses no

difficulties or inconsistencies.

A similar treatment may be undertaken for the Kelvin viscoelastic

model. In the Kelvin model, the stress is at all times the sum of the

stresses in an elastic and in a vit;cons element,

o=o•+ov

and the strains are equal,

CO = E v = E

with

de

O* = EkE and (,`kav _

dt I

C

independent of the history of the system. The equivalent mechanical system

p

of the Kelvin system is a spring and a dashpot in parallel. The wave j

equation for the Kelvin system may, after some minor manipulation, be cast

in the form of a wave equation for the strain,

YY

i

a 2 1 a2 CK d f 8 21 a 21

(ax e

E_

C^ 2 at 2 )

E

E ax 2 at ( 2 C ate) J

I

which manifestly describes the propagation of a damping wave with damping .

coefficient,

COCK EK

P

r

r
Page 72

Generalization to more ...mplex viscoelastic relations is straightforward. It can be seen on the basis of these calculations that we can traat

the propagation of tensioned and slack regions in a viscoelastic tether in

a tractable fashion as a perturbative extension of the method of characteristics.

2.2.5.4 Analytical Studies Status -

We have begun a series of calculations relevant to the slack tether

problem when tension has been lost. At present, this work has been concerned with the situation in which the tether is going from a tensioned to

a slack state. Further research will address dynamical problems of the

slack tether.

This program of calculations has begun from the application of the

method of characteristics to the dynamics of longitudinal stress waves for

a tether with a purely elastic constitutive relation. This work has first

of all verified analytically the approximate validity of tl,e constant velocity initial condition used in simulations of the slack tether in the

program SLACK2. Further, it has been demonstrated that for a tethe r with

an elastic constitutive relation, arbitrary patterns of slack and tensioned

tether segments can be built up by superposition of longitudinal stress

waves. In particular, the important result has been obtained that is possible to have a bounded, tensioned tether segment bounded by two slack

segments, the boundaries moving with the velocities of characteristics for

the hyperbolic equation for longitudinal stress waves.

i
^7.^... ••.4*F

l

Page 73

Work is underway developing the dynamics of a slackening tether exhib--

sting a viscoelastic constitutive relation, e.g, a Maxwell or Kelvin relation. As the tension of the tether diminishes, eventually viscous effects

must dominate the dynamics of the tether. This work has shown, as a perturbation upon the hyperbolic equation for the elastic constitutive relation, that the stress waves do not exhibit a nonlinear steepening for physically reasonable constitutive relations.

The work outlined above is of importance for defining the initial

conditions for a semi-analytic treatment of the slack tether problem. As

outlined by Landau and Lifshitz (1959), a slack tether can be modeled as a

rod with small (not necessarily positive) tension with inertial forces and

resistance to transverse bending. A taut string is the opposite limit in

which a positive tension dominates the dynamics. The advantages of this

type of mathematical treatment are several: 1

1) A relatively consistent formalism may be developed applicable to

both tensioned and slack tether segments.

2) The beginnings of loop formation may be treated in a manner similar

to the problem of a column buckling under compression. This would allow

identification of the scale size of loop formation in the slack tether.

3) No nonphysical mathematical singularities are introduced into the

behavior of the tether at zero tension.

4) Provided deviations from a linearly extended tether remain small,

tether dynamics will be accessible by relatively standard mathematical

techniques.
Page 74

References to Section 2.2.5

Courant, R., and Hilbert, D., Methods of Mathematical Physics, Interscience

Publishers, John Wiley & Sons, New York (1962).

Drucker, D.C., Introduction to the Mechanics of Deformable Solids,

McGraw-Hill Book Company, New York (1967).

Landau, L. and Lifshitz, D., Theory of Elasticity, vol. 7 of Course of

'Pheoretical Physics, Addison-Weslay Publishing, Inc., Reading, MA

(2959) .

Morse, P.M. and Feshbaek, H., Methods of Theoretical Physics, McGraw-Hill

Book Co., New York (1953).

2.2.6 Concluding Remarks -

During the reportinc,;, ?.nd SAO has substantially enhanced the simulation program SLACK2/3: The troatm.,yr.t of viscoelastic damping has peen

improved, and the program has been extended to full three dimensional capability. In the next reporting period we expect to finalize the coding of

SLACKS (simplified input procedure and a few other minor items) and use the

program to conduct a series of operational mission simulations relevant to

safety issues. We shall also investigate the scale length for critical

damping in the tether 'material and the implications for the bounce calculations in the SLACK3 model.

An efficient, high resolution (yet strongly restricted) program for a

partially slack ball-and-spring model has been developed. The program has

confirmed the initial conditions being used by SLACK2/3 (which assumes the

tether has become completely slack and hence cannot follow the initial

loss-of-tension). These initial conditions have also been confirmed analytically in the undamped case. In the next Period we hope to extend the

program to include damping. Other possible extensions include modeling
5r'r

Page 75

with the gravity gradient force, and 2 or 3 dimensional modeling with coriolis and drag forces; and proper treatment of the free end boundary condition from numerical analysis of the governing partial differential aquation. Perturbation (stability) studies might also be feasible.

Analytical studies of the slack tether have been initiated, with the

intent both of supporting the SLACK2/3 simulations and producing significant results not amenable to direct modeling. In the next period, we

expect these studies to provide a good idea of the scale of the wrinkles or

loops resulting from the initial buckling as a cut tether goes slack. This

result will be relevant to safety studies and will also suggest how initial

conditions for SLACK2/3 should properly depart from purely symmetric recoil, a consideration now handled by ash hac randomization.

A continuing frustration, one which is becoming critical for some of

the projected effort, is the lack of hard data on tether material properties beyond the simple axial stiffness (EA). We have asked Martin Marietta

Corporation to send us what data they have on tether properties, and we

also intend to request samples of tether from MC, upon which we may perform some s_ ple experiments.
Page 76

2.3 Development Of A Computer Code To Model The Electric Potential Around

A Severed Tether Immersed In A Plasma

2.3.1 Introductory Remarks -

In the Quarterly Report 01, we showed that the prolate spheroidal

model of the broken tether was very inaccurate at distances close to the

wire. We are now developing a numerical model which represents the wire as

a small cylinder. The program uses the Poisson's equation to calculate the

electric scalar potential on an axially symmetric three-dimensional grid.

The grid lines are closely spaced near the end of the wire but the spacing

increases geometrically as one moves away from the end. The lines parallel

to the wire are spaced 0.1 mm apart inside the wire and within one wire

radius. Then the spacing is increased geometrically as the distance from

the wire is increased. For example by choosing a geometric factor of 1.2

(each grid spacing is 207, larger than the preceeding one) and by adopting

for the first grid spacing a value of 0.1 mm, then 100 grid lines will span

a distance of 41.4 kilometers. The radial grid lines are, likewise,

closely spaced near the end of the wire, and their spacing increases as the

distance from the and increases.

While this approach could, in principle, be used to model the entire

tether, we decided to develop and use a more efficient model. Rather than

numerically modeling the entire length of the tether, the voltage at the

grid points within a region spanning a space of 10 to 100 meters from the

end of the wire is initialized by a prolate spheroidal model for the wire

immersed in an ot:ierwi:e uniform electric field created by the motion of

the wire through a uniform transverse magnetic field (VxB). The plasma

density is initialized at this potential and then the trajectories of a
Page 77

number of "computer particles" each of which represents many real particles

is ca:culated. At each time interval the electric charge density on the

grid is calculated and a discrete implementation of Poisson's equation is

used to numerically estimate an updated electron potential. In this manner, the motions of the electrons and Sons near the end of the wire can be

modeled. While we have made considerable progress in developing this computer program the code is not yet operational.

2.3.2 Generating The Coordinate System By Conformal Transformation -

Prolate Spheroidal Coordinates are the most useful coordinates for

calculating the electric field around a long wire. Moon and Spencer give a

most convenient: transformation for generating an orthogonal set of hyperbolas and ellipses;
Page 78

z = x+iy (2)

W = q+10 (3)

Given a planar grid in the w = 9+10 coordinate system this transformation

maps the straight lines q = constant and 0 = constant onto ellipses and

hyperbolas in the z = x+iy coordinate system as shown in Figure 2.3.1.

Rotation about the major axis of the ellipses gener&tes the prolate spheroidal coordinate system. The pur pose of this Section is to derive the

relationships required to perform calculations in this interesting coordinate system.

i

M

1

1

1

1

t
Figure 2.3.1 Prolate spheroidal coordinate system (a twodimensional view) obtained by conformal transformation.
Page 80

2.3.3 Elliptical Coordinates -

Substituting (2) and (3) into (1) we find that

x+iy c coshq + 10)

= c (en• 10 + e-n-10)

2

2 eq (cos0 + isin0) + e' q (cos0 - isin0)

= c

2 (e

q + e-n) cos0 + i (en - e''+) sin0

from which we see that

x = c coshq cos0

y = c •inhq sin0

Therefore,

X

3

+

y 1 = eos2 0 + sin 2 0 = 1 (7) c coshq c sinhq

Thus we have the equation for an ellipse where the major semi-axis is given

by

(4)

a = c coshq (8)

and the minor semi-axis is given by

1 '

;1
b = c sinhq

From the definition of sinhq Lind cosh" it follows that

i

M' cosh°q - sinh°q = 1

and so

i

Page 81

(9)

(1Q)

a° - b' = c 3

(11)

At this point it is useful to introduce two vectors which extend from the

two foci at x = fe to any point on the x-y plane; let their magnitudes be

ri =(x-c) 2 +y2 (12)

and

r7 = \1(x+c) 2 + y2 (13)

Recalling that

x = c coshq cos8 (5)

y = c sinh" sin8 (6)

It follows that
Page 82

r2 = (c coshn cosO -c) 2 + c z sinh 2 n sinaB

c%/cosh e s cos R O - 2 ci.shq cosO + 1 + sinh a q sinaB

cosh2 q cos z O - 2 coshq cosO + 1 + sinh a q (1-cos2B)

C VCOS20 (cosh 2 q - sinh a q) - 2 coshn cosh + coshzn

r; = c -V/cos 2 B - 2 coshq cosO + cosh2n

r 2 = c f cos8 - coshn) 2

rl = cicoso - coshn)

But coshn 7 1 1 cosO, therefore

r2 = c (coshn - cosO) (14)

Comparing the definition of r l and r2 one sees that the same logic used to

evaluate rl in terms of q and 8 results in

r 2 = o(coshn + cosB) (15)

Therefore

ri + r 2 = 2c coshn (15)

►1
r ' - r1 = 2c cosB

Page 83

(17)

Two new parameters are now introduced which simplify the equations

coshrl = ri + 12

2c

and

μ coaB

r2 - ri

= =

2c

It was shown earlier° that constant F = coshq defines an ellipse, Now we k

shall see that constant p = cosB defines a hyperbola

X y l^ = cote:, ; - sinh2q = 1 (20)

's

[c cosB, - [c sinBJ

G

r

The family of hyperbolas is orthogonal to the family of ellipses. This

results from the conformal transformation whi,.h translates squares in the

(q,B) coordinate system into curvilinear squares in the x-y coordinate system. This is easily verified by testing and mapping function to determine

if the Cauchy -Riemann conditions

i

I

1

2x __ 2y 3x _ a

8q 8B' 8B q

(21)

s

i

are met. Again recalling that f

1 i ^

I

i I X = c coshq cosB i

II i

I

i

i

-.FIT ^ M:. _..

(18)

(19)
Paye 84

y = c sinhq sing

it follows that

aq = c sinhq cosd

a9

= -c cosh* sin g (22)

ae 

= c Binh* coeB

aq 

= c cosh* sin g (23)

Thus we see that the Cauchy- Riemann conditions are met and the set of 1

hyperbolas are orthogonal to the set of ellipses as illustrated in Figure ')

u

2.3.1. Rotating this coordinate system about the major axis of the ellipses generates the prolate spheroidal coordinate system ) _ n 3 dimensions

μ, and t/i. It is conventional to rotate the ellipse with its major axis i

coincident with the z axis of the cartesian coordinate system. To accom-:

plish this one must replace x by z and y by R = x 2 + y= in the foregoing

equations. The equations I^

s

x = c cosh* toad (24)

+

y = c sinhq sind (25)

are replaced by

i

z = c cosh* cosd ( 26),

R = c sinhq sind (27)
Page 85

x = R coso (28)

y = R sinO (29)

Therefore the relationship between prolate spheroidal coordinates and cartesian coordinates is given by

x = c sinhq sin g coso (30)

y = c sinhq sing sinO (31)

z = c coshq cos8 (32)

or

x = c ` 2 --1 1-μ2 cosO (33)

y = c ` 2 --1 1-μ2 sinO (34)

Z = c ^ p (35)

where

F = q = cosh rl r2

2c

μ = cosB = r a - rl

2c

' t)

p

4

F rtl

j

1

1
Page 86

r, = VX 2 + yz + (z-7) 2

r2 Vx2 + y2 + (z+c) z

2.3.4 The Gradient, Divergence, And Laplacian In Prolate Spheroidal Coordinates -

2.3.4.1 The Metric Coefficients -

To perform field calculations in prolate spheroidal coordinates it is necessary to calculate the gradient. divergence and curl in that coordinate

systeu.. The easiest approach to this problem is to derive the gradient,

divergence and curl in generalized coordinates and then to evaluate the

metric coefficients h i , hz, and by using the transformation from prolate

spheroidal to cartesian coordinates derived in Section 2.3.3.

To calculate vector relationships in generalized coordinates one must

include the metric coefficients. In cartesian coordinates differential

length ds is given by

ds = V ( dx ) 2 + (dy) 2 + (dz) 2 (1)

In cylindrical coordinates ds is given by

ds = (dr) 2 + (rdO) z + (dz)

not h,-

I

da = (dr)- + (d B) 2 + (dz) 2 !

I

I

i
Page 8^

In any orthogonal system,

(ds) 8 = (hldgi) 3 + (h3dg3) 3 + (h3dg3) 3 (2)

where hl , h3 and h3 are the metric coefficients.

Now to determine the metric coefficients we consider the following total differentials:

dx °^ 41I dg1 +C3q _) ) 2 +(ag3)dg3

(3)

dy °ra `dql +/^1dg 2 +1^ldq^ (4)

dz

 =^a

Z d l + ZI ) d 2 + az ,d (5)

q l) q C q2 q q3. q3
Page 88

where q l , q2 , and q 3 are '^I:e generalized coordinates.

Squaring these terms results in:

(dx) 2 = [^—1 d1g1J2 + 1 42 dq z^ 2 + ^aq3 dq 2 +cross terms (6)

(dy)2 [a dql^ 2 + ra dg 2^ 2 + la dg 31 2 + cross terms (7)

ql q2 l q3

(dz) 2 rza dq l^ 2 

+ 1a42 dg

2]

2 + l

a 

z dq^ 2 + cross terms (8)

In any orthogonal coordinate system the Pythagorean theorem demands that the

cross terms must vanish. Thus, assuming that our generalized coordinate

system is orthogonal, it follows that

(ds) 2 = (( ax 2 + (x/ 2 + ( az 1 2 (d )2

8q1aql agl) ql

+ l`ag2I 2 + C4 ) 2

2 + \a42) 2 

I(dg2)2

F f(ag3) 2 + la ) 2 + \aQg) 

2 (d4 3 ) 2

(9)

I
Page 89

Comparing equation (9) to (2) we see that:

hi \a41) 2 + (6^) 2 +(aq ) 2

i i = 1, 2, 3 (10)

Now to . evaluate the metric coefficients for prolate spheroidal coordinates we

make use of equation (10) and the coordinate defining equations obtained by

the conformal transformation defined in Section 2.3.3.

x = c sinhn sin6 cosy o <_ n < 00 (11)

y = c sinhn sin6 sinr o S 6 < n (12)

z = c coshn cosh o < ^ < 2rr (13)

where

ql=n, q 2 =e . q3=^-

The metric coefficients are calculated as follows:

hl

( 2x l 2 2 2

°/ + lan/ +`aei

(c coshn sin6 cos^) 2 + (c coshn sin6 sinr) 2 + (c sinhn cose)2

= c2 cosh 2n sin 26 + c2 sinh2n cos26

= c 2 (1 + sinh 2n) sin 26 + c 2 sinh2n Cos 26

hlc(sin 3 26 + sinh2n) (14)

h2 = (ae! 

2 

+ (46 ) 2 + f NO 2

= (c sinhn cosh cosO 2 + (c sinhn cosh sin$) 2 + (c coshn sin6)2

= c2 sinh2n cos 26 + c2 cosh 2n singe
Page 90

h2 = hl (15)

h 3 = 1 a ^ 2 + ^ / 2 + ra h / 2

= (c sinhn sing sin^) 2 + (c_ sinhn sine cos^)2

h3 = c sinhn sine (16)

2.3.4.2 The Gradient, Divergence and Laplacian in Generalized Coordinates

The gradient in generalized coordinates is given by:

V^ = ul aL + u2 a^ + u3 9^

(17) hl 

aql h2 dq 2 h3 aq3

where u l , u2 , u3 are the unit vectors of the generalized coordinate system.

The divergence of a vector is defined as the limit of the total flux emerging

normally from a closed surface divided by the enclosed volume as the volume

approaches zero. Thus

V•A = lim 1 A•dS

OV+O AV Av

where

A Alu"1 + A 22 + A 3 u 3

Consider a differential curvilinear cube centered at (q l , q2 , q3 ) as shown in

Figure 2. 3.2. The area of a differential square in the q 2 , q3 plane at the

center of the cube is h 2h3Ag2Ag 3 . The area of the end of the cube located at

q l + Aq l is there?ore given by Taylor's theorem as

_Z

(18)

A

i
Page 91

+ "q ' q + "9z q + Aq3

2 2 2 3 2

A u

1 1

Lq 2

Aq3

n 2 - 2 , Q3 3

^2 , Q 3 + 23 ]

F

"r.

G

U

3

Figure 2.3.3 Curvilinear Cube Illustrating

Divergence in Curvilinear Coordinates

.;, _-
V1

_ ^mow-*^---^

Page 92

Aq

As = Ag2Ag3 Ch1h2 + ql (h2h3) 21

It follows that the normal flux through this surface is given by

` A(D - Al°AS

Aq

o pg2pg3 Ch1h2A1 + qa l (h2h3A1) 2 ]1

Adding to this the flux through the opposite side of the cube in the -ul

direction, we see that the total flux exiting normally in the,u1 direction

is given by

Aq l

A@1 

1 

2h3 + onl (h 2h3A1) 21 J 

Ag2Ag3

(( Aq

- [h 2h3 + aq (h2 h3A1) \ 2 1 j pg2pg3

1

But the volume of the curvilinear cube is

Av h1h2h3Ag1Ag2Aq3

a

Therefore the v i contribution to the divergence is given by

(v°A)1

AV 0 6Vragl (h2h3a1 ) Ag1pg2Ag3

1

Applying this procedure to the uther two axes and summing the results, we see

that the divergence of A in generalized coordinates is given by

v°A h1h2h3 C aq1 (h2h3A1) + aq2 (h1 h3A2 ) + aq3 (h1h2A3 )^ (21)

Now that we have the equation for the divergence and the gradient in generalized

coordinates, it is clear that the Laplacian is given by '

V2^ 0 1 C 8 h2h3 ^1 + 8 IL)

+ ^

hlh2 ^1^(22)

hlh2h3 Lagl hl aql aq1 \ h2 a

hlh3 q

.1 aq3 h g aq3

1.

(19)

(20)

-V i
(24)

(25)

A^

Page 93

i 2.3.4.3 The Gradient And Laplacian In Prolate Spheroidal Coordinates -

Substituting ( 14), (15) and (16) into ( 17) gives the gradient operator

in prolate spheroidal coordinates

fi a^ + e aQ + 2¢

(23)

0O o c^e + sinh z n c sinz e + sinh zrl N sinhh sine 3^

To calculate the coefficients of the Laplacian we utilize the following facts

h 

I 

h 2 h 3 = c 3 (sin 26 + sinh 29) sinhh sine

h3h

= h 3 = c sinhh sine

1

hhh 3 _ h

3 = c sinhh sine

2

h 1 h 2 _ c (sin 20 + sinh2n)

h3sinhn sine

(26)

r

w

l

E

(27)

Introducing the prolate spheroidal metric coefficients into the Laplacian in

generalized coordinates results in a cumbersome equation. It is clearer to

deal with each term of (22) individually:

as \ hhl3 3 )

= 3 (c sinhh sine an )

 an

z

= c sine (sinhh an1 + coshn an }

(hlh3

a3 h

 4

q, a0 c sinhh singa3

q2 z 

= c sinhh tsine a: + cose Be

(29)

-- ---

_

(28)
Page 94

2 3 ( h 132 L1 a 2 r (sin2 6 + sinh 2g) 2^ 1

2q ` h 3q / a^ l c sing sine ay 1

m c (sin26 + sinhzn) az

sinhn sine 8^ (30)

Substituting (24) through (30) into Laplace's equation (22) results in:

z

Qz e 1 sine (sinhn a + coshg 2¢

c z (sinz e + sinh z n) sinhn sine an Bn

(

a26

^.} (sinze + sinhzn) az^ 1

+ sinhn sing 5-6t+ core 86 + sinhg sine 9WT J

which simplifies to:

Oz^ - c sin e l+ sig

h TOUgz + cothg na + e3 + xc 6 ae

1 3Z¢

+ c sine adz (31)

i

II
nm^+B cos m* (5)

n negative integer.

C —

J,I

Page 95

2.3.5 Solution Of Laplace'S Equation For The Electri c Scalar Potential In

Prolate Spheroidal Coordinates -

It was shown in Section 2.3.4 that Laplace's equation in prolate

spheroidal coordinates is given by

Oz^ = c z ( sin z 6 + sinhzn) ^24 + cothn 8n + 36^ + cote a6

1 azb

+ c—z sinh z n sine 3Tr (1>,

We solve this equation by the technique of separation of variables. Assume

a solution of the form

^ = H(n)e ( 6) `P 3 (2)

Substituting ( 2) into ( 1) and multiplying by c z sinh zn sin z e / Hey results in

sinh z n sin 2 01 d 2 dH 1 d z0 d0 1 dz4'

sin 2 6+sinh z n H dr^j + cothn do + d6 +cote de t d7 = 0 (3)

The last term is a function only of 1P which does not appear in the rest of the

z

equation. In order for the equation to hold for any value of ^, dz must be

constant. By convention we choose this constant as -m z , therefore

z

d^z + M'T = 0

From the geometry of the problem we know that ( 4) must have a single valued

trigonometric solution; therefore, we deduce that

(4)
Page 96

Equation ( 3) may now be rewritten as

1 d z H

+ cothn dH + 1 0 + cote d01 _ msin + sinhzn) z ( z 6

0 (6) H dr^j do 0 (

dz

wde d9 / sinh n sin B

Now the rightmost term is easily shown to be

m2 l sinhl z n s n 9 + 1 J (^)

Therefore equation ( 6) is easily rewritten as

1 d z H dH mz 1 d z 0 d01_ m z

_

H (dn2 + cothn do - sinhz n + 0- (dBz + cote M sin z e - 0 (8)

Since the first and second terms are dependent only on n and the third and

fourth terms only on 6 we may reasonably assume that they are independent of

each other. Therefore each is equal to a constant which by convention is

chosen to be k(k+l)

z

l

z

H dn2 + cothn drxjl - sinh z n - R(k+l)

0 Cd8 

+cote 

d6) sin z e = - k(k+1)

These two equations may be rewritten as

dnH + cothn do - [k (k+1) + 

si^j H = 0

d 20 del + cotB de + Ck(k+1) - sin'e J 0 0

(9)

(10)

(11)

(12)
Page 97

Equations (11) and (12) may be simplified by the introduction of two new

variables

^ = coshn

and

P = cose

Considering first equation (11) we see that

(13)

(14)

.I

do = sinhn

d 2

dnz coshn

sinh z n = cosh 2n-1 = V-1

dH _ dH dC =dH

do d'^ do sinhn dE (15)

z

d = do (sinhn) d^ + sinhn did d^

= coshn dH + sinhn dE d dH d^ do do do

z

_ dE + sinhzn dam=

z z

ibstituting equations (15) and (16)into (11) results in

«z-1) do + 2E d4 - I R(R+1) + i-1 J

H - e

zis U Legendre's equation with which oe will deal later.

(16)

(17)
EF +J

Page 98

Applying the same method to equation (12) we see that

ARd0 = - sine

d2u

- cosB = μ

6322

sin g e = 1 - cos 2 0 = 1 - 112

d0d0 du = _ sins d0

de, = d{i d0 dμ

d 20 Cosa d0 - sine d d0 d0 2 = dμ d8 dip

- cosy dO + sins • du d d0

du d8 dp d i

_ - Cosa do + sin 28 du

(1-μ2) du0 - μ a0

Substituting equations (18) and (19) into equation (12) results in

(1-μ2) dd^0 - 2μ do + [R(2+l) - m1

μ21 

0 = 0

Multiplying (20) ' by -1 results in

(μ2-1) a + 2μ do - [R(R+1) + μms, 0 = 0

Thus the separated equations for H(TI) and 0(0) are of identical form and are

both Legendre's equations for which the solutions are associated Legendre

functions of the first kind Pm(&) and the second kind QQ(O.

(18)

(19)

(20)

(21)

- — 11
Page 99

Thus the separated equations for H(q) and 0(0) are of idential form and are

both Legendre's equations for which the solutions are associated Legendre

functions of the first kind Pe (£) and the second kind Qf (£) . To solve

equation (21) we ser m=o and solve for the zero order Legendre function.

The first few of whi=h are listed here:

Pe (£) = 1 (22)

Pi(£) _ £ (23)

P3(C) = 1 (3£ z -1) (24) ^,..

2

(

^ `f

Qo(£) = 2 log ^£+1J (25)

I

I

i

1

f+

i
V

Page 100

The Legendre functions of order m greater than zero are the Legendre associated functions which are derived from the zero order function by the

following relationships

Pt (e) = (£ Z -1)°^^ d°a££) (28)

and

Qt (t) = (^Z-1)°^a d ae

£( ) (29)

These functions are very useful for calculating the electric and magnetic

fields in and around many objects which are most conveniently described in

spheroidal coordinates.

!/ r

l
(3)

r1+r2

= 2c

1....... , , o

Page 101

2.3.6 Electric Scalar Potential Around A Solid Permeable Prolate Spheroid

Immersed In A Uniform Field -

2.3.6.1 Introduction -

2,3.6.2 Matching Boundary Conditions At The Surface Of The Spheroid -

It is shown in Section 2,3.5 that when Laplace°s equation is solved in prolate spheroidal coordinates by the method of separation of variables the

resulting ordinary differential equations in the var iables a coshn and

}i = cosO are Legendre equations. It is further shown there . that the eigenfunctions of Legendre ' s equation are the Legendre functions of the first

kind PR (V , the Legendre functions of the second kind Q R ( ^) and the associated

Legendre functions of the first and second kind PM ( E) and QQ ( E). The potential around a prolate spheroid polarized parallel to its major axis is

axially symmetric. Therefore only the Legendre functions of the first and

second kinds are necessary to construct a solution.

2.3.6.2 Matching Boundary Conditions at the Surface of the Spheroid

Two potential 'functions, ¢l inside the spheroid and ¢ 2 exterior to the

spheroid, must be :onstructed. All of the Legendre functions of the second

kind include the terms

QO(V 2 log \ ^-i (1)

or

QO(11) = 2 log (1

+u /

(2)

Note that both funccions in (1) and ( 2) approach - as their arguments

approach unity. Since the z -axis through the spheroid corresponds to 11=1,

the Qk (p) functions must be excluded from both ^l and ^ 2 . Similarly, the

z-axis inside the spheroid corresponds to ^=l, where
,.

Page 102

and therefore QR (e) must be excluded from 01.

Now that the pathological functions have been excluded we are ready to

construct well, behaved potential functions for ¢ 1 and 0 2 . Using all of the

allowed eigenfunctions of the axially symmetric Legendre equation we sec

that the general solution for Q 1 and ^2 is

W

^1 = E 

U 

ARPR(^)PR(P)

R=

0

2 E CBRPR (^) + CRQR (^ d PR(u)

R=0 -^

where

rl+r2

= coshn = `

c

r2 rl

u = cose = 2c

(4)

(5)

(6)

(7)

I

and r and r 2 are the magnitudes of

the spheroid to any point in space,

and n and 6 are the coordinates use

tem as shown in Section 2.3.2•

the vectors r and r2 from the two foci of

c is the semi focal length of the spheroid,

i to define the ellipsoidal coordinate sysThe unknown coefficients AR , BR and C will be ev ,,ivated by imposing three

boundary conditions:

1) At great distance from the spheroid the potential ^2 must approach the

potential of a uniform field. Since a potential function may 'have a constant

added without changing the associated field E because

— vO (s)

we are free to choose the origin of the spheroidal coordinate system as the

point of zero potential. The potential function which gives rise to a uniform

` field Eo in the z direction is

$ _ -zEo (9) j

r

t
Page 103

But from the transformation from prolate spheroidal to Cartesian coordinates

derived in Section 1, we see that equation (32) of Section 2.3,E is

z = c^u (10)

4 = cpl (0Pl (u) (11)

r

At great distances from the origin, Q2 (C) becomes negligibly small therefore

4 ( 5) simplifies to

W

0 2 = F B

k=0

RP R (O p i (11) = -cP l ( f)P l (P)E0

(12)

It is therefore clear that

B0=0

B1 = -CEO

(13)

BR = 0 for all R>l

2) At the interface between regions 1 and 2 ( the surface of the spheroid

where = g0 ) the potentials must be equal. From ( 4), (5) and ( 13) we see

that

m m

AtPR (Ypk (U) E CZQP. Q0 ) PR ( 11 ) - cEoPl (UX0

(14)

k=0 R=0

Since the Legendre Polynomials P R (p) are orthogonal over the interval -1 < u < 1

we can multiply ( 14) by a Legendre polynomial Pk (u) and integrate from -1 to +1

to separate the coefficients. This operation results in

A010 (c 0) - C0Q0 (c0 ) = 0 (15)

f A.
Page 109

A1^0 - C1Q1 Y _ -cEOEO (16)

Akpk (E 0 ) - CIA(

Y = 0 for all k>1 (17)

3) At the interface between the two regions 1 and 2 the normal derivatives

of the potentials multiplied by the relative nermeativit'ies of their respective

regions must be eq al. This is required to insure continuity of the no•,mal

component of flux across the boundary. Thus both Dirichlet and Neumann

l boundary conditions are imposed. Thus

a^ a^

C

r a^ a^ when l _ f0 (1&)

or from (4) and (5)

m m ;,

s r LE0 AtPX 0 )PR (U) = ky0 CBRPR(^0) + Cp4^(E0 )] PR (U) (19)

CS

e:

a

Again invoking the orthogonality of the Legendre polynomials over the interval

-1SUS1 and using (13) we can separate the coefficients with the result; . I

e r AOPO(C0 ) - COQ-(C 0) = 0 (20)

Cr Al - B 1P1(C0 ) - C1Qi(C0 ) 0 (21)

e r.AkPk(Y - CkQk(E 0 ) = 0 (22)

Equations (20) and (21) can be simplified by recalling that

P0 (0 = 1 PD(0 = 0 (23)

Pl (O = E P1M = 1 (24)
Y77

Pago ].05

Since PD(C)°0, we conclude that

CO ° 0 (25)

From (15) we see that A0.0

At this point we shall deal with all of the coefficients for k>1 recalling

(17) and (22)

kpk( Y - CkQk(C0) - 0 (17)

e r Aj'(^0 ) - Ckgao) - 0 (22)

Multiplying (17) by -'Qk(E 0) and (22) by Qk (C O ) and adding the two equacions,

we find that

Ak IF k(^0 )Qk (E 0) - Pk(EO)Qk-(EO )l = 0 (26)

If we can prove that the term in the brackets is never equal to zero we have

established that all A and C coefficients are zero for k>l. The Russian

mathematician N.N. Lebedev l shows that the Wronskian W(P R (V I VO ) which is

defined as

pk(O QRM

W 

IF 

v 'Qy (E] =

PR(E) Q,(9)

pk (mia) - PRmyc) (27)

1/(1-E 2 ). Therefore, equation (26) can be rewritten as

:. CEr -1)P k(C0 )Qk( g0 ) + U-J 0 (28)

0

^yl^
_

a

Page 106

Now C 0

>1 and the Q functions are always positive for E>l. Furthermore, all

of the roots of Pk(0 lie in the interval -1 < < 1. Thus Pk(E E 0 ) >0 for

^O>1. Since the term enclosed in the brackets is always greater than zero

the ,anly possible solution to (28) is

Ak - 0 for all k>1 (29)

from which we deduce from (17) that

C = 0 for all k>l. (30)

Thus we see that the only non-zero coefficients of (4) and (S) are A 1 , B1

and Cl . But it has already been shown by the first boundary condition that

B1 -cEo Therefore, equation (21) simplifies to

g er Al - C 1Qi(E0) = -cEo (31)

F

But

f 

Ql (E) QOM- Ell (32)

Substituting (32) into (31) results in

C^_

I

er Al - Cl 

QO - 2 _ cEo E=EO (33)

or

c

Cc Eo

1 _ F_- _ : Al

er CQO E z -1 Cr

E=E0 (34)

Substituting (34) into (16) {

i

l

C l z cEoE

e L QO - ' 2=1] - C1Q1 = - Er (eT-1) . E=E O

(35)

1

i

I

1

- M
But Page 107

Ql(4) a EQ0 (E) -1

Therefore (35) simplifies to

C l 1(erl)Q1 + -^ j _ EcEo(cr 1) E-EO (36)

Dividing (36) by (er -1) results in

cf,Eo

C1

1 - E-E0 (37)

Q l(E) + (Er -1) W7

Substituting (37) into (16) results in

cEoQl

Al - CEO + 1 E-90 (38)

Q1 + (er l jW-1)

which simplifies to

- CEO

A, E=E

1 + (Or-1)(EZ-1)Ql(E) 0

(39)

Thus we have evaluated A l and C1 in terms of c, Eo, er and E0 , and we can write

the solution for the potential in the two regions 1 and 2 as

^l A1E11 1

sE<E0 (40)

4 2 C 1Q 1 (E)u - CEO Ell

co 

<E<- (41)

But by the transformation which generated the prolate spheroidal coordinates

resulted in the relationship

z - c Eu (1-35)

Therefore we see that the field within the spheroid which is given by

E _ - VO1

is uniform and is given by

+ Al ,

E C z

(42)
a of a dielectric spheroid immersed in an otherwise

To specialize to the case of a conducting wire we

natant to go to infinity. Under this condition Al=0

e ,E,

O1

Q1 U.) (43)

(43) into ( 4) with the result

^^ = eE,I EOQ1(f - f] μ (44)

L Q1 (G) 1

18) and (19) from Section 2.3.2

ri = \/x 2 + y2 + (z-c) 2 (12)

r2 = X2 -+ Y 2 + (z+c) 2 (13)

rl + r2

2c

r2 - ri

μ 2c

s

- _R iatr^r

'J

++' Page 108

G, II

i

t^
Page 109

we see that for every point x, y, z in space ue can calculate the electric

potential 0. This equation is used to determine the boundary conditions

for the numerical solution of the electric field around the wire in the

presence of plasma. A plot of the electric field around a conducting

spheroid is shown in Figure 2.3.3. The electric flux lines were calculated

from an analog of the magnetic vector potential

D=VxA. (45)

A smaller plot for a thin wire is shown in Figure 2.3.4. Equation (44) is

also used to initialize all of the ,interior grid points before the iterative solution is begun. The details of this variable grid iterative solution will be described in the next report.

i
Page 110

Figure 2.3.3 Electric field around a conducting

spheroid.
OLE

R

Page 111

Figure 2.3.4 Electric field around a conducting

thin tether.
Page 112

REFERENCES

1. P. Moon and D.E. Spencer; Field Theory for Engineers, 1961; D. Van Sostrand

Company, Inc., p. 348.

2. M. Harrold; Calculation of Equipotentials and Flux Lines in Axially Symme-

: LGEL transactions an

3. P. Moon and D.E. Spencer; Field Theory for Engiceers, 1961; D. Van Sostrand

Company, Inc., Chapter 9.

4. John David Jackson; Classical Electrodynamics, 2nd Edition. 1975; John

Wiley and Sons, Inc., p. 87.

5. N.N. Lebedev; Special Functions and Their App Licatfons (translated by

Richard Silverman); 1972; Dover Publications; P. 167,

6. Ibid., p. 183.

7. Ibid., p. 167.

S. I.S. Gradshteyn and I.M. Ryzhik; Table of Integrals, Series, and Products

(translated by Scripts Technics, Inc.); 1965, Academic Press; p. 1C^5.

9. A. Erdely L, W. Magnus, F. Oberhettinger and F.G. Tricomi; Higher Transcendental Functions, Volume I; 1953; McGraw-Rill Book Company, Inc., p. 145.

10. Handbook of Mathematica.l Functions with Formulas, Graphs, and Mathematical

Tables; Edited by M. Abramowitz and I.A. Stegun; 1965; National Bureau of

Standards, p. 334.

11. N.N. Lebedev; Special Functions and Their Applications (translated by

Richard Silverman); 1972; Dover Publications, p. 185.

12. P. Moon and D.E. Spencer; Field Theory for Engineers, 1961; D. Van Nostrand

Company, Inc, p. 137.

13. C.N. Hardy; A Course of Pure Mathematics, 10th Edition; 1967; Cambridge

University Press. p. 229.

14. N.N. Lebedev; Special Functions and Their Applications (translated by

Richard Silverman); 1972; Dover Publications, p. 165.

15. P.M. Morse and H. Feshbach; Methods of Thec[etical Physics; 1953; McGrawHill Book Company, Inc., p. 401.

16. Ibid., p. 547,

17. L.J. Chu and J.A. Stratton; Elliptic and Spheroidal Wave Functions; Journal

of Mathematics and Physics Vol 20; 1941, pp. 259-309.

18. N.N. Lebedev; Special Func tlona and Their Appllcatio^ta (translated by

Richard Silverman); 1972; Dover Publications, p. 182.

19. Ibid., p. 178.

20. Ibid., p. 4.

21. Ibid., p. 12.

22. A.I. Markushevich; Theory of Functions, Vol. I; 1965; Chelsea Publishing

Company, p. 334.

23. Ibid., p. 404.

I

i

i

w

3

1

^+ f

o ^

t
Page 113

F

P N

2.4 Potential Build-Up Due To Plasma Contactor Failure

In an easterly orbit, with the satellite deployed upward, the satellite collects electrons and the Shuttle collects ions. Because of the low

thermal velocity of the ions, it is necessary to have ^ plasma contactor on

board the Shuttle in order to maintain it at plasma potential. If the

plasma contactor fails, the Shuttle will acquire a negative potential as a

result of the inability of the Orbiter to neutralize the electrons coming

down the wire from the satellite. If the satellite has a plasma contactor

It will be a better collector and the Shuttle will acquire a large potential.

Two simulations have been run in order to estimate the potential ac-'

quired by the Shuttle if its plasma contactor fails. These simulations are

relevant to the first electrodynamic mission with the Shuttle on an orbit 7

'i

i

of 28.50 inclination and 295 km altitude. The satellite is deployed upward j

is 1,

on a 20 km long tether. These simulations refine and confirm the value of f

-2kvolt for the Shuttle potential as it was estimated in the Quarterly

Report #1. In the first one the satellite has no plasma contactor and is',

modeled as a 1.5 meter dismmeter sphere. The Shuttle is modeled as a 12.4

meter diameter sphere. The wire is 20 km long and has a resistance of 4000

I

ohms. The Shuttle acquired a voltage of about -1195 volts and the satellite was at +1206 volts. The wire current was .253 amps. In the second

f

simulation, the satellite has a plasma contactor and the other parameters

are the same. The Shuttle acquired a potential of -2064 volts and the

satellite was within a volt of plasma potential. The current in the wire

was .337 amps.
Page 114

3.0 PROBLEMS ENCOUNTERED DURING REPORTING PERIOD

None.

4.0 ACTIVITY PLANNED FOR THE NEXT REPORTING PERIOD

During the next reporting period we will carry out the software modification of the force and torque package in our rotational dynamics computer code. The implementation of the reeling control laws will then allow

us to simulate actual maneuvers with a good fidelity. Our plotting program

will also be modified because presently it is not formatted to receive and

interpret the rotational dynamics variables.

The three-dimensional slack tether computer code SLACKS is now available and tested. During the next reporting period it will be used to

simulate actual cases of tether break and to investigate suitable Shuttle

avoidance maneuvers. In the mean time the analytic studies on the slack

tether will continue. The major goal is the evaluation of the wrinkle's

typical scale size at the onset of the compressional instability.

The development of the computer code to evaluate the electric field

around the severed tether is in progress. The next reporting period will

be devoted to the implementation of the electric potential variations due

to the plasma. Presently the program can evaluate the electric potential

around a charged body in vacuo.

i