Продраться через атмосферу.

Вопрос в том, какой режим оптимален для прохождения атмосферы. Чтоб вытащить в вакуум как можно больше полезной нагрузки.

Обозначения:
m — текущая масса ракеты;
m0 — конечная масса;
v — текущая вертикальная скорость;
v0 — скорость продуктов сгорания;
rho — плотность воздуха;
Cx — коэффициент аэродинамического сопротивления;
S — площадь сечения ракеты;
g — 9.81.
h — высота (путь).

Первое приближение: Cx не зависит от скорости, плотность атмосферы постоянная и скорость подъёма постоянная.

Вот какие силы действуют на ракету (кроме тяги):

Fc = - Cx * S * rho * v² / 2 — сопротивление

Fg = - m * g — вес

Предположим, что v0 очень большая, так что массу можно считать постоянной.

Тогда импульс, расходуемый на прохождение пути dh будет:

dp = dh / v * (Fc + Fg)

dp / dh = (Fc + Fg) / v = Cx * rho * v / 2 + m * g / v

Чем меньше израсходуем тяги на прохождение единичного пути, тем лучше.
Ищем минимум последнего выражения, для этого берём производную:

d/dv(dp / dh) = Cx * rho / 2 - m * g / v²

Ноль производной, и, соответственно искомый минимум получается при:

Cx * rho / 2 = m * g / v²

Домножим это на квадрат скорости:

Cx * rho * v² / 2 = m * g

Получается, что оптимальна скорость подъёма, это когда сила тяжести равна силе сопротивления.
А тяга, соответственно, равна удвоенному весу ракеты.

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

Так как плотность атмосферы и масса ракеты меняются, соответственно надо подкручивать тягу, чтоб от оптимума далеко не уходить, а на выходе из атмосферы надо будет резко добавлять газу, чтоб угнаться за оптимальной скоростью.

---

Симуляция.

Рассматривается ракета с удлинением 20.
Её сечение считается по формуле (всё в СИ):

D = (0.001 * M / 20)^(1/3)

S = pi / 4 * d² = pi / 4 * (0.001 * M / 20)^(2/3) = 0.001066 * M^(2/3)

Это соответствует средней плотности ракеты 1.25 г/см3, примерно в полтора раза меньше плотности топлива.

Если, например, массу взять 10 кг, то сечение получится 0.005 м2 = 50 см2 (80 мм в диаметре).

Ессно, M в формуле стартовая масса, а не текущая!

Плотность атмосферы берётся из стандартной таблицы и интерполируется.
Но можно взять и просто нечто среднее по простой формуле:

rho = exp(-h / 7000)

Ну, попробуем взять тягу просто равную утроенному весу. Это недалеко от оптимальной (удвоенной), так как ракете надо ещё и разгоняться.

Ft = 3 * m * g

Расход массы будет:

dm /dt = Ft / v0

Осталось написать симуляторную прогу:

m = 10
S = 0.001066 * m^(2/3)
Cx = 0.45
t = 0
h = 0
v = 0
g = 9.8
v0 = 2000 — удельный импульс (Н/кг)
dt = 0.01 — шаг по времени

while (v < 1000)
{
rho = exp(-h / 7000)
Fc = - Cx * S * rho * v² / 2
Fg = - m * g
Ft = 3 * m * g
a = (Fc + Fg + Ft) / m
dv = a * dt
dm = -Ft / v0 * dt
dh = v * dt
v = v + dv
m = m + dm
h = h + dh
t = t + dt
}

Ну вот у меня получается, что на 76 секунде от 10 кг останется 3.3, скорость будет 1 км/с, высота 37 км.

Возражайте.
Где я нулики потерял?
Чего не так сделал?

А откуда формула rho = exp(-h / 7000) ?
По твоей ссылке, в таблице, плотность на 50 км 1.15*10-3 кг/м3, а по формуле 0.79*10-3. Да и при h = 0 плотность не еденица Там наверно коэффициент забыл?

Предлагаю вычислять так:
//http://www.grc.nasa.gov/WWW/K-12/airplane/atmosmet.html
Получает данные как в той таблице, проверил


double ThermCels(double h)
{
if (h > 25000.0) return -131.21 + 0.00299 * h;
else if (h > 11000.0) return -56.46;
else return 15.04 - 0.00649 * h;
}

double PressKPa(double h)
{
if (h > 25000.0) return 2.488 * pow((ThermCels(h) + 273.1) / 216.6, -11.388);
else if (h > 11000.0) return 22.65 * pow(2.7183, (1.73 - 0.000157 * h));
else return 101.29 * pow((ThermCels(h) + 273.1) / 288.08, 5.256);
}

double Rho(double h)
{
return PressKPa(h) / (0.2869 * (ThermCels(h) + 273.1));
}