[PDF] - ( ) R O , i , j , k is an Earth centered equatorial PDF Document

SLIDE 1

Dynamics_revised.doc data creazione 12/ 07/ 2001 9.35.00 Pagina 1 di 18 data ultima revisione: 29/ 04/ 2002 8.45.00

Dynamic Model of the Spacecraft Position and Attitude

Basilio BONA, Enrico CANUTO Dipartimento di Automatica e Informatica, Politecnico di Torino Corso Duca degli Abruzzi 24 10129 Torino, Italy

tel. 011 564 7026, fax 011 564 7099

bona@polito.it

1 Reference Frames

1.1 Reference Frames Description

In this report the following reference frames will be considered:

1.1.1 Inertial Reference Frame – IRF

( )

J2000 J2000 J2000 J2000 J2000

, , , O R i j k

is an Earth centered equatorial inertial frame, with:

J2000

O

at the centre of the Earth.

J2000

i

along the intersection of the mean ecliptic plane with the mean equatorial plane, at the date of 01/01/2000; positive direction is towards the vernal equinox.

J2000

k

rthogonal to the mean equatorial plane, at the date of 01/01/2000; positive direction is

towards the north.

J2000

j

completes the reference frame.

1.1.2 Earth-fixed Reference Frame – EFRF

( )

E E E E E

, , , O R i j k

is an Earth fixed non-inertial frame, and shall be the most recently defined International Terrestrial Reference Frame (ITRF). We assume:

E

O at the centre of the Earth.

E

i along the intersection of the equatorial plane and the Greenwich meridian plane; positive

direction from the earth centre to the 0° longitude point on the equator.

E

k orthogonal to the equatorial plane; positive direction towards the north pole.

E

j completes the reference frame. 1.1.3 Orbit Plane Reference Frame – OPRF

( )

OP OP OP OP OP

, , , O R i j k

is a local non-inertial orbit plane reference frame, with:

OP

O

at the centre of the Earth.

OP

i

along the intersection of the orbit plane with the equatorial plane.

OP

k

btained by rotating

E

k positively around OP i

by the orbit inclination angle i .

SLIDE 2

Dynamics_revised.doc data creazione 12/ 07/ 2001 9.35.00 Pagina 2 di 18 data ultima revisione: 29/ 04/ 2002 8.45.00

OP

j

completes the reference frame.

1.1.4 Local Orbital Reference Frame – LORF

( )

O O O O O

, , , O R i j k

is a non-inertial local orbital reference frame, with:

O

O at the actual satellite centre of mass.

O

i (roll) parallel to the instantaneous direction of the orbital velocity vector = v r ; positive

direction towards positive velocity; O = v

i v

.

O

j (pitch) parallel to the instantaneous direction of the orbital angular momentum = × h r v , where r is the vector from the Earth centre to

O

O ;

O

× = × r v j r v

.

O

k (yaw) completes the reference frame:

O O O

× = = × × v h k i j v h

.

1.1.5 Spacecraft Reference Frame – SCRF

( )

SC SC SC SC SC

, , , O R i j k

is a local satellite non-inertial frame, with:

SC

O

at the centre of the launcher-satellite adaptor interface plane

SC

i

along the launch vehicle axis; positive direction is towards the launch vehicle nose.

SC

k

rthogonal to the satellite earth face; positive direction is towards nadir.
SC

j

completes the reference frame.

1.1.6 Attitude Control RF – ACRF

( )

AC AC AC AC AC

, , , O R i j k

is a local non-inertial satellite reference frame, with:

AC

O

at the actual satellite centre of mass (COM).

AC

AC AC

, , i j k

parallel to SC

SC SC

, , i j k

.

1.2 Transformation Matrices between Reference Frames

In the following we use the notations

( ) ( ) ( )

, , , , , α β γ R i R j R k

to represent, respectively, rotations around the X, Y, and Z axis. We remember that a rotation about fixed axes pre-multiplies, while a rotation about the moving axes post-multiply. a)

J2000 OP

⇒ R R

, given by matrix

J OP

R

:

( ) ( )

J OP

c s c s s , , s c c c s s c

i i i i i i

i

Ω Ω Ω Ω Ω Ω

⎡ ⎤ − ⎢ ⎥ ⎢ ⎥ = Ω = − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ R R k R i

(1.1)

SLIDE 3

Dynamics_revised.doc data creazione 12/ 07/ 2001 9.35.00 Pagina 3 di 18 data ultima revisione: 29/ 04/ 2002 8.45.00

notice that, while the inclination angle i is constant (see Table 2), the right ascension Ω varies at a rate of 0.9856 deg/day (

7

1.996425 10 ⋅

rad/s), to allow the satellite to be sun-syncronous b)

J2000 E

⇒ R R

, given by matrix

