Flat-Spin-Übergangsmanöver: Warum ist der Drehimpuls nicht konstant?

Kurze Antwort: Es handelt sich um ein ODE-Solver-Problem.

Beim Einstellen der relativen und absoluten Fehlertoleranzen werden bessere Ergebnisse erzielt. Beispielsweise optionskann die Variable im MATLAB-Code wie folgt definiert werden:

options = odeset('RelTol',1e-10,'AbsTol',1e-10);

Um den anfänglichen Drehimpuls statt (wie in dieser Antwort darauf hingewiesen wurde) zu lassen 3000, kann3000.0045 die Gleichung zur Berechnung des Drehimpulses auch wie folgt für die anfängliche Spinrate der Nebenachse gelöst werden:

omega_z = ( ( 3000^2 - ( Ix * omega_x )^2 ) / ( Iz^2 ) )^( 1/2 );

Um schließlich die Achsen des Drehimpulsdiagramms anzupassen, fügen Sie hinzu:

axis([0 1000 2999 3001])

Daraus ergibt sich folgende Abbildung:

Ich habe Ihr Matlab in Python umgeschrieben (was wirklich einfach ist und das ist kein Zufall!) und den Standard-SciPy-ODE-Integrator verwendet. Beim Überprüfen info['mused']sehe ich, dass immer der Adams-Algorithmus (nicht steif) verwendet wurde, sodass die Steifigkeit nicht das Problem ist.

Warum nicht einfach das „alte und kaputte“ Matlab fallen lassen und das „neue heiße“ Python übernehmen? ( MIB-Referenz )

Ich würde vermuten, dass Sie anfangen sollten, Ihre Matlab-Dokumentation zu lesen ode45und Ihre Fehlertoleranz niedriger einzustellen. Ich vermute, das ode45ist, was ich RK45 nennen würde und eine variable Schrittgröße 4/5 der Ordnung Runge-Kutta ist. Es sollte in Ordnung sein, aber wenn Sie keine Toleranz angeben, wird es eine für Sie auswählen.

Es ist möglich, dass Matlab zwischen der Veröffentlichung Ihrer Referenz und dem Ausführen Ihres Beispiels einige Standardeinstellungen geändert hat oder der Autor einige verwendet und vergessen hat, Ihnen dies mitzuteilen, oder Sie haben sie nicht bemerkt, die früher in dem, was Sie gelesen haben, erwähnt wurden .

Lange Rede kurzer Sinn, bei mir funktioniert es einwandfrei.

Der Anfangsdrehimpuls ist nicht 3000 unter Verwendung der Initialisierungswerte, sondern genau 3000.0045dort, wo die Simulation beginnt. Es driftet nach 3000.0011000 Sekunden ab. Das ist eine Änderung von 0,003. Die Abweichung von geraden 3000 (oder 3e3) wird im Diagramm angezeigt.

Ich kann das ein wenig verbessern, auf eine Gesamtdrift von nur etwa 0,001, indem ich rtol=1E-10oder rtol=1E-11in den Aufruf von einfüge ODEint.

def omega_dynamics(x_0, t, I, mu, T, J):

    # % Initial state vector
    omega_x, omega_y, omega_z, sigma_x, sigma_y, sigma_z = x_0[:6]

    # % Constants
    mu_x, mu_y, mu_z = mu[:3]
    Ix, Iy, Iz = I[:3]
    Tx, Ty, Tz = T[:3]

    # % Differential equations
    omega_dot_x = ((Iy - Iz)*omega_y*omega_z + mu_x*sigma_x + Tx)/(Ix - J)
    omega_dot_y = ((Iz - Ix)*omega_z*omega_x + mu_y*sigma_y + Ty)/(Iy - J)
    omega_dot_z = ((Ix - Iy)*omega_x*omega_y + mu_z*sigma_z + Tz)/(Iz - J)
    sigma_dot_x = -omega_dot_x - mu_x*sigma_x/J - omega_y*sigma_z + omega_z*sigma_y
    sigma_dot_y = -omega_dot_y - mu_y*sigma_y/J - omega_z*sigma_x + omega_x*sigma_z
    sigma_dot_z = -omega_dot_z - mu_z*sigma_z/J - omega_x*sigma_y + omega_y*sigma_x

    # % Return vector
    dx = np.array([omega_dot_x, omega_dot_y, omega_dot_z, sigma_dot_x, sigma_dot_y, sigma_dot_z])

    return dx

# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# %%%%% Example Rahn & Barba (1991) %%%%%
# %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

J        = 18.0;                                             # Principal moments of inertia of spherical propellant slug
I        = np.array([2000,   1500,   1000], dtype=float)   # Principal moments of inertia of spacecraft including slug
omega_0  = np.array([0.1224, 0,      2.99], dtype=float)   # Spacecraft body rates
sigma_0  = np.array([0,      0,      0   ], dtype=float)   # Relative rates between spacecraft body and propellant slug
mu       = np.array([30,     30,     30  ], dtype=float)   # Viscous damping coefficient
T        = np.array([0,      0,      0   ], dtype=float)   # Torques about principal axes

h_init   = np.sqrt(((I * omega_0 + J * sigma_0)**2).sum())  # % Calculate angular momentum

print "h_init: ", h_init

times        = np.linspace(0, 1000, 10001)
X0           = np.hstack((omega_0, sigma_0))

x1, info = ODEint(omega_dynamics, X0, times, args=(I, mu, T, J), full_output=True)
print x1.shape

h          = np.sqrt(((I * x1[:, :3] + J * x1[:,3:])**2).sum(axis=1))  # % Calculate angular momentum

if True:
    plt.figure()

    plt.subplot(2, 2, 1)
    plt.plot(times, x1[:, 2])
    plt.ylim(-2, 3)
    plt.title('Minor Axis Spin Rate')
    plt.xlabel('Time - seconds')
    plt.ylabel('rad/s')

    plt.subplot(2, 2, 2)
    plt.plot(times, x1[:, 1])
    plt.title('Intermediate Axis Spin Rate')
    plt.xlabel('Time - seconds')
    plt.ylabel('rad/s')

    plt.subplot(2, 2, 3)
    plt.plot(times, x1[:, 0])
    plt.title('Major Axis Spin Rate')
    plt.xlabel('Time - seconds')
    plt.ylabel('rad/s')

    plt.subplot(2, 2, 4)
    plt.plot(times, h)
    plt.title('Angular Momentum')
    plt.xlabel('Time - seconds')
    plt.ylabel('Nms')

    plt.show()