> restart:#"m13_p35"

Consider a cubic black-body satellite, of 0.5 m in size, three-axis stabilised in an equatorial orbit at 300 km (LEO), and the case it is at GEO. Find:

a) Orbital period.

b) Eclipse duration.

c) Solar input along the orbit.

d) Infrared input from the planet, assumed at a temperature of 288 K and with =0.6.

e) Albedo input along the orbit, assuming an albedo of 0.3 and a simple albedo model.

f) Periodic temperature evolution, assuming the body is isothermal, with a mass of 50 kg and a thermal capacity of 1000 J/kg•K).


> read`../therm_eq.m`:read`../therm_const.m`:read`../therm_proc.m`:with(therm_proc):with(plots):

> dat:=[alpha=1,L=0.5*m_,R=6370e3*m_,E=1370*W_/m_^2,H_LEO=300e3*m_,H_GEO=35.8e6*m_,G=6.67e-11*N_*m_^2/kg_^2,M=5.97e24*kg_,Tp=288*K_,epsilon=0.6,rho=0.3,m=50*kg_,c=1000*J_/(kg_*K_)];

[alpha = 1, L = `+`(`*`(.5, `*`(m_))), R = `+`(`*`(0.6370e7, `*`(m_))), E = `+`(`/`(`*`(1370, `*`(W_)), `*`(`^`(m_, 2)))), H_LEO = `+`(`*`(0.300e6, `*`(m_))), H_GEO = `+`(`*`(0.358e8, `*`(m_))), G = `...
[alpha = 1, L = `+`(`*`(.5, `*`(m_))), R = `+`(`*`(0.6370e7, `*`(m_))), E = `+`(`/`(`*`(1370, `*`(W_)), `*`(`^`(m_, 2)))), H_LEO = `+`(`*`(0.300e6, `*`(m_))), H_GEO = `+`(`*`(0.358e8, `*`(m_))), G = `...



Fig. 1. a) Lateral view of equatorial orbits. b) Orbit view from the North Pole, to show the angular potion origin along the orbit, and the eclipse period (in black). c) Cube face numbering.

Procedure p(phi) wrap around phi to the 0..2*Pi interval.

> dat:=op(dat),Const,SI2,SI1:p:=proc(phi) evalf(2*Pi*frac(phi/(2*Pi)));end;

proc (phi) evalf(`+`(`*`(2, `*`(Pi, `*`(frac(`+`(`/`(`*`(`/`(1, 2), `*`(phi)), `*`(Pi))))))))) end proc

a) Orbital period.

From orbital mechannics:

> Orbit_period:=T=2*Pi*sqrt(a^3/(G*M));a=R+H;T_LEO:=subs(SI2,evalf(subs(a=R+H_LEO,dat,Orbit_period)));T_GEO:=subs(SI2,evalf(subs(a=R+H_GEO,dat,Orbit_period)));

T = `+`(`*`(2, `*`(Pi, `*`(`^`(`/`(`*`(`^`(a, 3)), `*`(G, `*`(M))), `/`(1, 2))))))
a = `+`(R, H)
T = `+`(`*`(5423.985711, `*`(`^`(`*`(`^`(s_, 2)), `/`(1, 2)))))
T = `+`(`*`(86225.39395, `*`(`^`(`*`(`^`(s_, 2)), `/`(1, 2)))))

i.e. the LEO period is 5400 s (90 min), and GEO period 86200 s (23 h 58 min). Notice that the influence of H in LEO is not great, and that the GEO altitude can be computed from orbital mechanics knowing that its period is a sideral day (roughly 24 h).

> plot(subs(a=R+1000*H_km,dat,SI0,rhs(Orbit_period)),H_km=0..1000,T_s=0..6500);


Notice that we do not need to find G and M, since the GM-product can be found from g and R by the relations:

> eqG:=F=G*M*m/R^2;eqg:=F=m*g;eqGg:=g=G*M/R^2;eqGg_:=subs(dat,%);

F = `/`(`*`(G, `*`(M, `*`(m))), `*`(`^`(R, 2)))
F = `*`(m, `*`(g))
g = `/`(`*`(G, `*`(M)), `*`(`^`(R, 2)))
`+`(`/`(`*`(9.80665, `*`(m_)), `*`(`^`(s_, 2)))) = `+`(`/`(`*`(9.813440652, `*`(m_)), `*`(`^`(s_, 2))))

b) Eclipse duration.