J E

R :

( )

J E

c s , s c 1

δ δ δ δ

δ ⎡ ⎤ − ⎢ ⎥ ⎢ ⎥ = = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ R R k

(1.2) notice that the angle δ varies at a rate of 360 deg/day (

5

7.2921 10 ⋅

rad/s), due to daily rotation. c)

J2000 O

⇒ R R

, given by matrix

J O

R :

J O

⎡ ⎤ × × × ⎢ ⎥ = ⎢ ⎥ × × × ⎢ ⎥ ⎣ ⎦ v r v v r v R v r v v r v

(1.3) where v and r are defined in Section 1.1.4. d)

E OP

⇒ R R

, given by matrix

E OP

R

:

( )

E E J J J OP J OP E OP

c s 0 c s c s s c c s s c s c s c c c s s s c s s c 0 s c c c s s c c s s s c c c c s s s c c s 1 s c s c

i i i i i i i i i i i i i i i i δ δ δ δ δ δ δ δ δ δ δ δ δ δ δ δ Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω

′ = = = ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ − + − + − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − − = − + + − − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ R R R R R

(1.4) e)

OP O

⇒ R R

, given by matrix

OP O

R

:

( ) ( )

OP O

s c , 90 , 90 c s 1

α α α α

α ⎡ ⎤ − ⎢ ⎥ ⎢ ⎥ = + = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦

R

R i R j

(1.5) where α is the true anomaly. f)

O AC

⇒ R R

, given by matrix

O AC

R

:

( ) ( ) ( )

O AC

c c c s s s c c s c s s , , , s c s s s c c s s c c s s c s c c

ψ θ ψ θ ϕ ψ ϕ ψ θ ϕ ψ ϕ ψ θ ψ θ ϕ ψ ϕ ψ θ ϕ ψ ϕ θ θ ϕ θ ϕ

ψ θ ϕ ⎡ ⎤ − + ⎢ ⎥ ⎢ ⎥ = = + − ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ ⎣ ⎦ R R k R j R i

(1.6) where ϕ = roll angle, θ = pitch angle, ψ =yaw angle; for small angles the above matrix reduces to:

O AC

1 1 1 ψ θ ψ ϕ θ ϕ ⎡ ⎤ − ⎢ ⎥ ⎢ ⎥ = − ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ ⎣ ⎦ R

(1.7)

SLIDE 4

Dynamics_revised.doc data creazione 12/ 07/ 2001 9.35.00 Pagina 4 di 18 data ultima revisione: 29/ 04/ 2002 8.45.00

2 Dynamic State Space Model

2.1 Assumptions

The state space model derived in this report describes a) the dynamics of the S/C COM position with respect to the inertial frame

J2000

R

(orbital dynamics); b) the dynamics of the S/C reference frame

AC

R

with respect to the orbital frame

O

R

(attitude dynamics). The following assumptions hold: Assumption 1: the dynamics of the COM position is described by the differential equations of the

AC

R

rigin represented in

J2000

R

by the vector

x y z ′ ⎡ ⎤ = ⎢ ⎥ ⎣ ⎦ r

. The COM velocity origin is represented in

J2000

R

by the vector

x y z

v v v ′ ⎡ ⎤ = ⎢ ⎥ ⎣ ⎦ v

. Assumption 2: the attitude of

AC

R

with respect to

O

R

is given by the three roll-pitch-yaw angles

α = ϕ θ ψ ′ ⎡ ⎤ ⎢ ⎥ ⎣ ⎦ implicitely defined in eqn. (1.6).

Assumption 3: in order to transform the vectors represented in a generic referecne frame

A

R

into vectors represented in another reference frame

B

R

, the following relation holds:

B A

⎡ ⎤ ⎡ ⎤ = ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦

R R

p R p

(1.8) where

B

⎡ ⎤ ⎢ ⎥ ⎣ ⎦R p

represents a generic vector in the reference frame

B

R

and

A B

R is a orthonormal matrix with

positive unitary determinant (proper rotation). Assumption 4: inputs and/or measurements, defined or obtained in reference frames other than the reference frame

AC

R

, are assumed to be already represented in this reference frame. Assumption 5: the axes

AC AC AC

, , i j k

f the A/C reference frame are assumed to coincide with the principal

axes of inertia. Assumption 6: the S/C mass m and inertia matrix

x y z

J J J ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ J

(1.9) are assumed to be time-invariant, although this is not true, due to propellant consumption.

2.2 Plant Dynamic Equations

In order to model the S/C dynamics it is necessary to introduce first the orbital kinematic equations.

2.2.1 Kinematic equations

Assumption 7. The total inertial velocity vector ω is the sum of the orbital rotation velocity and the S/C attitude rate.

SLIDE 5

