Application of a satellite aerodynamics model based on normal and tangential momentum accommodation coefficients

Planer. Space Sci.. Vol. 39, No. IO, pp. 1405-1419. Printed in Great Britain.

1991 0

00324633/91 $3.00+0.00 1991 Pergamon Press plc

APPLICATION OF A SATELLITE AERODYNAMICS MODEL BASED ON NORMAL AND TANGENTIAL MOMENTUM ACCOMMODATION COEFFICIENTS P. MOORE and A. SOWTER Department of Computer Science and Applied Aston University, Aston Triangle, Birmingham

Mathematics, B4 7ET, U.K.

(Received 27 May 1991) Abstract-Tangential and normal momentum accommodation coefficients are reviewed as a parametrization of the gas-surface interaction for modelling aerodynamic forces on artificial Earth satellites. By reference to experimental data and other studies, a simple empirical relationship for the normal momentum coefficients, u’, is deduced whilst the tangential coefficient, (r, is taken to be invariant with the angle, 5, between the incident stream and the surface normal. Parameters within these coefficients are estimated by least-squares reduction of orbital elements of ANS-1 (1974-704). Derived results are r~ = 0.93 and u’ = 1.114). 17 set r, yielding a drag coefficient for a sphere of 2.05kO.03, which is less than 10% below the accepted values of 2.2 for heights near 300 km.

1. INTRODUCTION

Modelling of aerodynamics forces on low Earth satellites is dependent on knowledge of the gassurface interaction for estimation of drag and lift coefficients. Most studies utilize Schamberg’s model for the reflected gas distribution (Schamberg, 1959), although it is conceptually simpler and more direct to parametrize the forces by use of the normal and tangential momentum accommodation coefficients. In this study aerodynamic modelling is reviewed through a critical analysis of Schamberg’s model and the momentum accommodation coefficient approach, by reference to experimental results and other theoretical studies. It will be seen that the use of momentum accommodation coefficients has distinct advantages, particularly as the dependence on the angle of incidence can be deduced, at least in a qualitative sense. Furthermore, an attempt has been made to describe the aerodynamics of a satellite in a way that agrees with the accepted conventions of rarefied gas dynamics and surface science. Utilizing a relationship incorporating the angle of incidence, the momentum accommodation coefficients are estimated by a least-squares analysis of orbital elements of the Astronomical Netherlands Satellite (ANS-1) (1974-704). Previous studies of 197470A have been undertaken by Sehnal(1982a,b), based on Schamberg’s model, whilst Crowther and Stark (1991), in a parallel but independent study, utilize modified momentum coefficients as opposed to the classical definition employed within this study.

2. AERODYNAMIC

MODELLING

In order to model the aerodynamics of a satellite, it is important to initially identify the flow regime about a body. The dimensionless Knudsen number, fi, may be used for this purpose where, if h is the mean free path of the atmosphere and I is a typical dimension of the body :

hf. For Kn cc 1 standard continuum flow is encountered, whereas for Kn > 10 the gas is no longer a continuum but may be modelled as a highly rarefied and collisionless collection of molecules ; this type of flow is called free molecule and is the most common type for Earth satellites. An alternative definition of the flow regimes was suggested by Cook (1966), who, recognizing that gas particles reflected from the body of the satellite may hamper the flow, suggested that free molecular flow will occur if KC>;, T

where Vi and V, are the average speeds of particles incident to and reflected from the surface, respectively. For a rarefied gas in a Cartesian coordinate system 0x,xzx3, where the molecular velocity is given by :

c = (Cl, c2, 1405

CA

1406

P. Mooas and A.

the number of molecules, dN, with velocities (c,,c,,~~)and(c,+dc,,c,+dc,,c,+dc,),isgivenby:

between

SOWTER

and mm

dN = F(c) dc, dc, dcj,

ss0

where F is the distribution function of the gas. For free molecular flow through a uniform gas in equilibrium, Fis a Maxwellian distribution function of the form (Shidlovskiy, 1967) : F = n(27~RT)~” exp { -(c-V)~/~RT},

exp { -[(c,

- Vcos 5)’

+ (c2 + Vsin [)‘+c:]/2RT}. Thus, if we consider that each gas particle of m, the total momentum striking a flat unit area in the normal (along x,) and (along x2) directions per unit time, pi and tively, are given by :

has a mass surface of tangential zi, respec-

mc:F(c) dc, de2 dc,

Pi =

MT

=

2” cos *e-

s2c0s’c +J&+s’

co2 5)

x[l+erf(scosl)]}

-rn s -03 = pVsin

mc,c,F(c)

dc, dc, dc,

5 x cos <[ 1 + erf (s cos c)]},

(2)

where

where n is the number density of gas, Tits temperature, V the velocity of the gas relative to the satellite (the macroscopic velocity) and R the universal gas constant. If we consider the coordinate system of Fig. 1, the Maxwellian distribution function is given by : F = n(2nRT)“*

m

Zi =

(1)

erf (x) = 1Jo

6 e-‘* dt, s

p is the air density and s is the speed ratio given by :

For most satellite large due to the under the limit above equations

applications, s is assumed to be very high velocity of the platform. Thus, of s + cc (hyperthermal flow), the simplify to :

and

Strictly speaking, the limit should be that s cos 5 --) co (Cercignani, 1975) as, for grazing angles of incidence, the thermal velocity of the gas is still important. However, we will consider here that the regions for which s cos 5 << 1 are small. As the force on the body is defined as the change of the momentum per unit time, all that remains is for us to find the momenta of the reflected gas particles emitted per unit time, the difference with the incident momenta giving the force. Thus, if the momenta of the gas particles emitted per unit time are pr in the normal direction (this time pointing away from the surface) and 2, in the tangential direction, the normal momentum transfer and the tangential momentum transfer at the surface per unit time are given by : P=Pi+Pr

(3)

z = ti -z,,

(4)

and

FIG. 1. COORDINATE SYSTEM FOR A GAS PARTICLE STRIKING ANELEMENTALSURFACE.

respectively. The major problem in satellite aerodynamics is the adequate modelling of pr and z, in a way that will characterize the forces of the body. Historically this has been achieved by the modelling of the reflected distribution of the gas, the most popular of which is based on Schamberg (1959) and further enhanced by

Satellite aerodynamics model Cook (1965, 1966). We will briefly describe this model before introducing a model based on momentum accommodation coefficients.

