1
$\begingroup$

For some reason I need to build a Elliptical Orbits satellite dynamic model by myself, in here I am using MATALB.
The model its fine before (I think the J2 effect part works well) I add drag effect on it, the Altitude became ascent, it should descent by time right?
Here are some equation I am using:
Drag
And the simulator result:
enter image description here enter image description here
Drag part code:

% Calculate the drag acceleration
    v_rel_mag = norm(v_rel);
    rho = atmospheric_density(norm(r_icrf) - RE); % call to your atmospheric density function
    a_drag = -0.5 * Cd * A / m *  rho * v_rel_mag^2 * (v_rel / v_rel_mag)
    
    % Total acceleration
    a_total = a_gravity + a_J2 + a_drag;

    % Output derivative of state vector
    Ydot = [vvector; a_total];
end

function rho = atmospheric_density(altitude)
    % Simple exponential atmosphere model (or use a more complex model)
    rho0 = 3.614*10^-14; % kg/m^3 at sea level
    H = 8500; % scale height in meters
    rho = rho0 * exp(-((altitude-700) / 88.667));
end

[Full code][4] https://pastecode.io/s/pcw2mvjb
Can anyone help me? thank you!
$\endgroup$
3
  • $\begingroup$ What is the integration scheme? Matlab's builtin ode45 or similar? $\endgroup$
    – AJN
    Commented Dec 15, 2023 at 11:12
  • $\begingroup$ @AJN I am using ODE23s $\endgroup$ Commented Dec 15, 2023 at 11:32
  • 2
    $\begingroup$ I think there's a nearly-identical existing question on the site which makes it clearer that the problem is that the satellite is seeming to ascend because of atmospheric drag. I believe AJN's indication that integration precision issues can cause that was the resolution to that problem; I don't have time at the moment to figure out which question it was, but you'd be well-served looking for it on this site yourself. $\endgroup$
    – Erin Anne
    Commented Dec 15, 2023 at 21:31

1 Answer 1

3
$\begingroup$

When coding/scripting a new project, always independently test each of your functions. Don't just run the whole thing at once.

Also, the single scale height is useless at orbital altitudes. Above about 100 km the fall-off is MUCH slower.

atmospheric knee at 100 km

"rho0 = 3.614*10^-14; % kg/m^3 at sea level" seems way off. It's about 1.3 mg/cm^3 and so about 1 .3 kg/m^3. And two lines below that you don't even use the scale height H (which looks good) but instead use 88.667 and have no idea what that number is.

Your Matlab:

function rho = atmospheric_density(altitude)
    % Simple exponential atmosphere model (or use a more complex model)
    rho0 = 3.614*10^-14; % kg/m^3 at sea level
    H = 8500; % scale height in meters
    rho = rho0 * exp(-((altitude-700) / 88.667));
end

My Python:

import numpy as np
import matplotlib.pyplot as plt

def atmospheric_density(altitude):
    # Simple exponential atmosphere model (or use a more complex model)
    rho0 = 3.614E-14 # kg/m^3 at sea level
    H = 8500 # scale height in meters
    rho = rho0 * np.exp(-((altitude - 700) / 88.667))
    return rho

def atmospheric_density_fixed(altitude):
    rho0 = 1.3 # kg/m^3 at sea level
    H = 8500 # scale height in meters
    rho = rho0 * np.exp(-((altitude - 700) / H))  # DIVIDE BY H
    return rho

alti = np.linspace(0, 4E+05, 1000) # sea level to 400 km

density = atmospheric_density(alti)
density_fixed  = atmospheric_density_fixed(alti)

densities = density, density_fixed
labels = 'original', 'fixed'

fig, ax = plt.subplots(1, 1)

for den, label in zip(densities, labels):
    ax.plot(alti/1000, den, label=label)
ax.set_yscale('log')
ax.set_xlabel('altitude (km)')
ax.set_ylabel('density (kg/m^3)')
ax.set_ylim(1E-30, 10)
ax.legend()
plt.show()

Result:

enter image description here

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.