The eclipse duration depends on the inclination of the equatorial plane to the ecliptic. Maximum duration occurs when the Sun is in the equatorial plane, i.e. at the equinocces.

We only consider one orbit (phi=0..360º from subsolar point), and only this maximum eclipse duration (with Fig. 1).

> Eclipse_start:=phi[es]=Pi/2+arccos(R/(R+H));Eclipse_start:=phi[es]=Pi-arcsin(R/(R+H));Eclipse_fraction:=T[e]/T[o]=expand(1-rhs(Eclipse_start)/Pi);eqh:=h=(H+R)/R;LEOh:=subs(H=H_LEO,dat,eqh);GEOh:=subs(H=H_GEO,dat,eqh);LEO_es:=phi_es=evalf(subs(H=H_LEO,dat,rhs(Eclipse_start)));LEO_es_deg:=phi_es=evalf(rhs(%)*180/Pi)*'º';LEO_ee:=phi_ee=evalf(2*Pi-rhs(LEO_es));LEO_ee_deg:=phi_ee=evalf(360*º-rhs(LEO_es_deg));LEO_ef:=evalf(subs(H=H_LEO,dat,Eclipse_fraction));LEO_ed:=rhs(%)*rhs(T_LEO);GEO_es:=phi_es=evalf(subs(H=H_GEO,dat,rhs(Eclipse_start)));GEO_es_deg:=phi_es=evalf(rhs(%)*180/Pi)*'º';GEO_ee:=phi_ee=evalf(2*Pi-rhs(GEO_es));GEO_ee_deg:=phi_ee=evalf(360*º-rhs(GEO_es_deg));GEO_ef:=evalf(subs(H=H_GEO,dat,Eclipse_fraction));GEO_ed:=rhs(%)*rhs(T_GEO);

phi[es] = `+`(`*`(`/`(1, 2), `*`(Pi)), arccos(`/`(`*`(R), `*`(`+`(R, H)))))
phi[es] = `+`(Pi, `-`(arcsin(`/`(`*`(R), `*`(`+`(R, H))))))
`/`(`*`(T[e]), `*`(T[o])) = `/`(`*`(arcsin(`/`(`*`(R), `*`(`+`(R, H))))), `*`(Pi))
h = `/`(`*`(`+`(R, H)), `*`(R))
h = 1.047095762
h = 6.620094193
phi_es = 1.871857043
phi_es = `+`(`*`(107.2495084, `*`(`?)))
phi_ee = 4.411328265
phi_ee = `+`(`*`(252.7504916, `*`(`?)))
`/`(`*`(T[e]), `*`(T[o])) = .4041693977
`+`(`*`(2192.209038, `*`(`^`(`*`(`^`(s_, 2)), `/`(1, 2)))))
phi_es = 2.989956966
phi_es = `+`(`*`(171.3119150, `*`(`?)))
phi_ee = 3.293228342
phi_ee = `+`(`*`(188.6880850, `*`(`?)))
`/`(`*`(T[e]), `*`(T[o])) = 0.4826713848e-1
`+`(`*`(4161.853030, `*`(`^`(`*`(`^`(s_, 2)), `/`(1, 2)))))

i.e. the LEO eclipse, which starts 107º after the subsolar point, lasts 40% of the orbit (2200 s, 37 min), whereas the GEO eclipse, which starts 171º after the subsolar point, only lasts 5% of the orbit, but this amounts to 4160 s (69 min, 72 minutes

including the penumbra, nearly double than in LEO).

Notice that the above computations refer only to maximum values, which occur at 20 March and 22 September. LEO eclipse duration varies little with the season (i.e. with solar declination), but in GEO the influence is great: there are only GEO eclipses near the equinocces (when solar declination <8.7º, Fig. 1, which occurs from 27-Feb to 11-Apr and from 28-Aug to 11-Oct).

c) Solar input along the orbit.

The energy absorbed (or 'heat inpout') varies because the frontal area changes, and is zero in eclipse.

It is better to build a sunlit model for each face and then all together (we are forced to build separate models in the case each face has a different temperature).

