Gelöste Aufgaben/Kerb

Aufgabenstellung

Ein Satellit soll eine stabile Umlaufbahn um die Erde beschreiben.

Gesucht sind mögliche Lösungen.

Lösung mit Maxima

Bei unserer Suche geht es um kartesische- und polare-Koordinatensysteme, Gleichgewichtsbedingungen formuliert mit dem Prinzip der virtuellen Verrückungen und er Stabilität von Lösungen.

Header

Ein Satellit im Orbit ist einer Erdbeschleunigung g ausgesetzt, die nichtlinear von seinem Abstand zum Erdmittelpunkt abhängt. Dieses nichtlineare Anfangswertproblem lösen wir hier als Anfangswertproblem numerisch.

/*********************************************************/
/* MAXIMA script                                         */
/* version: wxMaxima 15.08.2                             */
/* author: Andreas Baumgart                              */
/* last updated: 2018-03-21                              */
/* ref: Kerb (TM-C, Labor 6)                             */
/* description: finds persiod solution for               */
/*              the tracectory of a satellite            */
/*********************************************************/

Declarations

Wir brauchen die Systemparameter 

Symbol	Bedeutung
R = 6378,137 km	Erdradius am Äquator
g0 = 9.81 m/s2	Erdbeschleunigung an der Erdoberfläche
T = 1 d	Periodendauer einer Erdumdrehung (d .. day)
Or = 35786 km	Höhe eines Satelliten über der Erdoberfläche in einem geostationären Orbit.

/*********************************************************/
/* declarations */
/* declare variational variables - see 6.3 Identifiers   */
declare("δW", alphabetic);
declare("δr", alphabetic);
declare("δu", alphabetic);
declare("δφ", alphabetic);

assume(R>0, g>0, rho>0, mu>0, T>0);

/* parameter */
params: [R = 6378137*m,    /* earth radius at equator */
         g = 9.81*m/s^2,   /* gravitational constant  */
         T = 24*60*60*s];  /* 24 hours                */
verify: [O=R+35786000*m];  /* geostat. orbit          */
dimles: [g=mu*R/T^2];      /* dimensionless parameter */

Kinematics

Wir führen ein kartesisches x, y, z - Koordinaten-System im Erdmittelpunkt (earth-centered inertial) ein, das nicht mit der Erde mitrotiert - wir setzen dieses als Intertialsystem an (obwohl wir wissen, dass die Erde um die Sonne rotiert und die Sonne selbst auch kein Intertialsystem in unserem Sinne ist).

Satellite in Orbit

Wir brauchen die Einheitsvektoren

${\vec {e}}_{x},{\vec {e}}_{y}$

entlang der x,y-Achsen, um den Ortsvektor zum Satelliten

${\vec {r}}=r(t)\cdot {\vec {e}}_{r}(t)$

zu erfassen.

Dafür verwenden wir

${\vec {e}}_{r}(t)=\cos(\varphi (t))\cdot {\vec {e}}_{x}+\sin(\varphi (t))\cdot {\vec {e}}_{y}$.

Die Bewegung des Satelliten können wir also auch in kartesischen Koordinarten

${\begin{array}{l}u_{x}(t)=r(t)\cdot \cos(\varphi (t))\\u_{y}(t)=r(t)\cdot \sin(\varphi (t))\end{array}}$,

anschreiben, was hier allerdings wenig anschaulich ist! Für die Variation von u_x, u_y erhalten wir

${\begin{array}{l}\delta u_{x}=\cos(\varphi (t))\cdot \delta r-r(t)\cdot \sin(\varphi (t))\cdot \delta \varphi \\\delta u_{y}=\sin(\varphi (t))\cdot \delta r+r(t)\cdot \cos(\varphi (t))\cdot \delta \varphi \end{array}}$

/******************************************/
/*  kinematics                            */
/******************************************/
kinem : [u[x] = r(t) * cos(phi(t)), u[y] = r(t) * sin(phi(t))];
varia : [r(t) = r(t) + epsilon*δr, phi(t) = phi(t) + epsilon*δφ];
varia :  [δu[x],δu[y]] = subst([epsilon=0],diff(subst(varia,subst(kinem,[u[x],u[y]])),epsilon));
kinem : append(makelist(lhs(varia)[i] = rhs(varia)[i],i,1,2), kinem);

Equilibrium Conditions

Die Bewegungsgleichungen schreiben wir mit dem Prinzip der virtuellen Verrückungen an, hier mit den Anteilen der D'Alembert'sche Trägheitskräfte und der Gewichtskraft

$\delta W=\underbrace {-m\cdot {\ddot {u}}_{x}(t)\cdot \delta u_{x}-m\cdot {\ddot {u}}_{y}(t)\cdot \delta u_{y}} _{\text{d'Alembert'sche Trägheitskräfte}}\;\;\;\underbrace {-m\,g(r)\cdot \delta r} _{\text{Gewichtskraft}}$,