3.SCHAMBERG'S

MODEL

REFLECTED

incidence,

1407 v, such that :

Bi, by the parameter

cos 8, = (cos f3,)“. The thermal energy accommodation coefficient, LX, describes the change in the kinetic energy of the gas stream before and after reflection, and is given by :

OFTHE

DISTRIBUTION

4mV.2-4mV_2 If we consider a body in hyperthermal free molecular flow, we can assume that all gas particles strike the surface of the body with the same velocity, Vi, at an angle of Oi with the surface. After striking the surface, Schamberg (1959) proposed that the gas particles are re-emitted from the surface with a uniform velocity, V,, in a conical beam of half-angle width $O about an axis that makes an angle of 0, with the surface (see Fig. 2). The angular distribution of reflected particles obeys Knudsen’s cosine law, where the number of particles, dN,, emitted per unit time between angles 4 and 4 + d4 is given by :

K may be derived from the assumption that the total number of particles emitted equals the number of particles, N,, incident on the surface per unit time, i.e. dN, =

dN,, s

where the integration is over the re-emitted Furthermore, N, is given by :

beam.

and E, = ZpRT, is the mean kinetic energy of a surface particle and T, is the surface temperature of the satellite exposed to the atmosphere. As the gas particles strike the surface with high velocity, we will assume that the ratio :

EW +rnVi2 is numerically approximation

resulting

negligible

such that we may use the

:

in the useful expression

:

;=Jln. We may now continue to find the reflected momenta of the gas particles in the normal and tangential directions, giving :

Ni = ni V, sin Bi, where II: is the number densitv of the incident beam. For the case of a conical three-dimensional beam :

and Z, =

pv,Zsin