Dynamics_revised.doc data creazione 12/ 07/ 2001 9.35.00 Pagina 5 di 18 data ultima revisione: 29/ 04/ 2002 8.45.00

To this end let us denote

1. with

( )

O O O O O

, , , O R i j k

the local orbital frame,

2. with

( )

O O O O O

, , , O ′ ′ ′ R i j k

the same frame rotated around the z-axis by the angle ψ,

3. with

( )

O O O O

, , ,

AC

O ′′ ′ ′ R i j k

the previous frame rotated around the y-axis by the angle θ,

4. with

( )

O,

, ,

AC AC AC AC

O R i j k

the attitude frame, obtained by rotating the previous frame around the x-axis by the angle ϕ,

5. with

2 O O O

ω × = = ω r v j r

the orbital rotation velocity. Then we can write

O O O O AC

ω ψ θ ϕ ′ = + + +

ω

j k j i

Upon transformation of the vector components in the attitude reference frame, it results

( ) ( ) ( ) ( )

, , , ,

x y O z

ϕ ω ω ϕ θ ψ ω ϕ θ ω ψ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = = − − − + − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦

ω

R i R j R k R i

and explicitly

s c s s s s c c c s s s c c s c c

x y O z

c s

ψ θ θ ψ θ ϕ ψ ϕ θ ϕ ϕ ψ θ ϕ ψ ϕ θ ϕ ϕ

ϕ ω ω ω ψ θ ω θ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = + + + ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦

(1.10)

The kinematics non linear equations are therefore

( ) ( ) ( ) ( ) ( ) ( )

2 2

1

x O y O z O y z O z y z O O z y O z y O x O

s s c c s c c s s s c c c s c s c s c s s s c c s c s c c c s c s c s s c c s c s s s c s c s s s s c c

θ ψ θ ϕ ϕ ϕ ϕ ψ θ ϕ ψ ϕ ϕ ϕ θ ϕ ϕ ϕ θ ϕ ψ θ ϕ ψ ϕ ϕ ϕ ψ θ ϕ ϕ ϕ ϕ ψ ψ θ ϕ ψ ϕ ϕ ϕ ψ θ θ θ ϕ ϕ ψ θ ψ θ θ

ϕ ψ ω ω θ θ ω ω ψ ω ψ ω θ ω ω ω ψ ω ω ω ω ω ψ ω ω ω ϕ ω ω ω ω ω = + − + = − + − − + + − = − − = + − − − − = + − = + − + −

(

)

x z y O

s s c s c c

ψ θ ϕ ϕ θ θ

ω ω ω ω = + + −

(1.11) The final form of non linear equations is the following

( ) ( )

1

x z y O y z O z y O

s s c s c c c s c c s s s c

ψ θ ϕ ϕ θ θ ϕ ϕ ψ ϕ ϕ ψ θ θ

ϕ ω ω ω ω θ ω ω ω ψ ω ω ω = + + − = − − = + −

(1.12)

SLIDE 6

Dynamics_revised.doc data creazione 12/ 07/ 2001 9.35.00 Pagina 6 di 18 data ultima revisione: 29/ 04/ 2002 8.45.00

Because attitude will be controlled, we will adopt a quasi linear form, by neglecting terms of the third order and more:

x z O y z O z y O

ϕ ω θω ψω θ ω ϕω ω ψ ω ϕω θψω = + − = − − = + −

(1.13)

2.2.2 Orbital Dynamics

The differential equations describing the S/C COM dynamics are given by the Newton equation:

( ) ( ) ( )

total t

m t m t + − =

F

g a

(1.14) where

total k k

= ∑ F F is the vector sum of all disturbance and command forces, except those due to gravity

and inertia, acting on the S/C COM, g is the local gravity acceleration vector acting on the S/C COM,

=

a

v and v are the S/C COM linear acceleration and velocity respectively, and m is the S/C mass.

We note that acceleration

( )

t

a

can be decomposed into two terms: a) Centripetal acceleration

( )

c t

a

, due to the S/C orbital angular velocity

2

O O

T π ω =

, with

,nominal

5370 s.

O

T

This acceleration is given by

( )

c O O

= × × ω ω a r and gives origin to a

centrifugal force

( )

c c

t m = − F a .

The vector

r = − ρ r

gives the S/C position relative to the Inertial Reference Frame

J2000

R

, where

ρ is the radial versor and nominal

S

r R h

⊕

= + 6628km;

O

ω has been previously defined.

For a circular orbit

O

= − ρ k .

b) Residual acceleration

( ) ( ) ( )

res c

t t t = −

a

a a

. We further assume that the external forces acting on the S/C consist only in:

1. Drag forces

( )

d t

F

, due to aerodynamic and magnetic drag; we assume to possess the drag forces computed from GMV simulations and represented in

O

R

