> restart:#"m13_p37"

Find the electrical power produced by an spherical satellite of 0.5 m in diameter, fully covered by solar cells of an efficiency =15% and a packaging factor Fpg=0.8, in a low Earth orbit without eclipses, and set the thermal balance, assuming an absorptance and emissivity of ==0.75 for the solar cells, a thermal capacity of C=30 kJ/K for the satellite, and that the electrical dissipation is only important during 15 minutes of the orbit, and can be considered constant in that period.


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

> dat:=[D=0.5*m_,eta=0.15,Fpg=0.8,H=300e3*m_,R=6370e3*m_,E=1370*W_/m_^2,To=90*60*s_,alpha=0.75,epsilon=0.75,C=30e3*J_/(kg_),tdis=15*60*s_];

[D = `+`(`*`(.5, `*`(m_))), eta = .15, Fpg = .8, H = `+`(`*`(0.300e6, `*`(m_))), R = `+`(`*`(0.6370e7, `*`(m_))), E = `+`(`/`(`*`(1370, `*`(W_)), `*`(`^`(m_, 2)))), To = `+`(`*`(5400, `*`(s_))), alpha...
[D = `+`(`*`(.5, `*`(m_))), eta = .15, Fpg = .8, H = `+`(`*`(0.300e6, `*`(m_))), R = `+`(`*`(0.6370e7, `*`(m_))), E = `+`(`/`(`*`(1370, `*`(W_)), `*`(`^`(m_, 2)))), To = `+`(`*`(5400, `*`(s_))), alpha...

> dat:=op(dat),Const,SI2,SI1:


Fig. 1. Evolution of input loads (Qele_dis+Qsolar) and satellite temperature (T) after an initial value of 300 K. Detail of the periodic temperature evolution and its orbit mean.

a) Find the electrical power produced...

The power produced by the solar cells is:

> Wele_gen:=eta*Fpg*E*Afrontal;eqA:=A=Pi*D^2;Afrontal:=A/4;Weke_gen_:=subs(dat,evalf(subs(eqA,dat,Wele_gen)));Wele_dis:=Wele_gen*To/tdis;Wele_dis_:=subs(dat,evalf(subs(eqA,dat,Wele_dis)));

`*`(eta, `*`(Fpg, `*`(E, `*`(Afrontal))))
A = `*`(Pi, `*`(`^`(D, 2)))
`+`(`*`(`/`(1, 4), `*`(A)))
`+`(`*`(32.279864515635125526, `*`(W_)))
`+`(`/`(`*`(`/`(1, 4), `*`(eta, `*`(Fpg, `*`(E, `*`(A, `*`(To)))))), `*`(tdis)))
`+`(`*`(193.67918709381075315, `*`(W_)))

i.e. 32 W are being produced all the time during the typical 90 minutes of a LEO period. If this electrical energy is to be consumed in 15 min, the rate of consumption must be 32•90/15=194 W.

The thermal balance takes the form:

> eqEBth:=C*diff(T(t),t)=Qele_dis+Qsolar-Qout;eqEBth:=C*diff(T(t),t)=Qele_dis+alpha[th]*E*Afrontal-epsilon*A*sigma*T(t)^4;Qele_dis:=piecewise(To*frac(t/To)<tdis,'Wele_dis',0);

`assign`(eqEBth, `*`(C, `*`(diff(T(t), t))) = `+`(Qele_dis, Qsolar, `-`(Qout)))
`assign`(eqEBth, `*`(C, `*`(diff(T(t), t))) = `+`(Qele_dis, `*`(`/`(1, 4), `*`(alpha[th], `*`(E, `*`(A)))), `-`(`*`(epsilon, `*`(A, `*`(sigma, `*`(`^`(T(t), 4))))))))
piecewise(`<`(`*`(To, `*`(frac(`/`(`*`(t), `*`(To))))), tdis), Wele_dis, 0)

Notice that a periodic step function must be built (equal to 0 except for a 15 min period during the orbit) if more than one orbit is to be simulated (what is always needed if the initial value problem does not happen to coincide with the boundary value problem in the periodic behaviour).

The solar thermal absorption, and the satellite infrared emission are:

> Qsolar:=alpha[th]*E*Afrontal;Qsolar:=(alpha-eta*Fpg)*E*Afrontal;Qsolar_:=subs(dat,evalf(subs(eqA,dat,%)));Qout_:=evalf(subs(eqA,dat,epsilon*A*sigma*T(t)^4));eqEBth_:=evalf(subs(eqA,dat,SI0,lhs(eqEBth)=Qele_dis+Qsolar_-Qout_));

`assign`(Qsolar, `+`(`*`(`/`(1, 4), `*`(alpha[th], `*`(E, `*`(A))))))
`assign`(Qsolar, `*`(`*`(`/`(1, 4), `+`(alpha, `-`(`*`(eta, `*`(Fpg))))), `*`(E, `*`(A))))
`assign`(Qsolar_, `+`(`*`(169.47, `*`(W_))))
`assign`(Qout_, `+`(`/`(`*`(3.34, `*`(`^`(10, -8), `*`(kg_, `*`(`^`(m_, 2), `*`(`^`(T(t), 4)))))), `*`(`^`(s_, 3), `*`(`^`(K_, 4))))))
`+`(`*`(0.30e5, `*`(diff(T(t), t)))) = `+`(piecewise(`<`(`+`(`*`(5400., `*`(frac(`+`(`*`(0.18518518518518518519e-3, `*`(t))))))), 900.), 193.67918709381075315, 0.), 169.46928870708440901, `-`(`*`(0.33...

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 (independent on T0) remains.

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

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

> N:=9;Dt:=100;TE[1]:=300;To:=subs(dat,SI0,To);Np:=round(N*To/Dt);for i from 1 to Np-1 do TE[i+1]:=TE[i]+evalf(subs(t=i*Dt,eqA,dat,SI0,Dt*(Qele_dis+Qsolar_-epsilon*A*sigma*TE[i]^4)/C));od:plEuler:=plot([seq([Dt*i,TE[i]],i=1..Np)]):sol_:=dsolve([eqEBth_,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 9th 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 277 K and 282 K with a period average of 279 K (i.e. between 4 ºC and 9 ºC, with a mean of 6 ºC).

Notice that spacecraft temperature increases during electrical dissipation and decreases otherwise (when infrared emission surpasses solar input).