> Qs0:=alpha*A*E;eqAface:=A=L^2;Qs0_:=subs(eqAface,dat,Qs0);Qs1:=Qs0*piecewise((p(phi)<Pi/2 or p(phi)>3*Pi/2),cos(phi),0);Qs2:=Qs0*piecewise(p(phi)>(2*Pi-phi_es),-sin(phi),0);Qs3:=Qs0*piecewise(p(phi)<phi_es,sin(phi),0);Qs4:=0;Qs5:=0;Qs6:=0;Qs:='Qs1+Qs2+Qs3+Qs4+Qs5+Qs6';plot(subs(LEO_es,eqAface,dat,SI0,[Qs1,Qs2,Qs3,Qs4,Qs5,Qs6,Qs]),phi=0..2*Pi,Qs_W_m2=0..500);N:=100:Qs_max_:=max(seq(evalf(subs(LEO_es,phi=2*Pi*i/N,eqAface,dat,SI0,Qs)),i=1..N))*W_;plot(subs(GEO_es,eqAface,dat,SI0,[Qs1,Qs2,Qs3,Qs4,Qs5,Qs6,Qs]),phi=0..2*Pi,Qs_W_m2=0..500);N:=100:Qs_max_:=max(seq(evalf(subs(GEO_es,phi=2*Pi*i/N,eqAface,dat,SI0,Qs)),i=1..N))*W_;