dabei ist m die Masse des Satelliten - die sich aber sofort herauskürzt.

In die Gleichgewichtsbedingung

$\delta W{\stackrel {!}{=}}0$

setzen wir alle Größen ein und finden

$\delta W=\displaystyle -{\frac {\left(g_{0}\,m\,{{R}^{2}}+m\cdot {r^{2}}\,{\ddot {r}}(t)-m\cdot {r^{3}}\,{\dot {\varphi }}^{2}\right)\,\delta r+\left(2\,m\cdot {r^{3}}\,{\dot {\varphi }}\,{\dot {r}}+m\cdot r^{4}\,{\ddot {\varphi }}\right)\,\delta \varphi }{r^{2}}}$.

Das Trennen nach den virtuellen Verrückungen liefert dann die beiden Bewegungsgleichungen

${\begin{array}{cl}\displaystyle -r\cdot {\dot {\varphi }}^{2}+{\ddot {r}}+g_{0}\cdot \left({\frac {R}{r}}\right)^{2}&=0,\\-2\,r\cdot {\dot {\varphi }}\cdot {\dot {r}}-r^{2}\cdot {\ddot {\varphi }}&=0\end{array}}$.

Für die dimensionslose Schreibweise setzen wir

${\begin{array}{ll}t&=\tau \cdot T\\r(t)&=\rho (t)\cdot R\\\end{array}}$

an und verwenden die Abkürzungen

${\begin{array}{ll}\displaystyle {\frac {d}{d\tau }}\varphi &=\Phi \\\displaystyle {\frac {d}{d\rho }}\varphi &=P\end{array}}$,

wobei das P für das große griechische ρ steht. Damit - und mit dem dimensionslosen Parameter μ statt

$\displaystyle g={\frac {\mu \,R}{T^{2}}}$

lesen sich die beiden Bewegungsgleichungen so:

${\begin{array}{ll}\displaystyle {\frac {d}{d\tau }}P&\displaystyle ={\frac {\Phi ^{2}\,\rho ^{3}-\mu }{\rho ^{2}}}\\\displaystyle {\frac {d}{d\rho }}\Phi &=\displaystyle -2{\frac {\Phi \,P}{\rho }}\end{array}}$.

/******************************************/
/* define ODE in dim'less coordinates     */
/******************************************/
/* principle of virtual work              */
PvV:  δW = - m*'diff(u[x],t,2)*δu[x] - m*'diff(u[y],t,2)*δu[y]  - m*g*(R/r(t))^2*δr;
PvV: subst(kinem, PvV);
PvV: trigsimp(expand(ev(PvV,nouns)));
/* equations of motion */
eom: makelist(ev(coeff(expand(subst(PvV,δW)), [δr, δφ][i]),nouns),i,1,2);
/* transformation to dimensionless coordinates */
trafo: [diff(  r(t),t,2) = R*'diff(Rho,tau)/T^2, diff(  r(t),t) = R*Rho/T,
        diff(phi(t),t,2) =   'diff(Phi,tau)/T^2, diff(phi(t),t) =   Phi/T,
               r(t)      = rho*R,                        phi(t) = phi    ]; 
eom : ratsimp(solve(subst(trafo,eom), ['diff(Rho,tau),'diff(Phi,tau)]))[1];
eom : ratsimp(subst(dimles,eom));

Stationary Solution

Für den stationären Orbit ist

$\displaystyle {\frac {d}{d\tau }}P^{*}=0{\text{ und }}P^{*}=0$,

also

$\displaystyle \Phi ^{*}={\sqrt {\frac {\mu }{\rho ^{3}}}}$.

Hier wählen wir als stationäre Lösung den geostationären Orbit - der Satellit dreht sich dann in T = einem Tag einmal um den Erdmittelpunkt und erscheint so von einem Betrachter auf der Erde als "fest"-stehend. Die dazugehörige Drehgeschwindigkeit in dimensionsloser Koordinate ist

$\displaystyle \Phi ^{*}{\stackrel {!}{=}}{\frac {2\pi }{1}}$

und damit

$\displaystyle \rho ^{*}={\sqrt[{3}]{\frac {\mu }{4\pi ^{2}}}}$

Dieser Wert stimmt gut mit bekannten Werten für die Höhe des geostationären Orbits (O = 35.786 km s.o.) überein.

/******************************************/
/* stationary solution (dphi/dt = 2*%pi/T)*/
/* (geostationary orbit) for ....         */
/******************************************/
geostat : solve(rhs(eom[1])=0,Phi)[2];
geostat : [geostat,solve(subst(geostat, Phi = 2*%pi),rho)[3]];
geostat : append(geostat, subst(params,solve(dimles,mu)));
test: float(subst(params,subst(geostat,rho)=subst(verify,O)/R));