; it is threfore necessary to represent them in the inertial reference frame, as follows:

( ) ( ) ( ) ( ) ( ) ( ) ( )

J2000 O

, J O , , dx dx GMV d dy dy GMV dz dz GMV

F t F t t F t F t F t F t ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦

R R

F R

(1.15) where

J O

R is given in eqn. (1.3).

2. Command forces due to actuation of nine thrusters

0, 1, , 8 j = …

whose thrusts are collected in a vector

( )

t t

u

:

SLIDE 7

Dynamics_revised.doc data creazione 12/ 07/ 2001 9.35.00 Pagina 7 di 18 data ultima revisione: 29/ 04/ 2002 8.45.00

,AC 1 AC ,AC ,AC 8

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

t x t y t t t z t

u t F t u t t F t t F t u t ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = = = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎢ ⎥ ⎣ ⎦

F

V u V

(1.16) where

t

V is the action matrix to be defined later.

Assumption 8: we assume that the command forces are represented in the attitude control reference frame

AC

R

, so it is further necessary to transform them into inertial frame as follows

( ) ( )

J AC AC

( ) ( ) ( )

x y z

F t t F t t F t ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ = = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ F R F

(1.17) where

J J O AC O AC

= R R R

(1.18) The gravity vector ( )

t g

, due to zonal harmonics induced by Earth oblateness, can be approximately written as:

( ) ( ) ( )

2 2 2 3 3

3 , 5sin 1 2 sin 2

res

R J r r r r r μ μ δ δ δ

⊕ ⊕ ⊕ ⊕

⎛ ⎞ ⎛ ⎞ ⎟ ⎜ ⎟ ⎟ ⎜ ⎜ ⎟ ⎟ ⎜ ⎜ − − − − = − + ⎟ ⎟ ⎜ ⎜ ⎟ ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠ ⎟ ⎜ ⎝ ⎠

g

r r k r g

(1.19) where

⊕

k is the earth pole axis.

Finally, eqn. (1.14) can be rewritten as:

( ) ( ) ( ) ( ) ( )

3 d res c res

t t m m t m t m t r μ⊕ + − + − − = F F r g a a

(1.20) Assuming a nominal orbital motion with angular velocity

O

ω , recalling the triple product rule

( ) ( )

× × = ⋅ − ⋅ a b c a c b a b c

(1.21) where

⋅ a b denotes the scalar product ′ ⋅ = a b a b , and noting that ⊥ r

Ο

ω

, it follows that the centripetal force:

( )

2 3 c c O O O

m m m m r μ⊕ = − = − × × = = ω ω ω F a r r r

(1.22) where

2

S

T π =

Ο

ω

and

3

2

S

r T π μ⊕ =

, we have a nominal equilibrium between

c

F and the first right-side

term of (1.19) Then, eqns. (1.14) and (1.20) can be rewritten as

SLIDE 8

Dynamics_revised.doc data creazione 12/ 07/ 2001 9.35.00 Pagina 8 di 18 data ultima revisione: 29/ 04/ 2002 8.45.00

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

, , , , , ,

( )

d res res dx x res x res x dy y res y res y dz z res z res z

t t m t m t F t F t a t g t F t F t m a t m g t F t F t a t g t + − + = ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ + − + = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ F F a g

(1.23) so that

, , , , , ,

1 1 1

res x dx dy dz x res x res y dx dy dz y res y res z dx dy dz z res z

a F F F F g m a F F F F g m a F F F F g m ⎡ ⎤ = + + + + ⎢ ⎥ ⎣ ⎦ ⎡ ⎤ = + + + + ⎢ ⎥ ⎣ ⎦ ⎡ ⎤ = + + + + ⎢ ⎥ ⎣ ⎦

(1.24)

2.2.3 Attitude Dynamics

The differential equations describing the S/C attitude dynamics are given by the Euler equation:

( ) ( ) ( )

( )

( ) ( ) ( )

k k k a k k

t t t t t t + × − − × =

∑ ∑

ω ω C b F Ja J

(1.25) where

dk k

∑C

is the vector sum of all disturbance and command torques, except those due to gravity and inertia, acting on the S/C COM,

k

b is the vector of application of

dk

F , J is the inertia matrix, ω is the S/C

inertial angular velocity, and

a

a is the S/C reference frame angular acceleration. Contrary to the orbital

dynamics equation, here all vectors are represented in

AC

R

(see also Section 2.2.1). The following assumptions hold: Assumption 9: the only torques acting on the S/C are due to aerodynamic drag, magnetic drag and thrusters

actuation. Gravity gradient torques are assumed negligible.

Assumption 10: gravitational forces, centrifugal forces

( )

c t

F

and drag forces

( )

d t

F

