What formulas do I use to calculate the gravity and drag forces on an object ascending from Earth's surface?

Question

I'm trying to model a slingshot launch to space, literally, and I'm using 2 basic formulas. Simplifying my model, a object with a (really big) initial velocity is launched from Earth's surface on the equator, and it is affected by air resisance and gravity. I am using the following formulas:

Air resistance:

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

A = Reference area C = Drag coefficient F = Drag force, N. V = Velocity, m/s. ρ = Density of fluid (liquid or gas), kg/m3.

And : universal gravitational force:

$$F = G \frac{m_1 m_2}{r^2}$$

but I'm getting some wierd results.

I want to know if i'm using the right formulas for a object that is really low mass in comparison to Earth.

More info:

I'm using a Python program to plot a graphic that would be the X and Y position. For the method I'm using (Odeint) I need the differential equation (don't know if the name is correct in english) for it, as F = ma:

$$F/m_1 = G \frac{m_2}{r^2}$$

if $m_1$ = my object mass

so my Speed Variation/Time Variation = $G*M2 / r^2$ (being M2 = earth mass)

(I think it makes sense because a object's gravity acceleration does not depend on that object's mass.)

I'm almost sure I'm doing a basic and terrible mistake here, but I'm in my first year in college and I really need help, thanks.

Answer 1

One thing that may be tripping you up is that the d term in the gravitational formula is the distance between the centers of mass of the objects, not the altitude above Earth's surface.

The other thing to keep an eye on is your units. The big G gravitational constant is ~6.67 x 10^-11 m³ kg^-1 s^-2; if you're using that value, make sure you're consistently using kg and m rather than grams and cm, for instance.

You're correct that the mass of an object factors out when you consider its acceleration under gravity. Near Earth's surface you should be getting acceleration of ~9.8 m s^-2.

An unpowered ballistic trajectory that can reach space from Earth's surface is going to incur a terrifying amount of drag. I know maximum dynamic pressure for a launching rocket is usually in the 20-30KPa range, but your case will be much higher, because maximum speed -- the moment of launch) -- will coincide with maximum air density, which is not the case for rockets.

Answer 2

This might help to get you started. It's just a 1D radial solver, but you can play with the math. I included a 3D version of the derivative function to show one way to make the acceleration a vector in NumPy. Don't forget to add the rotation of the earth.

The idea about starting on top of a mountain like Mt. Kilimenjaro is worth exploring for sure.

I just chose a value for the scale height, the 1/e length of the decrease in density with altitude. It's not really that simple, and of course drag will be much more complicated and much worse! than this simple equation, but it's at least a start, and you can add/modify math as you learn more about supersonic drag.

You can see the dramatic loss of velocity while still in the atmosphere.

# FROM: https://space.stackexchange.com/questions/16581/what-formulas-do-i-use-to-calculate-the-gravity-and-drag-forces-on-an-object-asc

def deriv(X, t):

    r, v = X
    acc_grav = -GMe / r**2
    acc_drag = -0.5 * rhocalc(r) * v**2 * Cd * A / mass

    return [v, acc_grav + acc_drag]

def rhocalc(r):

    altitude = r - re
    rho      = rho_0 * np.exp(-altitude / h_scale)

    return rho

def deriv3D(X, t):

    r, v = X[:3], X[3:] # 3D state vector np.array([x, y, z, vx, vy vz])

    acc_grav = -GMe * r * ((r**2).sum(axis=0))**-1.5
    acc_drag = -0.5 * rhocalc(r) * v*np.sqrt((v**2).sum()) * Cd * A / mass

    return np.hstack((v, acc_grav + acc_drag))

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

GMe = 3.986E+14   # m^3/s^2

mass = 1000       # kg
Cd   = 1.0
A    = 0.1        # m^2
rho_0   = 1.3     # kg/m^3  or  mg/cm^3
h_scale = 7600    # m   scale height 

re  = 6371000.    #m     radius of earth
alt = 0.          #m     launch altitude
v0  = 3000.       #m/s   launch velocity

X0 = [re+alt, v0]  # 1D state vector

print deriv(X0, 0)

t = np.linspace(0, 90, 1000)

tol = 1E-9

answer, blob = ODEint(deriv, X0, t, rtol=tol, atol=tol, full_output=True)

plt.figure()

r, v = answer.T

plt.subplot(4, 1, 1)
plt.plot(t, r-re)
plt.title('altitude (m)')
plt.subplot(4, 1, 2)
plt.plot(t, v)
plt.title('velocity (m/s)')
plt.subplot(4, 1, 3)
plt.plot(t, rhocalc(r))
plt.title('density (kg/m^3)')
plt.subplot(4, 1, 4)
plt.plot(t[:-1], blob['hu'])
plt.title('step size (sec)')
plt.yscale('log')
plt.show()