eom : append(['diff(rho,tau)=Rho,'diff(phi,tau)=Phi],eom);

Solving

Die dimensionslosen Zustandsgrößen des Satelliten sind

${\underline {q}}=\left({\begin{array}{c}\rho \\\varphi \\P\\\Phi \end{array}}\right)$.

Die Bewegungsgleichungen können wir nun als Anfangswertproblem numerisch lösen. Wir wählen zwei Tage als Integrationsdauer, als Löser ein Runge-Kutta-Verfahren sowie drei verschiedene Anfangsbedingungen. Davon ist die erste (blau) die geostationäre Lösung.

/******************************************/
/* numerical solution of IVP              */
/******************************************/

/* state variables */
Q     : [rho,phi,Rho,Phi];
states: [r,p,v,o];

/* smapling characteristics */
times : subst([t0 = 0, tmax = 2, dt = 0.001],
                    [t, t0, tmax, dt]);
/* initial values */
init : [[subst(geostat,rho),0,0,2.0*%pi],
        [6,0,0,2.0*%pi],
        [7,0,0,2.0*%pi]];

/* transformation to nomenclature with is "digestible" for maxima */
newnom: makelist(Q[j]=states[j],j,1,4);

/* right-hand-side for the integration routine */
dydt : subst(newnom,subst(params,makelist(rhs(subst(solve(dimles,mu),eom[i])),i,1,4)));
/* test: does that work? */
test: subst(makelist(states[j]=init[1][j],j,1,4),dydt);

/* integrate for all sets of initial values */
for i: 1 thru length(init) do
   num[i] : rk(dydt,states,init[i],times)$

Post-Processing

Wir sehen, dass der Satellit bei der blauen - geostationären - Lösung in der Tat eine Kreisbahn um den Erdmittelpunkt beschriebt. Werden Anfangsbedingungen näher am Erdmittelpunkt (rot) oder weiter entfernt davon (grün) gewählt, so ergeben sich weiterhin periodische Lösungen, allerdings keine kreisförmigen Lösungen mehr um den Erdmittelpunkt.

Orbit	Animation
Orbits (nicht maßstäblich).	Video laden Lokale Datei Lokale Datei sammelt möglicherweise persönliche Daten. FortfahrenAusblenden Trajectory animation

In den Phasen-Diagrammen für die Lösungen kann man gut erkennen, wie sich dabei die Geschwindigkeiten über euinen Umlauf ändern.

Phasen-Diagramme

Analog können wir die Winkelgeschwindigkeit Φ über der Radialgeschwindigkeit P auftragen:

Phasengeschwindigkeit

Und wir können uns sowohl die Koordinate ρ also auch die Winkelgeschwindigkeit Φ im Zeitbereich ansehen:

Weg-Zeit-Diagramm

/******************************************/
/* post-processing                        */
/******************************************/
for i: 1 thru length(init) do
	(T[i] : makelist(num[i][j][1],j,1,length(num[i])),
	 r[i] : makelist(num[i][j][2],j,1,length(num[i])),
	 p[i] : makelist(num[i][j][3],j,1,length(num[i])),
	 v[i] : makelist(num[i][j][4],j,1,length(num[i])),
         o[i] : makelist(num[i][j][5],j,1,length(num[i])))$
/* set print-precision */
fpprintprec: 4$

/* plot results */
/* - orbits */
orbits:  makelist([discrete, r[i]*cos(p[i]), r[i]*sin(p[i])],i,1,length(init))$
legends: append([legend], makelist(sconcat("start @: ",string(float(init[i]))),i,1,length(init)))$
plot2d(orbits,
       [xlabel,"u/R →"], [ylabel,"v/R →"], [legend, "geostat.", "var. 1", "var 2"], [title, "Orbits"], same_xy)$

/* - phase diagrams */
phases :  makelist([discrete, r[i], v[i]],i,1,length(init))$
plot2d(phases,
       [xlabel,"r/1 →"], [ylabel,"P/1 →"], legends,
       [title, "Phase Digram for radial coord."])$

/* - speed diagrams */
v_o :  makelist([discrete, v[i], o[i]],i,1,length(init))$
plot2d(v_o, [xlabel,"v/1 →"], [ylabel,"o/1 →"], legends)$


/* - time - radius / angular speed */
ro_t : append(makelist([discrete, T[i], r[i]],i,1,length(init)),
              makelist([discrete, T[i], o[i]],i,1,length(init)))$
plot2d(ro_t, [style, [lines,1,1], [lines,1,2], [lines,1,3], [lines,2,1], [lines,2,2], [lines,2,3]],
       [xlabel,"tau/1 →"], [ylabel,"r/1 →"], [legend, "r geostat.", "r var. 1", "r var 2","o","o","o"])$

---

Links

• Kerbala Space Program

Literature

• ...