Какие формулы я использую для расчета сил гравитации и лобового сопротивления объекта, поднимающегося с поверхности Земли?

Одна вещь, которая может сбить вас с толку, заключается в том, что dтермин в гравитационной формуле — это расстояние между центрами масс объектов, а не высота над поверхностью Земли.

Еще одна вещь, за которой нужно следить, — это ваши юниты. Большая гравитационная постоянная G составляет ~6,67 x 10 ^-11 м ³ кг ^-1 с ^-2 ; если вы используете это значение, убедитесь, что вы постоянно используете кг и м, а не граммы и см, например.

Вы правы в том, что масса объекта не учитывается, если учесть его ускорение под действием силы тяжести. Вблизи поверхности Земли вы должны получить ускорение ~9,8 мс ^-2 .

Баллистическая траектория без двигателя, которая может достичь космоса с поверхности Земли, будет иметь ужасающее сопротивление. Я знаю, что максимальное динамическое давление для стартующей ракеты обычно находится в пределах 20-30 кПа, но в вашем случае оно будет намного выше, т.к. максимальная скорость -- момент старта) -- будет совпадать с максимальной плотностью воздуха, что не так. для ракет.

Это может помочь вам начать работу . Это просто одномерный радиальный решатель, но вы можете играть с математикой. Я включил 3D-версию производной функции, чтобы показать один из способов сделать ускорение вектором в NumPy. Не забудьте добавить вращение Земли.

Идея о том , чтобы начать с вершины горы, такой как гора Килименжаро , безусловно , заслуживает изучения.

Я просто выбрал значение шкалы высоты , 1/е длины уменьшения плотности с высотой. На самом деле все не так просто, и конечно драг будет намного сложнее и намного хуже! чем это простое уравнение, но это по крайней мере начало, и вы можете добавлять/изменять математику, когда узнаете больше о сверхзвуковом сопротивлении.

Вы можете увидеть резкую потерю скорости еще в атмосфере.

# 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()