have application points coinciding with the S/C COM, so they do not give origin to torques. With these assumptions, eqn. (1.25) reduces to

d a

+ − − × = ω ω C C Ja J

(1.26) where C is the command torque vector due to thruster actuation ( )

t u

(already represented in

AC

R

, as in Assumption 8):

( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )

1 8 t x t y tc t tc z t

u t C t u t t C t t C t u t ⎡ ⎤ ⎢ ⎥ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = = = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎢ ⎥ ⎣ ⎦

C

V u V

(1.27) where the matrix

tc

V will be defined later.

As with the drag forces, the drag torques

d

C are obtained from GMV simulation, and represented in

O

R

, so now it is necessary to trasform the as follows:

SLIDE 9

Dynamics_revised.doc data creazione 12/ 07/ 2001 9.35.00 Pagina 9 di 18 data ultima revisione: 29/ 04/ 2002 8.45.00

( ) ( ) ( )

( )

( ) ( ) ( ) ( ) ( ) ( )

AC O O

,GMV O AC ,GMV ,GMV

c c s c s c s s s c s s s c c c s c s c +s s s s c c s c c

dx dx dx dy dy dy dz dz dz

C t C t C t C t C t C t C t C t C t

ψ θ ψ θ θ ψ θ ϕ ψ ϕ ψ θ ϕ ψ ϕ θ ϕ ψ θ ϕ ψ ϕ ψ θ ϕ ψ ϕ θ ϕ

⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ′ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = = − + ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦

R R R

R

(1.28) Then, eqn. (1.26) can be rewritten as:

1 a d − ⎡

⎤ = + − × ⎢ ⎥ ⎣ ⎦ ω ω a J C C J

(1.29) Since

( ) ( ) ( )

,

y z y z z y x x z x y y z x x z y x z z x y x y

J J J J J J J J J ω ω ω ω ω ω ω ω ω ω ω ω ω ω ω ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ − − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ × = − = − − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ω ω J

(1.30) it follows:

( ) ( ) ( )

1

c c s c s c s s s c s s s c c c s c s c +s s s s c c s c c

y z y z ax dx x ay dy y z x x z az dz z x y x y

J J a C C a C C J J a C C J J

ψ θ ψ θ θ ψ θ ϕ ψ ϕ ψ θ ϕ ψ ϕ θ ϕ ψ θ ϕ ψ ϕ ψ θ ϕ ψ ϕ θ ϕ

ω ω ω ω ω ω

−

⎧ ⎡ ⎤ ⎪ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ ⎡ ⎤ − − ⎪ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎪ ⎪ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎪ = − + + + − ⎨ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎪ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎪ − ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎣ ⎦ ⎩ J ⎫ ⎪ ⎪ ⎪ ⎪ ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎭

(1.31)

2.2.4 State Space Equations

The components

i

x of the system state are assumed to be the 6 linear positions and velocities of S/C COM

and the 6 angles and angular velocities around the three

AC AC AC

, , i j k

axes:

1 2 3 4 5 6 7 8 9 10 11 12

;

x y z

x x x y x z x x x y x z x x x x x x ϕ θ ψ ω ω ω ⎡ ⎤ = ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎣ ⎦

x

(1.32) From (1.13), (1.24), (1.31) and (1.32) we can write the differential equations in state variables form. Time dependence is omitted for brevity.

SLIDE 10

Dynamics_revised.doc data creazione 12/ 07/ 2001 9.35.00 Pagina 10 di 18 data ultima revisione: 29/ 04/ 2002 8.45.00

1 4 2 5 3 6 4 , 5 , 6 , 7 10 8 12 9 8 11 7 12 9 12 7 11 8 9 10 8 9 8 9 8

1 1 1 1 c c c s s

dx dy dz x res x dx dy dz y res y dx dy dz z res z O O O dx dy dz x

x x x x x x x F F F F g m x F F F F g m x F F F F g m x x x x x x x x x x x x x x x x C C C J ω ω ω = = = ⎡ ⎤ = + + + + ⎢ ⎥ ⎣ ⎦ ⎡ ⎤ = + + + + ⎢ ⎥ ⎣ ⎦ ⎡ ⎤ = + + + + ⎢ ⎥ ⎣ ⎦ = + − = − − = + − = + −

(

) ( ) ( ) ( ) ( ) ( ) ( )

11 12 11 7 8 9 7 9 7 8 9 7 9 7 8 10 12 12 7 8 9 7 9 7 8 9 7 9 7 8 10 11

1 s s c c s s s s c c s c 1 c s c +s s c s s s c c c

y z x dx dy dz z x y y dx dy dz x y z z

J J x x C x C C C J J x x C J x C C C J J x x C J ⎡ ⎤ + − + ⎢ ⎥ ⎣ ⎦ ⎡ ⎤ = − + + + + − + ⎢ ⎥ ⎣ ⎦ ⎡ ⎤ = + − + + − + ⎢ ⎥ ⎣ ⎦