v, cog e, eicos ei@(&JKcose.

I

1- (7-@w

K,!!

where @(4,,) is the beam width function conical beam, by :

27~ 1- (7c/2~Jz sin 4O The angle of reflection,

O,, is related

to the angle of

1- (2&/@ @(‘O) = l-4(24,/+

Axis of ... re-emitted

recognizing

given, for a

.i sin 24~~-Wolf) sin $o-

(240/71) '

that :

0

@ ;

=:

and CD(O)= 1. Thus, the normal surface are :

and tangential

p = pv;’ sin2 ei FIG.~. SCHAMBERG'SGAS-SURFACE INTERACTIONSCHEME.

forces on our unit

Is I1 [1+@(c$~)+

1408

P. Mooas and A.

and r=p-_2sin0icos0i

l--@($J~~

1 .I.

V, cos 8,

If we assume c(, $,, and v do not vary over a changing incidence angle, the drag coefficient for a sphere, CFHERE, is given by :

SOWTER

p = pVz sin’ and

z = pV;2 sin Bi cos

cSPHERE D

[l -JLZ],

=2

model gives :

p=

and

s

ei

giving

The quasi-diffuse

pv;‘sin2ei[i+fJiY]

and

I

l,(v) =

ei [i +&ii]

x&l

-x”‘)(l

-x2)

dx.

r = Vi* sin

0

Note that this is a very complex model for drag, involving the three coefficients c(, $,, and v, of which only t( has been adequately defined and experimentally observed. In order to use these expressions we must perform some further simplification of the problem. This is usually achieved by defining the special cases of specular and diffuse reflection. Specular reflection occurs when the characteristics of the gas particle are unchanged after reflection, except that its normal velocity component has reversed (i.e. V, = Vi and 6, = 0,). In this case :

cc=0 v=l & = 0. For diffuse reflection, the gas particle is completely absorbed by the surface, only to be reemitted at some later date, in some random direction (with mean about the surface normal), in thermal equilibrium with the surface of the body. Here, u=l

r&

=

;.

A common method for simplifying the analysis would be to measure the force for purely specular reflection and again for purely diffuse reflection, assuming that all possible configurations fall between these two extremes. However, in this case both “extremes” result in C2PHERE= 2. We could attempt a somewhat artificial generalization by allowing 40 and v to take their specular and diffuse values, resulting in quasispecular and a quasi-diffuse models. For the quasispecular model (v = 1 and &, = 0), we find that :

ei cos 8.19

resulting in CgpHERE= 2(1+&/i=& The quasi-diffuse model has remained very popular as it allows non-diffuse values of tl to be adopted in an easy analytical form. Indeed, this model formed the basis of the papers of Cook (1965, 1966), which have become standard references for satellite drag coefficients (e.g. Cook’s results have been used in the generation of atmospheric models such as CIRA, 1972, produced by the COSPAR Working Group in 1972). The main criticism levelled at Schamberg’s model and its interpretation lies with the definition and use of the specular and diffuse scattering models, particularly when it is assumed that they represent extremes of the reaction. For instance, most analyses based on Schamberg assume that, if Or = Bi for specular reflection and 6, = n/2 for diffuse reflection, then tii < 0, < n/2 for all possible cases. This is explicit in the common assumption that v > 1, but is disputed by Goodman (1969), who claims that for the majority of such cases this is incorrect, indeed, in a later paper (Goodman, 1971) he points out that v c: 1 is a characteristic of the reflection at very high velocities, citing many experimental observations (e.g. Calia and Oman, 1970 ; Romney and Anderson, 1969; Miller and Subbarao, 1970). Goodman’s results are further enhanced by observations made by Nocilla (1963). In its most general form, Schamberg’s model is a very useful and adaptive model of the gas-surface interaction at satellite velocities. However, if we require numerical results, we find that we have to prejudice the model to the constraints of specular and diffuse reflection which, although intuitive, do not appear to adequately describe the reaction at the satellite surface. If we are to look objectively at the problem of satellite aerodynamics, we would perhaps be wise to take a step back and attempt to measure the inter-

1409

Satellite aerodynamics model action using another, perhaps simpler, model based on measurable parameters, preferably a model that parametrizes the force of the body without the extra effort of modelling the scattered distribution, if possible. Such a model exists, is standard in many textbooks on rarefied gas dynamics and is based on momentum accommodation coefficients.

4. MOMENTUM

ACCOMMODATION

COEFFICIENTS

We may model the momentum transfer at the interface between a gas and a surface by the normal and tangential momentum accommodation coefficients, cr’ and cr,respectively, which are defined, using the earlier notation, as :

p = pV2(2-a’)

co? {

and r=pV*asin
5

(5)

and FL = -ApV2(2-d-o)cos*

lsin t.

(6)

Similarly, if we define the drag and lift coefficients of the plate, Cn and CL respectively, by : FD = +pV*A’C,

oI _ Pi-Pr Pi -Pw

gives us :

the hyperthermal approximation

and

and

FL = - ipV*A’C,,

2,-r, Zi_Z,’ wherep, and r, are the mean momentum components of the surface particles in the normal and tangential directions, respectively. Explicitly, pw = bpJ_

and r, = 0. Thus, we can find the reflected momenta components, these values into equations (3) and (4), using equations (1) and (2), to find (Schaaf and Chambre, 1958) :

where A’ = A cos 5 is the cross-sectional we see that : CL)= 2[o+(2-a’-(r)

area, then

cos2 r]

and CL =2(2-a’-a)cos
pr and r,, and substitute

However, before we can make any assumptions we must first examine the experimental and simulated results for the behaviour of the momentum accommodation coefficients with angle of incidence.

+f&rcost][l+erf(rcost)]j and apV’sin5 r=

2&s

_2 {e ’ ws25+&s cos {[l +erf(s cos c)]}.

If we further assume that : & N - 0, Pi

FIG. 3. DRAG AND

LIFTDIRECTIONSONA

FLATPLATE.

1410

P. MOORE and A. 5. THE DEPENDENCE

OF THE MOMENTUM

ACCOMMODATION ON ANGLE

COEFFICIENTS

OF INCIDENCE

Experimental observations of Q’ and 0 have remained elusive for many years ; when they were discussed by Schaaf and Chambri in 1958, results obtained by Millikan (1923) were used as there were no contemporary observations available. Since then, an improvement in experimental technique has led to a substantial increase in the observation and measurement of the angular dependence of normal and tangential momentum transfer at a surface of the type experienced by a satellite. Tangential momentum transfer has been perhaps the easiest of the coefficients to measure as it accounts for a frictional effect on the surface which could, for example, be measured by observing the torque on a rotating disk. Variance with angle of incidence has been measured by many authors (e.g. Thomas and Lord, 1974 ; Lord, 1977) showing much inconsistency in the observed trends. Indeed, although Doughty and Schaetzle (1969) show a marked decrease with increasing (i.e. away from the normal) angle of incidence, this contrasts with the results of Steinheil et al. (1976) where only slight variance is shown, depending on the particular gas-surface combination. This, along with the results of Liu et al. (1979) and Seidl and Steinheil(1974) seem to conclude that, in general, a constant value of e will suffice with only a small margin of error. The behaviour of the normal momentum accommodation coefficient with angle of incidence has shown much more consistency in its behaviour from experiment to experiment, decreasing markedly with increasing incident angle. This is demonstrated in the results of Seidl and Steinheil (1974); Knechtel and Pitts (1969) ; Doughty and Schaetzel(l969) ; Liu et al. (1979) ; and Moskal (1971). Under the assumption that all of the surfaces used in the experiments are “cold” (i.e. pw = 0 ),Knuth (1980) has suggested the following empirical relationship : c-J’= 1

forO<5<30

0’ = 1- &

for 30” < 5 < 90”.

Empirically, this does seem to describe behaviour of u’, although we would have culty in justifying such a marked change 5 = 30”. Many attempts to model the behaviour surface interaction by the simulation of lisions amongst variously shaped elastic

the general some diffiin value at of the gassimple colmolecules

SOWTm

have generally failed to match the results obtained in experiment, although qualitatively displaying many of the trends. A computer simulation of the problem using elastic spheres was tried by Goodman (1967) who suggested that 0’ was suitably described by the relationship : 0 = c&--o’, set 5, where ai and a’, are constants. This seems to correspond very well with experimental behaviour. Indeed, an inspection of Liu et al. (1979) suggests that a’,, x 0.84, a’, z 0.12 and cr w 1, giving CrHERE = 2.32, whereas the paper of Seidl and Steinheil (1974) suggests that ob z 1.09, a’, FZ0.35 and c z 0.72, giving CgpHERE= 2.1. A model of satellite aerodynamics based on the equations : f~’ = ah-o’,

set <

(7)

and 0 = constant

(8)

has many advantages: firstly, it represents a simple model based on measurable parameters ; secondly, we do not need to prejudice the model with definitions of specular and diffuse reflection in order to simplify it and finally, we can better relate these parameters to the work of rarefied gas dynamicists and surface scientists. In order to test the adequacy of such a model, we need to test whether it may be applied to the orbit of a real satellite with the generation of plausible results. The remainder of this paper will be concerned with modelling the perturbations of an orbit due to these parameters and applying these to the orbits of the satellite ANS- 1 (1974-704).

6. THE AERODYNAMIC OF FIXED ORIENTATION

FORCE ON A FLAT PLATE IN AN OBLATE

ATMOSPHERE

In order to resolve the aerodynamic force on a satellite, it is analytically more convenient to work in a satellite-based coordinate system Oxyz where the yaxis lies along the velocity vector of the orbiting body, the x-axis along the outward normal to the orbit and the z-axis along the normal to the orbital plane. Thus, a unit vector in the direction of velocity is given by :

y= iI 0

1 .

0

If we consider that the satellite is a flat plate of unit normal given by :

Satellite aerodynamics model

11

F, = -pVZA,[o’,

n,

n=

n*

1411 cos5+(2-a~-~)~0~~

F2 = --pV*A,[o

3

cos [+o’,

co?

(11)

co? 51,

(12)

5

+(2-06-o)

3

then a unit vector in the lift direction (i.e. in the plane of n and y perpendicular to y) is given by :


and

l=YxwY)

F, = -pV2A,[a’,

cos <+(2-&-o)

cos’ (In,.

Inx YI

(13)

[I

These force components may then be susbstituted into Lagrange’s planetary equations to find the perturbation upon the semi-major axis, a, the eccentricity, e, and the inclination, i, of the orbit in a form given as (Smart, 1953) :

111

;3 .

=*

Now,

n* = n-y = cos

which results in the condition

(9)

da _=_ dE

that

n, > 0

5+(2-a’-0)cos3

where A, is the ratio of the plate turned towards the mass of the satellite. and substituting for n2 that :

F, = -pV*A,[o

the total area of the side of the incident atmosphere and Using equations (7) and (8) from equation (9), we find

51 1 0

and

cos 5+(2-oh-a)cos*

51

force on the satellite

.F, F2 .F3

(15)

{cos w cos E-w

sin w sin E -ecos

w]F,,

(16)

respectively, where r is the distance of the satellite from the Earth’s centre, E the eccentric anomaly, o the argument of perigee and

I 3

where G is the gravitational constant and ME is the mass of the Earth. We will consider here that the normal to the flat plate, II, points towards some fixed position on the celestial sphere throughout the period of a single orbit, making a fixed angle 0 with the normal to the orbital plane. The projection n’ of the normal to the plate on the xy plane makes an angle 4 with the x-axis. Thus,

[I

+(2-a;-a)cos3

F=FD+F,=

E

P=GMr,

0

where

+F,msin

di ar _=___ dE pJi3

cos* 5 sin 51,

Thus, the total aerodynamic given by :

(14)

1

E-e)

e+f(cos

cos S+rr; co? 5

F, = -~V*.4,[o’,

fF* J I*

{]y

and F, = -pV2A,(2-a’-~)

P

(10)

must be satisfied for the side of the plate that has normal, n, to be exposed to the incident air stream. Using equations (5) and (6), the drag and lift forces in the satellite coordinate system are given by : FD = --pV2A,[ocos

2ra2V

is

L

cos e

_I

From Fig. 4 we see that : 4

=fo-.I-+$,

where ti is the angle of climb, given by (King-Hele, 1987) :

1412

P. M~~REand A. Sowran Hence, the normal to the side of the flat plate facing the incident atmosphere is given by :

where &= -1 FIG.4. THEFLAT

PLATE NORMAL PLANE

PROJECTED ONTO THE ORBITAL

OF THJZ SAllZLLZTB.

forf, < E
1

&=

elsewhere.

Thus, equations (14~(16) may be rewritten using equations (1 I)-( 13) up to first-order in e as : da z = -2paZA,(1

+2e cos E)

X [u; sin2 0 sin* (Jo-E)+ecr +&(2-&--u)

and 1 sin tj = -.._C_ e sin f, V XL-a(1 -e2) . .

& = -p&,,[o’~

sin B sin (&-E) (21)

sin’ B sin3 u,--E)]

sin’ 8(2(1 i-e

cos

E)

cos

E

J

.

^_

.

x sin’ (fc - E) + sin E sin (fO -E) cos U;, - E)f .

where f is the true anomaly ot the orbit. -bus, ‘I we see that :

f&(2-&--e)sin38{2(1-i-ecosE)cosE xsin3 (fO--E)+sinEsin’

cos#d

v

(fO-E)cos(So-E))

+ 2Et~ sin a( 1+ e cos E) cos E sin 0s - E)]

(22)

di - -paA,(l+ecosE)(cos(w+E)-ecoswj dE-

and I sintP=-i;

L a(*_e2) {sin (fO-f)+e I---

Using the standard equations for Keplerian motion, which give : rsinf

x[~‘, sinOcos5sin

sinf,).

=aesinE

(jb-E)+e(2-cr~-o)

x sin’ 8 cos 8 sin’ (SO-E)] .

In an oblate atmosphere, the density at a position in the orbit is given by (King-Hele, 1987) : p = Q,

r COSf = a(cos E-e)

ezcosE[I+ccos

2(o+E)

- 2ce sin 2(0 + E) sin E

and

+ic2{1 -l-COS~(CU+E))+O(CK, c2e,e2)J, I&-_.-!!a(t-e’)

l+ecos E 1-ecosE’

we can show that : sin4 = sin(fo-E)+O(e2)

1rK c=-_Lsin2i 2H

O(G).

Thus, the condition shown in equation (10) becomes, to order e’, sin0sin(&-E)>O, which, since 0 < 9 < rr, reduces to : sin (fO -E)

> 0.

(24)

where pP is the atmosphere density at the perigee position, fq = e-z-ccos20,

and cos Q,= cos (fo -E)+

(23)

and *,Qe H’ where H is the density scale height, rp is the distance of perigee from the Earth’s centre and ZC( % 0.~335) is the ellipticity of the Earth.

1413

Satellite aerodynamics model If ai is one of our orbital elements, a, e or i, we can find the change, Au,, in ai over a single orbital period, by integrating the appropriate expression from equations (21), (22) or (23) to find : Aai =

s

2nda. 2

dE

o

(1989). Here, we reproduce the special case of a spherically symmetric atmosphere (c = 0) as : Aa = -2p,k,a’A,

dE. 16eZ, -8Zo]a’, sin’ 0-- ~[{8eS2+8S’,}

The form of the integrand has many variables, each of which will vary with E over a single orbit. However, we can simplify this by noting that the change in all parameters, excepting E, over a single period is very slight. Therefore, we can integrate the terms that are explicitly functions of E, holding the other parameters constant. The values given to the remaining elements will be their average over the orbiting period. Thus, if we consider that one of Lagrange’s planetary equations is given as : da-’

_

posE

a, + f (a, cos nE+ b, sin nE) “=I

dE

+a

co+

f (c, cos nE+d, n=,

xcosfo+{8e(C;-C&)+8C;}

+{8e(C;+C;)+8C;} x cosfo -{24e(C;+ x (2-crb-0)

II

x cos

i: a,Z,,(r)+

f

sin 3f,-{24eS,+24S,} Ck)+24C;}

- ~[{8e(Z~+2Z2+Zo)+8Z,+24Z,}

2fo - 16e(Z2 + Z,) - 32Z,]a’, sin’ 0

-{24e(S3+S,)+40&}

COSfo

+{8e(C;+2C;+Cr)+8C;+24C;}

c,[Z,(z)+C&>fo)]

- {24e(C; + C;) f40C;

f,

sinfo]

+~[{8e(S,+2S,+S,)+8S,+24S,}cos3fo

II=0

-

cos 3fo

sin3 f?+O(e’)

Ae = -ppk,aAm

sin nE)

sin’f,]a

x sin 0+ G[{8e(S,+S,)+8SS}

where the a,, b,, c, and d, are constants (i.e. not explicitly functions of E) over an orbit, then :

fi=O

cos 2fo

- i [{8e(Z, +Z,)+8Z2}

4SMo))

sin 3fo

+ 56Ch) sinfo]

x(2-crb-a)sin’0-i[{4e(S,+S,)+8S2}

3

x cosfo + {4e(C; + 3C;) + SC; + SC;} sinfo]

where *2n

I .O

cosnEezcoSEdE=2~Z,,(z),

n=0,1,2,...

sin nE ezcosE dE = 0,

n = 1,2,3,...

x fl sin 0 + 0(e2)

2n

JO

and , Ai = -ppklaA,

[Z,(z) is the modified Bessel function and imaginary argument],

s s

- i[{4e(Z3-Z,)+8Zz}

sin@-fo)

of the first kind - 81, sin (w +fo)]a’, sin 0 cos 0

fO+x

GWO) = i n

and

cos nE ezcasE dE

n = 0, 1,2,

.

fo

-2S,

fO+”

sz(z,fo>= f

sin nE ercosE dE

3L

+ -[{4e(S,-S,)+ 16

n = 1,2,3,.

.

71 fo

The expressions that this creates for the elements a, e and i in an oblate atmosphere are large and the computer algebra system REDUCE was used in their derivation. The full expressions are contained in Sowter

12S3} sin (w-2f,)

sin 5w+{4eS2-8S,}

sin (w+2fo)

-{8eS2-16S,}sinw+{8C;-4eC;) x cos(m-22fo)+{4e(Cb-C;)+8C’,) cos(o+2f,)+{8e(C~--C~)+16C~} x (2-ai-a)

sin’ Bcos e+0(e2)

cos 01

P. Mooan and A.

1414 where Ifi = L(z)

n=0,1,2,...

S, = S?&,fo)

n=

c: = C,(z,fo)-Z”(Z)

n = 0,1,2,.

1,2,3,...

and

‘7.THE

ANS-1 SATELLITE

ma1 for the Sun-orientated plates is described explicitly by the angles B,, = B. and f,, =fo, as follows. In the Earth-centred system based on the equation and equinox of date, the Sun direction is specified by the right ascension, coo, and declination, 6,. Rotating this system through angles R, i and o, the latter identified with E = 0,the required angles are defined uniquely by :

cos O. = sin 6, cos i-cos

(1974-70.43

The Astronomical Netherlands Satellite, ANS-1, was launched on 30 August 1974 to undertake ultraviolet and X-ray observations (Wakker, 1978). The nominal pre-launch orbit was Sun-synchronous and near-circular with perigee near 510 km. A launch malfunction caused the actual orbit to deviate, placing it in an eccentric orbit with e N 0.064 and perigee of 266 km. Despite this a large number of scientific objectives were achieved before the operational life ended on 11 December 1975. The satellite decayed on 14 June 1977. Structurally, the satellite consisted of a rectangular prism densely packed with instrumentation. Protruding from opposite sides were two solar panels of fixed orientation to the body. The satellite was able to maintain Sun orientation by three-axis stabilization, which also defined the geometry of the telescopes perpendicular to the Sun-direction, resulting in full sky coverage over one year. For the purpose of this analysis we follow Sehnal (1982a,b) and assume the simplified form of Fig. 5, where the areas of each panel are P, = 0.517 m*, P, = 0.898 m*, P3 = 0.750 m* and P, = 0.445 m*. ANS-1 had mass 129.6 kg. Faces P, and P, are hereafter assumed to be Sunorientated, which for simplicity is approximated by the Earth-Sun direction. Using equation (20), the nor-

!hWTER

sin B. cosfo

6, sin ices (a,--Cl)

= cos 6, cos (c(o -Q)

cos 0

+cos6,cosicos(a,-Q)sinw+sin6,sinisinw sin8,sinf,=cos6,cosicos(a,-R)cosw +sin6,sinicoso-cos6,cos(cr,-R)sinw. Corresponding expressions e,(@,) for 0 and f2(f4) for fin equation (20) can be derived for side panel P,(P,). Introduce a coordinate system Ox’y’z’ with Oz’ Sunpointing by rotating the conventional satellite based system Oxyz through (n/2) - 4. (r$o = f. - E) about the z-axis, followed by a rotation 8, about the new x-axis, Ox’. The normal to P, relative to Ox’y’z’ can be written :

cos p 4 =

sin /3

t

0 i

for some angle p, which defines the attitude. to the Oxyz frame yields : sin e2 cos 42 sin

n; =

(

e2sin 42 co~e2 1 sin~,cosfi+cos~, =

By direct evaluation

I FIG. 5.Tm DIMENSIONSOF ANS-1

cos8,= sin

e2sinf,

cosB,sinp

.

-cos~,cos/?+sin~,cos0,sinp

-sin

i

1.23m

Rotation

-sine,

= -cosSo

e.sin

fi

i

: sin fi cos /J+sinf,

sin B, cosf2 = sinfo cos /J+cosf,

cos B. sin fi cos

e,sin p. (26)

Similarly e4, f4 for panel P4 are defined by replacing /J by /3 - (n/2) in equation (26). It is unrealistic to consider each panel in its entirety as certain portions of the surface will be shadowed from the incident atmosphere. Given the lack of attitude information, only a coarse estimation of shadow-

Satellite aerodynamics model ing was undertaken, which manifests itself in a reduction in total area for a panel used in the calculation for aerodynamic perturbations. Sehnal (1982) showed that 8o rarely exceeded 20” which implies that the greatest shadowing is of P, by P, with a smaller amount on P2 by P,. Details of the shadowing are given by Sowter (1989) with the conclusion that the total Sun-orientated panel should be reduced by 10% whilst P2 should be reduced by 5%.

(i) Orbital perturbations of ANS-1 The ESOC orbital elements for ANS- 1 were kindly supplied by L. Sehnal of the Czechoslovak Academy of Sciences. For analysis of aerodynamic effects, other perturbations due to atmospheric rotation and gravitational attraction require consideration. Perturbations due to the zonal harmonics in the geopotential and third body attraction of the Sun and Moon were removed by numerical integration of the equations of motion, utilizing the PROD program (Cook, 1972). In addition, ANS-1 passed through 15: 1 resonance with the longitude-dependent part of the Earth’s gravity field around MJD 42800 ; the effect of this is clearly visible in the inclination. For the semimajor axis and eccentricity no such perturbation is visible and thus resonance is discussed only as applicable to the inclination. The theory of resonance (Allan, 1973) has been developed in terms of lumped harmonics, i.e. linear combination of individual harmonics. To estimate the resonance perturbation, gravity field coefficients were taken from the GEM-1OB global model (Lerch et al., 1981), whilst lumped harmonics were available from the satellite derived coefficients (King-Hele and Walker, 1987) and from ANS-1 itself (Klokocnik, 1982). Utilizing coefficients from the above sources a decrease of 0.0065”~.0083” was estimated for the inclination. Despite the broad agreement, the discrepancy between the results is large compared with the perturbation anticipated for atmospheric lift. Thus, for ANS-1, only inclination values prior to resonance were considered, with analysis restricted to the period MJD 42290-MJD 42700. Given negligible effect on a and e, all elements up to the final values on MJD 42988 were accepted. Atmospheric rotation is also pertinent for the inclination, which is particularly significant for a near polar orbit. Perturbations were evaluated assuming a constant drag coefficient over an orbit. For ANS-1 where c, the oblateness coefficient, exceeds 0.2 and e < 0.2, the required expression is given by Boulton and Swinerd (1983) with the change in the inclination, Ai, over one orbit given by :

1415 A sin i A Ai = ----AT,

6J-F B where A is the rotation rate; AT the change in the orbital period over one revolution ; A and B functions of the orbital elements, the density scale height Hand c, whilst

uA being the rotational velocity of the atmosphere and u the satellite velocity relative to the air. Secondary orbital perturbations are caused by the solid earth and ocean tides and solar radiation pressure. For the latter, a spherical satellite with the same orbit and area-to-mass ratio as ANS-1 revealed perturbations in the inclination of amplitude 10p4’ or less. This is significantly less than the assumed accuracy of the ESOC data which is around 1.5 x 10p3” for the inclination. Removal of the aforementioned perturbations yielded the values of i, a and e plotted in Figs 68, where the variation can be attributed solely to aerodynamic effects of drag and lift. Evaluation of lift and drag effects necessitates estimation of atmospheric density at perigee by utilizing one of the numerous thermospheric models available to the analyst. The lifetime of ANS-I spanned the minimum of the 1 l-year cycle of solar activity. Models such as DTM (Barlier et al., 1978) or MSIS79 (Hedin et al., 1979) were based on data from periods of moderate solar activity with extrapolation to more extreme conditions questionable. In comparison 577

-l.OOF 42.29

a Time (MJDI

x103

42.70

FIG. 6. THE ORBITAL INCLINATION OF ANS-1 CLEARED OF ZONALHARMONICS,ATMOSPHERICROTATIONANDLUNI-SOLAR ~RA~ITA~~NALEFFECTSWITH~~RE~CALFIT(CONTINUOLJS CURvE)USING 577 DENSITYMODELANDA = 1.10.

1416

P. MOORE and A. SOWTER

Values of {xi}:= ,, were first determined for each of the orthogonal plates, corresponding to the Sun-orientated plate and the two side panels with allowance for shadowing, the orientation of the side panels being defined by the angle /3. Given no prior knowledge of attitude, we simply averaged over values of /I = O”, 30”, 60”, 90”, 120” and 150”. No further values in this sequence were necessary, given the symmetry of geometry. Using averaged values of xi, i = 0, 1,2,3, the parameters cr, ah and a’, were sought by leastsquares adjustment of the observed values.

6600: 42.29 x103

42.99

Time (MJD) FIG. 7. THE SEMI-MAJOR FIT (CONTINUOUS

AXIS OF

CURVE)

ANS-1

USING

WITH

THEORETICAL

577 DENSITY MODEL.

(Jacchia, 1977) was developed from data spread over almost two complete solar cycles. J77 differs from MSIS79 and DTM for low solar activity but later MSIS models (Hedin, 1983, 1987) appear to confirm the 577 neutral air-densities. Given this confirmation we have used 577 in the subsequent analysis, with values relative to MSIS79 derived for comparison. For the elements a, e and i, expressions for the change in the element over one revolution were evaluated by use of equations (21)-(23) with integration undertaken by use of the algebraic manipulation computer package, REDUCE. The output yielded the change in an arbitrary element, x, of the form : Ax = xO+x,a+x&,+x,a’,.

-0.3ot 42.29

’

s

’

Time (MJD)

FIG. ‘8. THE HARMONICS THEORETICAL

ECCENTRICITY

AND

OF

LUNI-SQLAR

FIT (CONTINUOUS

ANS-1

CLEARED

GRAVITATIONAL CURVE)

MODEL.

( 42.99

x103

USING

OF ZONAL

EFFECTS WITH

577

DENSITY

(ii) Analysis of the orbital inclination In consequence of the analytic form of Ai, equation (23), only a’, and the linear expression (2 - 06 -a) can be estimated by least-squares adjustment. In addition a shift, Ai,, in the initial value of the inclination was sought from least-squares analysis of the inclination values between MJD 42290 and MJD 42700. Relative to A = 1.1 and the 577 atmospheric model, the derived results and their correlation coefficients are summarized in Table 1. For comparative purposes, the values relative to J77 and A = 1.15 were -0.13f0.07, 0.17f0.01 and 1.44 x IO-‘f 1.6 x 1O-4 for the three parameters. Using values from Table 1, the theoretical values of the inclination are compared against observed values in Fig. 6. The fitted curve is far from ideal, although it demonstrates the general behaviour. The discrepancy may be a consequence of a relaxation in the Sunorientation of the solar panels, which was necessary to achieve the mission objectives. Furthermore, the attitude assumption of side panel orientation, through averaging over the angle /I, may also be erroneous. Sehnal (1982a,b) recognized the difficulties in obtaining a good fit when attempting an aerodynamic analysis using Schamberg’s model, based on a quasi-constant thermal accommodation coefficient, tl. Sehnal proposed a linear increase in c( from 0.1 at MJD 42290 to 1.O at MJD 42550 to yield an improved fit. Sehnal postulated that the variation could be explained by a degradation in the surface reflection characteristics with the reflection mode becoming more diffuse with time. Although plausible, we prefer to limit ourselves to invariant momentum parameters, given the uncertainty in the attitude and the lack of any serious evidence for surface degradation and associated gassurface reflection changes over a protracted time period. The effect of thermospheric mis-modelling is not severe for the inclination, as can be seen from Table 2, where the estimation procedure utilized MSIS79 and A = 1.1 for atmospheric rotation.

1417

Satellite aerodynamics model TABLE1. VALUESFORTHEINCLINATION : 577 ANDA = 1.1

Parameter

Least-square value

2--86--a

-0.04+0.07

I

i:, (deg)

1.26x IO-3+1.6x 0.17+0.01

1O-4

-0.64 0.62

(iii) Analysis of the semi-major axis

The perturbation of the semi-major axis by aerodynamic forces is given by equation (21). In contrast to the inclination the structure of equation (21) gives the possibility of solving for cr, (~6and u’, separately, however this is illusory. The least-squares procedure for 5’77yielded the normal equations :

g

_iliiii

~;;~~~) (il)

= (_i/z‘iz)

( with the inevitable consequence of unacceptable correlations. Given the obvious similarity in the three equations linking cr, & and o’,, and in order to derive a meaningful result, a single equation was constructed by addition of equations, namely : CT--1.179&+ 1.688~; = -0.089. Combining this equation with the corresponding inclination values from Table 1 yields

c; = 1.11+0.03

-0.16 -

Ai,

-0.64 -0.16 -

o; = 1.33kO.03 a; = 0.1.5~0.01. The consequence of thermospheric mis-modelling is apparent within the right-hand sides of the two equations derived from analysis of the semi-major axis. In addition to the previous argument in favour of 577, note that at a normal incidence angle ( ti = 0), g’ > 1 for MSIS79, which is unacceptable. (iv) Analysis of the orbital eccentricity Analysis of the eccentricity parallels that of the semi-major axis, as equation (22) also leads us to understand that all coefficients may be solved for. However, in complete analogy the results are highly correlated and we again seek a representative equation from the least-squares procedure. For J77 the analysis gives : Q- 1.1840;+ 1.609a; = -0.217, with the inclination

fJ = 0.88 + 0.04

a; = 0.17f0.01. Analogously for the MSIS79 thermospheric model the single equation derived from the semi-major axis was :

which, when combined yielded :

0.62

which, when combined from Table 1, yielded :

c = 0.93 * 0.04

g-- 1.218a~+1.733~;

I 01

2%oh-0

= -0.603,

oh = 1.16kO.03 c’, = 0.17~0.01. In contrast the MSIS79 model yielded : u - 1.2300; + 1.659a; =

with values from Table 2, which with Table 2 gives : CT= 0.72kO.03 -0.714,

Q = 0.75 + 0.03

TABLE2. INCLINATION VALUES : MSIS79 ANDA = I. 1 Correlation coefficients Parameter

2-4-u J ii, (deg)

Least-square value

-0.08*0.06 1.36x 10-3f 0.15+0.01 I.5 x lo-”

-0.60 0.58

0.58 -0.05 -

-0.60 -0.05 -

values

P. Mooan and A. Sowma

1418

@; = I .36 rf:0.03 cr’,= 0.15_to.o1. The similarity of solution between the semi-major axis and eccentricity is reassuring, given that both are perturbed predominantly by the along-track aerodynamic force and thus, the two analyses are qualitatively equivalent. Graphs of the observed and theoretical values of the semi-major axis and the eccentricity are given in Figs 7 and 8. Finally, we note the insensitivity of the derived values with the atmospheric rotation rate, A. Solutions with either thermospheric model, but using A = 1.15, yielded agreement at the 20 level. The level of insensitivity is important given that the final values are heavily dependent on the inclination. (v) Drag coe~c~en~s for a sphere The preferred results for the momentum modation coefficients are : 0 = 0.93

and

(r’= 1.11-0.17

accom-

eccentricity and, even in the presence of many possible errors, the widely scattered observed values of the inclination bear reasonable resemblance to the theoretical fit. The preferred solution of: d = 0.93,

CT’= 1.11-0.17 set 5,

yields a drag coefficient of 2.05 for a sphere, a value that is within 10% of the accepted value at heights near 300 km. Otherwise the feasibility of the solution is difficult to quantify, although rr’ < 1 is consistent. Furthermore, laboratory results for 0 are generally high, often exceeding unity (a characteristic of back scatter) but seldom falling below 0.5. The derived values are thus not inconsistent but, given the level of uncertainty within the methodology, must be treated with caution until validated by similar studies.

set 5,

as derived from the inclination and semi-major axis relative to 577 and atmospheric rotation rate A = l,l. Utilizing these expressions, the drag coefficient for a spherical body is calculated to be :

which is less than 10% below the minimum value set by Cook (1966) using Schamberg’s model from which CgHERE= 2.2 is the accepted value for heights near 300 km. 8. CONCLUSIONS Within this study an attempt has been made to describe the aerodynamics of a satellite in a way that agrees with the accepted conventions of rarefied gas dynamics and surface science. Deficiencies in models utilizing the scatter distribution of the gas have been circumvented by consideration of momentum accommodation coefficients which to a first approximation obey the empirical formulae : cf = const,

As far as the orbital fit is concerned, Figs 68 are visually good fits with the theoretical values exhibiting the general behaviour of the observed elements. This is particularly true for the semi-major axis and the

CT’= &--a’, sect,

although alternative variations with the incident angle can be introduced. Momentum coefficients of the above form have been used to predict the aerodynamics of satellite motion and applied to the orbit of ANS- 1. The question to,be raised here is whether the application has been successful in terms of the degree of fit to the orbital elements and also in terms of the feasbility of the values of the parameters describing the gas-surface interaction.

REFERENCES

Allan, R. R. (1973) Satellite resonance with longitude dependent gravity-III. Inclination changes for close satellites. Planet. Space Sci. 21,205. Barlier, F., Berger, C., Falin, F. L., Kockarts, G. and Thuillier, G. (1978) Atmospheric model based on satellite drag data. Ann. Geophys. 3b, 9. Boulton. W. J. and Swinerd. G. G. (1983) The change in satellite orbital inclination due to a rotating oblate airnosphere with a diurnal variation in density. Proc. R. Sot. A 389,349. Calia, V. S. and Oman, R. A. (1970) Scattering cross-section measurements for epithermal Ar on Ag(lll) surfaces. J. them. Phys. 52,6184. Cercignani, C. (1975) Theory and Applications of the Boltzmann Equation. Scottish Academic Press, Edinburgh. CIRA (1972) COSPAR International Rejkence Atmosphere 1972. Akademie, Berlin. Cook, G. E. (1965) Satellite drag coefficients. Planet. Space Sci. 13, 929. Cook, G. E. (1966) Drag coefficients of spherical satellites. Ann. Geophys. 22,53. Cook, G. E. (1972) Basic theory for PROD, A program for computing the development of satellite orbits. Gel. Mech. 7, 301. Crowther, R. and Stark. J. (1991) The determination of the gas&face interaction from satellite orbit analysis as annlied to ANS-1 (1974-704). Planet. Soace Ski. 39.729. Dot&y, R. 0. and‘Schaetzle; W. J. (1969) Experimental dete~ination of momentum accommodation coelhcients at velocities un to and exceeding earth escape velocity. In Rarefied Gas’ Dynamics (Edited by Triiling, L. -and Wachman. H. Y.), Vol. II. Academic Press, New York. Goodman, F. 0. (1967) A three-dimensional hard spheres theory of scattering of gas atoms from a solid surface. NASA CR-933. Goodman, F. 0. (1969) Empirical representation of the vel-

Satellite aerodynamics ocity distribution density function of gas molecules scattered from a solid surface. In The Structure and Chemistry of Solid Surfaces (Edited by Somorjai, G. A.). Wiley, New York. Goodman, F. 0. (197 1) Review of the theory of the scattering of gas atoms by solid surfaces. Surf. Sci. 26, 327. Hedin, A. E. (1983) A revised thermospheric model based on mass spectrometer and incoherent scatter data: MSIS83. J. geophys. Res. 88, 10170. Hedin, A. E. (1987) MSIS-86 Thermospheric model. J. .qeoph ys. Res. 92, 4662. Hedm, A. E., Reber, C. A., Spencer, N. W. and Brinton, H. C. (1979) Global model of longitude/UT variations in thermospheric composition and- temperature based on mass spectrometer data. J. geophys. Res. 84, 1. Jacchia, L. G. (1977) Thermospheric temperature, density and composition : new model. SAO Special Rep., p. 375. King-Hele, D. G. (1987) Satellite Orbits in an Atmosphere, Theory and Applications. Blackie, Glasgow. King-Hele, D. G. and Walker, D. M. C. (1987) Geopotenfial harmonics of order 15 and 30, derived from the orbital resonances of 25 satellites. Planet. Space Sci. 35, 79. Klokocnik, J. (1982) A readjustment of lumped coefficients from inclination of 1974-704. Bull. Astron. Inst. Czech. 33, 259. Knechtel, E. D. and Pitts, W. C. (1969) Experimental momentum accommodation on metal surfaces of ions near and above earth-satellite speeds. In Rarefied Gas D_vnamics (Edited bv Trillina. L. and Wachman. H. Y.). Academic Press. New York.-’ Knuth, E. L. (1980) Free-molecule normal-momentum transfer at satellite surfaces. AZAA J. 18, 602. Lerch, F. J., Putney, B. H., Wagner, C. A. and Klosko, S. M. (1981) Goddard Earth Models for oceanographic applications (GEM 10B and IOC). Marine Geodesy 5, 145. Liu, S.-M., Sharma, P. K. and Knuth, E. L. (1979) Satellite drag coefficients calculated from measured spatial and energy contributions of reflected helium atoms. AZAA J. 17, 1314. Lord, R. G. (1977) Tangential momentum accommodation coefficients of rare gases on polycrystalline metal surfaces. AZAA Progress in Astronautics and Aeronautics, Rarefied Gas Dynamics, Vol. 51, Pt 1 (Edited by Potter, J. L.). New York, pp. 531-538. Miller, D. R. and Subbarao, R. B. (1970) Scattering of 0.06 2.5-eV neon and argon atoms from a silver( 111) crystal. J. them. Phys. 51,425. Millikan, R. A. (1923) Coefficients of slip in gases and the

model

1419

law of reflections of molecules from the surfaces of solids and liquids. Phys. Rev. 21,217. Moskal, E. J. (1971) An investigation of the normal momentum transfer for gases on tungsten. Univ. of Toronto Institute for Aerospace Studies, UTIAS Rep., p. 166. Nocilla, S. (1963) The surface re-emission law in free molecule flow. In Rarefied Gas Dynamics, Vol. VI, p. 327. Academic Press, New York. Romney, M. J. and Anderson, J. B. (1969) Scattering of 0.05-5-eV argon from the (111) plane of silver. J. them. Phys. 51,249O. Schaaf, S. A. and Chambre, L. (1958) The flow of rarefied gases. Fundamentals of Gas Dynamics (Edited by Emmons, H. W.). Princeton University Press, Princeton, New Jersey. Schamberg, R. (1959) A new analytic representation of surface interaction for hyperthermal free molecule flow with application to neutral-particle drag estimates of satellites. Rand Corporation Reports, RM-23 13. Sehnal, L. (i982a) Change of the orbital inclination of satellite 1974-704. Bull Astron. Inst. Czech. 33. 16. Sehnal, L. (1982b) Determination of satellite drag coefficient from the orbital analysis of the ANS satellite (1974-704). Bull. Astron. Inst. Czech. 33,244. Seidl, M. and Steinheil, E. (1974) Measurement of momentum accommodation coefficients on surfaces characterized by auger spectroscopy, sims and leed. In Rarefied Gas Dynamics (Edited by Becker, M. and Fieberg, M.), Vol. II. DFVLR-Press, Part-Wahn, Germany. Shidlovskiy, V. P. (1976) Introduction to the Dynamics of RareJied Gases. American Elsevier, New York. Smart, W. M. (1953) Celestrial Mechanics. Longmans, London. Sowter, A. (1989) Drag coefficients with application to satellite orbits. Ph.D.Thesis, Aston University. Steinheil, E., Scherber, W., Seidl, M. and Rieger, H. (1976) Investigations on the interaction of gases and well-defined solid surfaces with respect to possibilities for reduction of aerodynamic friction and aerothermal heating. Proc. Tenth Znt., Symp. on Rarefied Gas Dynamics, Aspen, Colorado, p. 589. Thomas, L. B. and Lord, R. G. (1974) Comparative measurements of tangential momentum and thermal accommodations on polished and on roughened steel spheres. In RareJied Gas Dynamics (Edited by Karamcheti, K.), pp. 405412. Academic Press, New York. Wakker, K. F. (1978) Orbit prediction for the Astronomical Netherlands Satellite. J. Br. Znterplanet. Sot. 31, 387.