Gelöste Aufgaben/Kerb

Aufgabenstellung

Ein Satellit soll eine stabile Umlaufbahn um die Erde beschrieben.

Berechnen Sie 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.

tmp

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            */
/*********************************************************/

tmp

Declarations

Wir brauchen die Systemparameter

Symbol	Bedeutung
R = 6378,137 km	Erdradius am Äquator
g₀ = 9.81 m/s2	Erdbeschleunigung an der Erdoberfläche
T = 1 d	Periodendauer einer Erdumdrehung (d .. day)
O_r = 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).

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);

tmp

Equilibrium Conditions

..

/******************************************/
/* 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));

tmp

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}}$ .===Stationary Solution=== ..

XX

tmp

Solving

..

XX

tmp

Post-Processing

..

XX

Links

...

Literature

...

Trajectory animation