(1.33)

where

( )

c cos , s sin

i i i i

x t x t = =

.

3 Thruster command

3.1 Line of actions

Let us denote with

tj

v the unit vector of the mean line of action of a thruster j, defined as the mean thrust

direction. It is defined by the following data in case of the ion thruster (ITA):

α x y β vt0 α x π/2+β vt β y

SLIDE 11

Dynamics_revised.doc data creazione 12/ 07/ 2001 9.35.00 Pagina 11 di 18 data ultima revisione: 29/ 04/ 2002 8.45.00

Figure 1 – Angular coordinates of the mean lines of actions of ITA and MTA

cos cos sin cos sin 0.049 0.009 (3 ) 0.009 (3 )

t

α β α β β α α β β α σ β σ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ = ± + Δ = + Δ Δ ≤ Δ ≤ v

(1.34) and by the following data in case of micro thrusters (MTA):

sin cos cos cos sin 0.35 , 1,4 0.35 , 2,3 0.35 , 5,8 0.35 , 6,7 0.35, 0.005 (3 ) 0.005 (3 )

j j tj j j j j j j j j j j j j j j j j

j j j j β α β α α α π α α π α α π α α α α α α β β α σ β σ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ = − + Δ = − + Δ = = + + Δ = = − + Δ = = + Δ = = = + Δ Δ ≤ Δ ≤ v

(1.35) Remarks.

Note that the direction uncertainty in case of micro thrusters has not yet been frozen.

In the design case, the action matrix

t

V assumes the following form:

SLIDE 12

Dynamics_revised.doc data creazione 12/ 07/ 2001 9.35.00 Pagina 12 di 18 data ultima revisione: 29/ 04/ 2002 8.45.00

sin cos sin cos

t

c V s c c c c c c c c s s s s s s s s s c s c α α α α ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ = ± − − − − ⎢ ⎥ ⎢ ⎥ − − − − ⎢ ⎥ ⎣ ⎦ = = = =

(1.36)

3.2 Application points

Denote with

j

a the application point vector of thruster j. In case of ion thruster it holds 2.1 2.1tan

x x x y x y z z x y z

a a a a a b a a a a a a a α δ δ δ ⎡ ⎤ = − + Δ = − + Δ ⎢ ⎥ ⎢ ⎥ = = + Δ = + Δ ⎢ ⎥ ⎢ ⎥ = Δ ⎢ ⎥ ⎣ ⎦ Δ ≤ Δ ≤ Δ ≤ ∓ ∓ a

(1.37) and in case of micro thrusters it holds

, , 0.003 (3 ) 0.003 (3 ) 0.003 (3 ) 1.76m 2.08m 0.55m 0.18m

xj xj xj xj j yj yj zj zj xj yj zj

a a a a b a a e a a d a a a a a b c e σ σ σ ⎡ ⎤ = + Δ = − + Δ ⎢ ⎥ ⎢ ⎥ = = ± + Δ ⎢ ⎥ ⎢ ⎥ = ± + Δ ⎢ ⎥ ⎣ ⎦ Δ ≤ Δ ≤ Δ ≤ = = = = a

(1.38) In the design case, the application matrix

tc

V assumes the following form:

tc

a a a b b a a b b b e e e e e e e e d d d d d d d d ⎡ ⎤ − − − − − ⎢ ⎥ ⎢ ⎥ = − − − − ⎢ ⎥ ⎢ ⎥ − − − − ⎢ ⎥ ⎣ ⎦ ∓ V

(1.39)

3.3 Command matrix

The command matrix

t

B relating the thrust vector

t

u to force/torque vector F in spacecraft coordinates

holds

SLIDE 13

Dynamics_revised.doc data creazione 12/ 07/ 2001 9.35.00 Pagina 13 di 18 data ultima revisione: 29/ 04/ 2002 8.45.00

1 8 1 1 8 8 1 8 t t t t t tc t t t x y t z t t t t x y t z

F F u F u C C u C ⎡ ⎤ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ = = ⎢ ⎥ ⎢ ⎥ × × × ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = = = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦

V

v v v B V a v a v a v F B u B

(1.40) In the design case, it holds

t

c s c c c c c c c c s s s s s s s s dc es dc es dc es dc es dc es dc es dc es dc es as as bs bs as as bs bs ac ac bc bc ac ac bc bc ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ ± − − − − ⎢ ⎥ ⎢ ⎥ − − − − ⎢ ⎥ = ⎢ ⎥ − − + − + − − + − − − + ⎢ ⎥ ⎢ ⎥ − − − − ⎢ ⎥ ⎢ ⎥ − − − − ⎢ ⎥ ⎣ ⎦ B

(1.41)