`assign`(Qs0, `*`(alpha, `*`(A, `*`(E))))
`assign`(eqAface, A = `*`(`^`(L, 2)))
`assign`(Qs0_, `+`(`*`(342.50, `*`(W_))))
`assign`(Qs1, `*`(alpha, `*`(A, `*`(E, `*`(piecewise(`or`(`<`(`+`(`*`(6.28, `*`(frac(`+`(`*`(.16, `*`(phi))))))), `+`(`*`(`/`(1, 2), `*`(Pi)))), `<`(`+`(`*`(`/`(3, 2), `*`(Pi))), `+`(`*`(6.28, `*`(fra...
`assign`(Qs2, `*`(alpha, `*`(A, `*`(E, `*`(piecewise(`<`(`+`(`*`(2, `*`(Pi)), `-`(phi_es)), `+`(`*`(6.28, `*`(frac(`+`(`*`(.16, `*`(phi)))))))), `+`(`-`(sin(phi))), 0))))))
`*`(alpha, `*`(A, `*`(E, `*`(piecewise(`<`(`+`(`*`(6.283185308, `*`(frac(`+`(`*`(.1591549430, `*`(phi))))))), phi_es), sin(phi), 0)))))
`+`(Qs1, Qs2, Qs3, Qs4, Qs5, Qs6)
`+`(`*`(484.1291387, `*`(W_)))
`+`(`*`(484.1291387, `*`(W_)))

i.e. one frontal face gets Q0=343 W, but the maximum takes place at phi=45º, Qs_max=sqrt(2)*Q0=484 W. Notice the difference in solar loads from LEO to GEO.

d) Infrared input from the planet, assumed at a temperature of 288 K and with =0.6.

If the planet is assumed isothermal, the infrared input is constant along the orbit. The view factor from each face to a sphere is found in View factor tabulations.

Notice that face 1 is not seen the planet, but all the four lateral faces get some radiation.

> Qpi:=alpha*'A'*Fip*epsilon*sigma*Tp^4;F1p:=0;F2p:=(1/Pi)*(arctan(1/sqrt(h^2-1))-sqrt(h^2-1)/h^2);F3p:=1/h^2;F4p:='F2p';F5p:='F2p';F6p:='F2p';Qp:='alpha*A*(F1p+F2p+F3p+F4p+F5p+F6p)*epsilon*sigma*Tp^4';F2p_LEO:=evalf(subs(LEOh,F2p));F3p_LEO:=evalf(subs(LEOh,F3p));Qp_LEO:=subs(dat,evalf(subs(LEOh,eqAface,dat,Qp)));F2p_GEO:=evalf(subs(GEOh,F2p));F3p_GEO:=evalf(subs(GEOh,F3p));Qp_GEO:=subs(dat,evalf(subs(GEOh,eqAface,dat,Qp)));

`*`(alpha, `*`(A, `*`(Fip, `*`(epsilon, `*`(sigma, `*`(`^`(Tp, 4)))))))
`/`(`*`(`+`(arctan(`/`(1, `*`(`^`(`+`(`*`(`^`(h, 2)), `-`(1)), `/`(1, 2))))), `-`(`/`(`*`(`^`(`+`(`*`(`^`(h, 2)), `-`(1)), `/`(1, 2))), `*`(`^`(h, 2)))))), `*`(Pi))
`/`(1, `*`(`^`(h, 2)))
`*`(alpha, `*`(A, `*`(`+`(F1p, F2p, F3p, F4p, F5p, F6p), `*`(epsilon, `*`(sigma, `*`(`^`(Tp, 4)))))))
`+`(`*`(126.8637154, `*`(W_)))
`+`(`*`(1.507479871, `*`(W_)))

i.e. IR inputs are 127 W at LEO and 1.5 W at GEO.

e) Albedo input along the orbit, assuming an albedo of 0.3 and a simple albedo model.

Maximum albedo input occurs at subsolar position and is Qa_max:=alpha*A*Fbp*rho*E.

For the variation of albedo input with orbit angle, a simple model may be that it is proportional to the lit area seen (i.e. the crescent) squared (to account for the fact that the smaller the crescent, the smaller the irradiance got per normal surface). However, as the fraction of lit area seen depends on orbit altitude, a simpler model may be to compute the crescent area from afar (i.e. its orthogonal proyection) and multiply it (its square) by a parabolic function that ends at the eclipse; namely:

> Fcrescent:=(1+cos(phi))/2;Fa:=Fcrescent^2*piecewise(p(phi)<phi_es,(1-(p(phi)/phi_es)^2),p(phi)>(2*Pi-phi_es),(1-((2*Pi-p(phi))/phi_es)^2),0);

`assign`(Fcrescent, `+`(`/`(1, 2), `*`(`/`(1, 2), `*`(cos(phi)))))
`*`(`^`(`+`(`/`(1, 2), `*`(`/`(1, 2), `*`(cos(phi)))), 2), `*`(piecewise(`<`(`+`(`*`(6.283185308, `*`(frac(`+`(`*`(.1591549430, `*`(phi))))))), phi_es), `+`(1, `-`(`/`(`*`(39.47841761, `*`(`^`(frac(`+...

The view factors from a face to the planet are already computed (Fip).

> Qai:='alpha*A*Fip*rho*E*Fa';Qa1:=0;Qa2:=subs(Fip=F2p,Qai):Qa3:=subs(Fip=F3p,Qai):Qa4:=subs(Fip=F4p,Qai):Qa5:=subs(Fip=F5p,Qai):Qa6:=subs(Fip=F6p,Qai):Qa:='Qa1+Qa2+Qa3+Qa4+Qa5+Qa6';plot(subs(LEOh,LEO_es,eqAface,dat,SI0,[Qa1,Qa2,Qa3,Qa4,Qa5,Qa6,Qa]),phi=0..2*Pi,Qa_W_m2=0..250);N:=100:Qa_max_:=max(seq(evalf(subs(LEOh,LEO_es,phi=2*Pi*i/N,eqAface,dat,SI0,Qa)),i=1..N))*W_;plot(subs(GEOh,GEO_es,eqAface,dat,SI0,[Qa1,Qa2,Qa3,Qa4,Qa5,Qa6,Qa]),phi=0..2*Pi,Qa_W_m2=0..3);N:=100:Qa_max_:=max(seq(evalf(subs(GEOh,GEO_es,phi=2*Pi*i/N,eqAface,dat,SI0,Qa)),i=1..N))*W_;

`*`(alpha, `*`(A, `*`(Fip, `*`(rho, `*`(E, `*`(Fa))))))
`+`(Qa1, Qa2, Qa3, Qa4, Qa5, Qa6)
`+`(`*`(222.7793783, `*`(W_)))
`+`(`*`(2.647214198, `*`(W_)))

i.e. the satellite gets a maximum of 223 W when at the subsolar point in LEO. Notice that planet effects (IR and albedo) can be neglected at very high orbits like GEO.

f) Periodic temperature evolution, assuming the body is isothermal and has a mass of 100 kg with an average thermal capacity of 1000 J/(kg·K)..

The energy balance is dE/dT=W=Q=0+Q, since there is no mechanical or electrical work. For instance, in the LEO case:

> eqEB:=dE/dt=Q[net,inp]-Q[net,out];eqEB:=m*c*diff(T(t),t)=Qs+Qa+Qp-6*A*sigma*T(t)^4:eqEB_:=evalf(subs(LEOh,LEO_es,eqAface,dat,SI0,eqEB));

`assign`(eqEB, `/`(`*`(dE), `*`(dt)) = `+`(Q[net, inp], `-`(Q[net, out])))
`+`(`*`(50000., `*`(diff(T(t), t)))) = `+`(`*`(342.50, `*`(piecewise(`or`(`<`(`+`(`*`(6.283185308, `*`(frac(`+`(`*`(.1591549430, `*`(phi))))))), 1.570796327), `<`(4.712388981, `+`(`*`(6.283185308, `*`...
`+`(`*`(50000., `*`(diff(T(t), t)))) = `+`(`*`(342.50, `*`(piecewise(`or`(`<`(`+`(`*`(6.283185308, `*`(frac(`+`(`*`(.1591549430, `*`(phi))))))), 1.570796327), `<`(4.712388981, `+`(`*`(6.283185308, `*`...
`+`(`*`(50000., `*`(diff(T(t), t)))) = `+`(`*`(342.50, `*`(piecewise(`or`(`<`(`+`(`*`(6.283185308, `*`(frac(`+`(`*`(.1591549430, `*`(phi))))))), 1.570796327), `<`(4.712388981, `+`(`*`(6.283185308, `*`...
`+`(`*`(50000., `*`(diff(T(t), t)))) = `+`(`*`(342.50, `*`(piecewise(`or`(`<`(`+`(`*`(6.283185308, `*`(frac(`+`(`*`(.1591549430, `*`(phi))))))), 1.570796327), `<`(4.712388981, `+`(`*`(6.283185308, `*`...

a first-order ordinary differential equation which must be solved with some initial conditions (e.g. T(0)=300 K) until transients decay and a periodic solution remains.

This can be done by Euler's method or better by some Runge-Kutta method.

Let us run the numerical simulation for N=5 orbits, with Dt=100 s, starting with T0=300 K. The number of discrete t-points are N*Tperiod/Dt.

> To:=subs(SI0,rhs(T_LEO));N:=9;Dt:=100;TE[1]:=300;Np:=round(N*To/Dt);for i from 1 to Np-1 do TE[i+1]:=TE[i]+evalf(subs(phi=2*Pi*Dt*i/To,alpha=1,LEOh,LEO_es,eqAface,dat,SI0,Dt*(Qs+Qa+Qp-6*A*sigma*TE[i]^4)/(m*c)));od:plEuler:=plot([seq([Dt*i,TE[i]],i=1..Np)]):sol_:=dsolve([evalf(subs(phi=2*Pi*t/To,SI0,eqEB_)),T(0)=300],type=numeric,range=0..N*To);plRunge:=odeplot(sol_,color=black):display([plEuler,plRunge]);

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error

We see the 5th orbit is already periodic and plot it aside (notice that the initial and final T-values must coincide to be periodic, but this is not the mean temperature, Tmean.

> i:='i':Tmean:=sum(TE[i],i=round(Np*(N-1)/N)..Np)/(Np+1-round(Np*(N-1)/N));Tnax:=max(seq(TE[i],i=round(Np*(N-1)/N)..Np));Tnin:=min(seq(TE[i],i=round(Np*(N-1)/N)..Np));plmean:=plot([[0,Tmean],[N*To,Tmean]]):display([plEuler,plRunge,plmean],view=[(N-1)*To..N*To,250..300]);


i.e, once in periodic state, the body temperature oscillates between 262 K and 275 K with a period average of 269 K (i.e. between -11 ºC and 2 ºC, with a mean of -4 ºC).

In the case of GEO

> To:=subs(SI0,rhs(T_GEO));N:=5;Dt:=1000;TE[1]:=300;Np:=round(N*To/Dt);for i from 1 to Np-1 do TE[i+1]:=TE[i]+evalf(subs(phi=2*Pi*Dt*i/To,GEOh,GEO_es,eqAface,dat,SI0,Dt*(Qs+Qa+Qp-6*A*sigma*TE[i]^4)/(m*c)));od:plEuler:=plot([seq([Dt*i,TE[i]],i=1..Np)]):sol_:=dsolve([evalf(subs(phi=2*Pi*t/To,GEOh,GEO_es,eqAface,dat,SI0,eqEB)),T(0)=300],type=numeric,range=0..N*To);plRunge:=odeplot(sol_,color=black):display([plEuler,plRunge]);i:='i':Tmean:=sum(TE[i],i=round(Np*(N-1)/N)..Np)/(Np+1-round(Np*(N-1)/N));Tnax:=max(seq(TE[i],i=round(Np*(N-1)/N)..Np));Tnin:=min(seq(TE[i],i=round(Np*(N-1)/N)..Np));plmean:=plot([[0,Tmean],[N*To,Tmean]]):display([plEuler,plRunge,plmean],view=[(N-1)*To..N*To,200..300]);

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if `<`(1, nargs) then error

Notice how pronunced the effect of eclipse may be on geostationary orbits (fortunately, there are only GEO eclipses near the equinocces, contrary to LEO, where eclipse duration varies little with season).