Atmospheric drag effect

Question

While propagating the satellite motion faced the strange effect. Without the atmosphere results (x,y,z coordinates in meters and vx,vy,vz velocities in meters per second) are [3.30971e5, -6.55755e6, 2.57717e6, -1261.46, -2762.24, -6871.88] but including the atmosphere drag effect i obtain [-2.90003e5, -7.04404e6, -8.97772e5, -1268.31, 996.889, -7305.6]. Notice that difference between satelite positions is 300+ km. Can the atmospheric drag effect be so big or it is the implementation and modelling issues? The formula for atmospheric drag was taken from GMAT documentation:$$a=\frac{1}{2}\rho v^2_{rel}\frac{C_dh}{m_s}\hat{v_{rel}}$$ and the barometric formula(from Wikipedia) to calculate the density: $$\rho=\rho_d\exp\left[ \frac{-g_0M(h-h_0)}{ RT_b}\right ].$$ The $h_0,T_b$ and $\rho_d$ constants were taken from http://www.braeunig.us/space/atmos.htm. Here is the code:

using DifferentialEquations
jd  = Dates.datetime2julian(DateTime(2100,01,01,0,0,0)) * 86400
jd2 = Dates.datetime2julian(DateTime(2100,01,11,0,0,0)) * 86400
y = [-976.3107644649057E+03, -4835.627052558522E+03, -5031.728586125443E+03, -0.7944031487816871E+03, 5.474532271429767E+03, -5.094496750907486E+03]
GMe = BigFloat(398600.4415E+9)
req = BigFloat(6378136.3) #m
mass = 850
A = 15.0 #m^2
rho_b  = 6.65E-14 #kg/m^3,
g   = 9.80665 #m/s^2,
M = 0.0289644 #kg/mol,
h_0 = 650000 #m,
R = 8.3144598 #N·m/(mol·K),
T_b = 1011.5365 #K
omega = [0.0, 0.0, 7.292115078468551e-5]
v_wind = [0.0,0.0,0.0]
C_d = 0.47#from GMAT

#× vector product
function atmospheric_drag(y)
    v = [y[4],y[5],y[6]]
    r = [y[1],y[2],y[3]]
    h = sqrt(y[1]^2+y[2]^2+y[3]^2) - req
    rho = rho_b*exp((-g*M*(h-h_0))/(R*T_b))
    #println("rho = $rho")
    v_rel = v - omega×r + v_wind
    v_rel2 = v_rel[1]^2+v_rel[2]^2+v_rel[3]^2
    a = -0.5*rho*(C_d*A/mass)*v_rel2*v_rel
end

println(atmospheric_drag(y))

function f2(dy,y,p,t)

    date=t/86400
    re3=(y[1]^2+y[2]^2+y[3]^2)^(3/2)
    atm_dr = atmospheric_drag(y)
    dy[1] = y[4]
    dy[2] = y[5]
    dy[3] = y[6]
    dy[4] = - GMe*y[1]/re3 - atm_dr[1]
    dy[5] = - GMe*y[2]/re3 - atm_dr[2]
    dy[6] = - GMe*y[3]/re3 - atm_dr[3]
end

prob = ODEProblem(f2,y,(jd,jd2))
solution = solve(prob,Vern9(),abstol=1e-13,reltol=1e-13)
sol = solution[end]
println(sol)
#order 0 degree 0 no Sun, no Moon
GMATresult = [330.97011851316,-6557.5472439174,2577.1624246640,-1.2614569607713,-2.7622388770690,-6.8718847034529]

a = 0
println("Error:")
for i = (1:6)
  b = sol[i] - GMATresult[i]*1000
  println("$b m")
	if i < 4
	 a = a + b^2
	end
end
println("Distance error = $(sqrt(a)) m")

I've noticed that even little values like 1e-6 being added to f2 affect significantly (70+ km) on the result. The values of the acceleration provided by atmospheric_drag function are of 1e-5 order. Maybe that is the reason?

Answer 1

Back of the spherical cow's envelope:

It's moving at $v$ = 7500 m/s. After ten days it's off by 3E+05 meters.

$$x = \frac{1}{2} a t^2$$

gives an acceleration of 8E-08 m/s^2. With a mass $m$ = 850 kg that's a force of 7E-04 Newtons.

Using a simple drag equation for $A$ = 10 m^2 and $C_D$ = 1, plugging that force into

$$F_D = \frac{1}{2} \rho v^2 C_D A$$

and solving for $\rho$ gives 2.5E-12 kg/m^3.

Looking that up in this answer or the reference there gives an altitude of only 450 km, which means it's ballparkish reasonable.

However, I don't like your simple scale height-ish atmosphere expression for density versus height. The exponential factor $ R T / (g M) $ is about 30 km with a temperature of 1000 Kelvin and I'm not sure why that works even as well as it does because you are using an $h_0$ of sea level. With such a low slope, I think you are meant to use a $\rho_d$ and $h_0$ that's already very high in the atmosphere.

Oh! If you are just using that from Wikipedia, then it's just plain wrong. You'll need to find a model for drag that's more appropriate for space!

So instead you could of course just interpolate density from your own link http://www.braeunig.us/space/atmos.htm or even just use their scale height values at say 600 km:

altitude:      600       km
scale height:  74.8      km
rho_0          9.89E-14  kg/m^2

and use $$ \rho = 9.89 \times 10^{-14} \exp\left(\frac{r - (6378.137 + 600) km}{74.8 km}\right). $$

That should give you a roughly factor of 10 smaller difference in drag vs no-drag, maybe a 10 to 50 km instead of 300 km.

In a 0th order, spherical cowish kind of way, you could even use an interpolator on the data I've shown there plus the drag equation plus $a = F/m$ to get a better acceleration than you're using now at least.

Good luck!

below: ugly sketch from here to illustrate that you can't use a single, simple scale height formula from sea level to LEO.