4 Orbital Parameters

Let now define the Kepler orbital parameters in order to be able of derive from them the six position and velocity components and vice-versa.

4.1 Kepler Orbital Parameters

Six parameters define an ellyptical orbit according to the Kepler laws:

1. the semi-major axis a of the orbit ellipse;
2. the orbit eccentricity

2 2

1 b a ε = −

, where b is the semi-minor axis;

3. the inclination of the orbital plane with respect to the equatorial plane i;
4. the angle of right ascension of the ascending node (also known as RAAN) Ω , i.e. the angle on the

equatorial plane from the vernal equinox to the ascending node. The ascending node lies on the intersection between the equatorial plane and the orbital plane, in correspondence to the passage of the orbit from southern to northern emisphere;

5. the argument of the perigee ω , i.e. the angle from the ascending node to the perigee, measured on

the orbital plane:

6. the epoch of passage from the perigee 0

t , i.e. a reference starting time.

Fig. 1 depicts some of these parameters.

SLIDE 14

Dynamics_revised.doc data creazione 12/ 07/ 2001 9.35.00 Pagina 14 di 18 data ultima revisione: 29/ 04/ 2002 8.45.00

Fig. 1: the orbital parameters.

It is helpful to define now an auxiliary orbital plane reference

P

R

systems wich differs from

OP

R

since its x axis is aligned with the perigee and varies according to its rotation:

P

R

is defined as follows

( )

P P, , ,

O R p q w is a local non-inertial orbit plane reference frame, with:

P

O at the centre of the Earth.

p points toward the perigee of the orbit plane.
q is perpendicular to P

i and corresponds to a true anomaly angle ν (see Fig. 2) of 90 .

w completes the reference frame.

The transformation

J2000 P

⇒ R R

is given by the following matrix

( ) ( ) ( ) ( )

J J P OP

, , , , c c s c s s c s c c s s s c c c s s s +c c c c s s s s c c

i i i i i i i i i

i

ω ω ω ω ω ω ω ω ω ω

ω ω

Ω Ω Ω Ω Ω Ω Ω Ω Ω Ω

= Ω = = ⎡ ⎤ − − − ⎢ ⎥ ⎢ ⎥ ⎡ ⎤ = = + − − ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ R R k R i R k R R k p q w

(1.42)

4.1.1 From Kepler Parameters to Position and Velocity

The position r can be computed from the Kepler parameters in the following way:

( )

2

cos 1 sin a E a E ε ε = − + − r p q

(1.43) where E is the eccentric anomaly angle (see Fig. 2), obeys the differential equation

( )

3

1 cos GM E E a ε

⊕

− =

(1.44)

which has a solution

SLIDE 15

Dynamics_revised.doc data creazione 12/ 07/ 2001 9.35.00 Pagina 15 di 18 data ultima revisione: 29/ 04/ 2002 8.45.00

( ) ( ) ( )

3

sin GM E t E t M t t a ε

⊕

− = + −

(1.45) Fig.2 : the orbital plane and its parameters where

M is the mean anomaly at an arbitrary instant of time t . Notice that the mean anomaly M changes

by 360° during one revolution, but, in contrast to the true and eccentric anomalies, increases uniformly with time. The velocity

= v r can be computed as:

( )

2

sin 1 cos GM a E E r ε

⊕

= − + − v p q

(1.46) where r = r .

4.1.2 From Position and Velocity to Kepler Parameters

We assume to know the components of the versor w from (1.42):

x y z

w w w ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ = ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ w

The first two keplerian parameters i and Ω are computed as:

( )

2 2

atan2 , ; atan2 ,

x y z x y

i w w w w w = + Ω = −

(1.47)

SLIDE 16

Dynamics_revised.doc data creazione 12/ 07/ 2001 9.35.00 Pagina 16 di 18 data ultima revisione: 29/ 04/ 2002 8.45.00

where

1

atan2 , tan 90 if 0; 90 180 if 0; 180 90 if 0; 90 if 0; y y x x x y x y x y x y θ θ θ θ θ

− ⎛ ⎞

⎟ ⎜ ⎟ = = = ⎜ ⎟ ⎜ ⎟ ⎜ ⎝ ⎠ ⎧ ⎪ ≤ ≤ ≥ ≥ ⎪ ⎪ ⎪ ≤ ≤ ≤ ≥ ⎪ ⎪ ⎨ ⎪− ≤ ≤ − ≤ ≤ ⎪ ⎪ ⎪ − ≤ ≤ ≥ ≤ ⎪ ⎪ ⎩

(

)

(1.48) Next the semi-major axis a is computed as:

1 2

2 v a r GM

− ⊕

⎛ ⎞ ⎟ ⎜ ⎟ ⎜ = − ⎟ ⎜ ⎟ ⎜ ⎟ ⎝ ⎠

(1.49) where

2 2

v = v

. The eccentricity follows:

1 p a ε = −

(1.50) where the semi-latus rectum p is given by:

2

h p GM⊕ =

(1.51) and h =

× r v . The eccentric anomaly is now:

( )

2

atan2 ,1 E a n r a = ⋅ − r v

(1.52) where we have indicated

3

GM n a

⊕

=

(1.53) The mean anomaly therefore is

( ) ( ) ( )

sin M t E t E t ε = −

(1.54) with t the epochs of r and v . In order to find the remaining orbital parameter ω it is necessary to compute the argument of latitude u; this is the angle between r and the line of nodes, i.e. u

ω ν = +

:

( )

atan2 ,

z x y y x

u r r w r w = − +

(1.55) and

( )

2

atan2 1 sin ,cos E E ν ε ε = − −

(1.56) Finally

u ω ν = −

(1.57)

SLIDE 17

Dynamics_revised.doc data creazione 12/ 07/ 2001 9.35.00 Pagina 17 di 18 data ultima revisione: 29/ 04/ 2002 8.45.00

4.1.3 Latitude and Longitude

Let r be the vector representing the S/C position with respect to

E

R

, i.e. the origin of

O

R

represented in

E

R

. For nominal motion

S S

R R h

⊕

= = + r

. After some computation, the following relation holds:

( ) ( ) ( ) ( )

E

c c s s c c s s s c s s s c c s c s s c c c c s s s

i i x y i i S i z

r r R r

δ δ α δ δ α δ δ α δ δ α α Ω Ω Ω Ω Ω Ω Ω Ω

⎡ ⎤ ⎡ ⎤ + − − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = = − − + + ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎣ ⎦

R

r

(1.58) from which it is possible to compute the latitude and longitude of the S/C local vertical, as follows:

( )

latitude arcsin in degrees longitude atan2 , in degrees

z y x

r r r = = r

(1.59)

4.1.4 Small Eccentricity Orbits

For an orbit with small eccentricity ε , like the present one, the argument of the perigee ω is not well defined, since small changes in the orbit can cause large changes in the perigee location. It follows that the two reference frame vectors p and q in (1.42) are not well defined as they can vary significantly. Similarly the anomaly (true or eccentric) can vary considerably and even be non-monotonically increasing (with little physical sense): E in (1.52) can vary from

/ 2 π

to

/ 2 π −

since 1

r a −

.

In this case it is better to assume as x-axis vector for anomaly measurement the eccentricity vector e , defined as

r μ × = + v h r e

(1.60) Recalling that

= × h r v , the eccentricity vector e is well defined, since, recalling the triple product in

(1.21), results

( )

v

r μ

≈

= − ⋅ − e r v r v r

(1.61) Once e is normalized, it is possible to determine the y-axis as

× e w e 4.2 Starting Conditions for the Orbital Simulation

According to GMV and Alenia documents, it is necessary to assume the perigee of the orbit on the equatorial plane (i.e. on or near the ascending node) since in this situation the worst case for drag forces is considered. GMV documents fix the eccentricity between 1 and 2 × 10-3, but their simulations are obtained with

3

1 10 ε

−

= ×

. The inclination angle is 96.5° and the RAAN angle is

293 Ω =

assuming a launching date

n October 14, 2004.

The semi-major axis is assumed to be the equatorial earth radius plus the satellite mean altitude, i.e.

6628.1 km. a =

Since the perigee coincides with the equatorial radius, the line of nodes coincides with the direction of the perigee, and

ω =

follows.

The numerical values of the orbital parameters are summarized in Table 1.

SLIDE 18

Dynamics_revised.doc data creazione 12/ 07/ 2001 9.35.00 Pagina 18 di 18 data ultima revisione: 29/ 04/ 2002 8.45.00

parameter symbol value unit semi-major axis a 6628.1 km eccentricity

ε

3

1 10− ×

inclination i 96.5 deg RAAN

Ω

293 deg argument of the perigee

ω

deg Initial epoch t0 TBD s Table 1: numerical values of the orbital parameters The numerical values of the S/C mass and inertia parameters are summarized in Table 2. parameter value unit mass 920 kg Jxx 117 kgm2 Jyy 1357 kgm2 Jzz 1330 kgm2 Table 2: numerical values of the S/C mass and inertia parameters With these parameters, the initial position and velocity vector components are:

2587.2 = 6095.1 km 0.80890 0.34336 km/s 7.7127 ⎡ ⎤ ⎢ ⎥ ⎢ ⎥ − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ ⎡ ⎤ − ⎢ ⎥ ⎢ ⎥ = − ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣ ⎦ r v

(1.62) If the RAAN angle is assumed to be

Ω =

the initial conditions are