Решение системы дифференциальных уравнений движения

[Phil Smith]

Доброго Времени Суток!
При решении уравнений движения ракеты возникли сложности с решением системы.
Для начала привожу необходимые формулы и уравнения.

В переделах активного участка $i$-ой ступени пакеты интегрирование системы уравнений удобнее производить по аргументу $\mu_i$, который связан со временем полета ракеты следующим соотношением:

$$\mu_i = \frac {\dot m_i(t - t_{ei-1})} {m_{0i}}$$

где:
$\dot m_i$ - средний секундный расход топлива при работе двигателей $i$-ой ступени;
$t$ - текущее значение времени с момента старта (есть величина шага интегрирования);
$t_{ei-1}$ - момент выключения двигателей $(i - 1)$-й ступени. Для первой ступени $t_{ei-1} = 0$;
$m_{0i}$ - начальная масса $i$-й ступени.

Приведем уравнения движения на активном участке траектории для первой ступени.

Зададимся вспомогательными коэффициентами:

$$A_1 = \frac {g_0} {1-\mu_1}(P_{sp}^{sl} + (P_{sp}^v - P_{sp}^{sl})[1-\frac {p_h}{101360}])$$

$$A_2 = g_0\frac {70952} {P_{m1}}(\lambda_{01}P_{sp}^{sl})\frac {M^2\frac {p_h}{101360}} {1-\mu_1}$$

$$A_3 = \lambda_{01}P_{sp}^{sl}$$

где:
$g_0$ - ускорение свободного падения;
$P_{sp}^{sl}$ - Удельный импульс двигателя на уровне земли;
$P_{sp}^v$ - Удельный импульс двигателя в вакууме;
$p_h$ - плотнось воздуха на высоте $h$;
$101360$ - плотность воздуха на уровне земли [н/м2].
$P_{m1}$ - начальная нагрузка на мидель ракеты;
$M$ - число Маха. $M =V/343.2$, где $V$ - текущая скорость, $343.2$ - скорость звука [м/с];
$\lambda_{01}$ - начальная тяговооруженность первой ступени.

Тогда система уравнений движения первой ступени примет вид:

$$\begin{cases} \frac {dV_x}{d\mu_1} = A_1\cos\varphi - A_2C_x(M)\frac {Vx} V - A_2C_y^\alpha(M)\alpha\frac {Vy} V - A_3\frac {g_0R^2} {r^3}x\\ \frac {dV_y}{d\mu_1} = A_1\sin\varphi - A_2C_x(M)\frac {Vy} V + A_2C_y^\alpha(M)\alpha\frac {Vx} V - A_3\frac {g_0R^2} {r^3}(R + y)\\ \frac {dx}{d\mu_1} = A_3V_x\\ \frac {dy}{d\mu_1} = A_3V_y\\ V^2 = V_x^2 + V_y^2\\ h = r - R\\ \alpha = \varphi - \arccos\frac {Vx} V\\ \varphi = \varphi_{program}(t)\\ \end{cases} $$

где:
$V_x$ - приращение горизонтальной составляющей скорости;
$V_y$ - приращение вертикальной составляющей скорости;
$x$ - приращение дальности полета;
$y$ - приращение высоты полета;
$V$ - скорость полета;
$R$ - радиус Земного шара;
$r$ - радиус-вектор;
$h$ - высота полета;
$\alpha$ - угол атаки;
$\varphi$ - угол тангажа;
$\varphi_{program}(t)$ - программа изменения угла тангажа.
Программа изменения угла тангажа следующаяя:

$$ \varphi=\begin{cases} 10 sec \qquad 0.3833^0/sec \qquad до 49^0\\ 120 sec \qquad 0.4167^0/sec \qquad до 30^0\\ 185 sec \qquad 0.15^0/sec \qquad до 15^0\\ 295 sec \qquad 2.1667^0/sec \qquad до 0^0\\ \end{cases} $$

где первый столбец - начало времени, с которой начинается данная строка (секунды с момента старта);
второй столбец - угловая скорость изменения угла тангажа (градусов/сек);
третий столбец - предельный угол тангажа, до которого продолжается изменения в данной строке.

$C_x(M)$ и $C_y^\alpha(M)$ - аэродинамические коэффициенты. Задаются зависимостью от числа Маха (т.е. Скорости деленной на скорость звука):

$$ C_x = \begin{cases} 0.29\qquad\qquad 0 \leqslant M < 0.8\\ M - 0.51\qquad\qquad 0.8\leqslant M < 1.068\\ 0.091 + 0.5M^{-1}\qquad\qquad M\geqslant1.068\\ \end{cases} $$

$$ C_y^\alpha = \begin{cases} 2.8\qquad\qquad 0 \leqslant M < 0.25\\ 2.8 + 0.447(M - 0.25)\qquad\qquad 0.25\leqslant M < 1.1\\ 3.18 - 0.660(M - 1.1)\qquad\qquad 1.1\leqslant M < 1.6\\ 2.85 - 0.35(M - 1.6)\qquad\qquad 1.6\leqslant M < 3.6\\ 3.55\qquad\qquad M\geqslant3.6\\ \end{cases} $$

Итак, известна начальная и конечная массы ракеты, удельные импульсы в пустоте и на уровне земли,
массовый расход топлива, зависимости плотности воздуха от высоты, коэффициентов лобового и
бокового сопротивления от скорости, закон изменения угла тангажа.
Вопрос: как мне организовать решение данной системы через EXCEL (задав шаг ингегрирования допустим 0.1 сек)? Возникли сложности с численным решением данной системы. С какого преобразования стоит начать ее решать численно?
Заранее благодарен за ответ!

[Deggial]

Тема перемещена из форума «Помогите решить / разобраться (М)» в форум «Карантин»
Причина переноса: не приведены попытки решения, не совсем четкая формулировка

Приведите попытки решения, укажите конкретные затруднения, мешающие Вам решить задачу. Если Вы пытались решать задачу численно, напишите, какие численные методы Вы пытались использовать.

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

Замечу также, что формулы на форуме следует оформлять $\TeX$ом. Инструкции по оформлению формул здесь или здесь (или в этом видеоролике).

После исправлений сообщите в теме Сообщение в карантине исправлено, и тогда тема будет возвращена.

[GAA]

Тема перемещена из форума «Карантин» в форум «Помогите решить / разобраться (М)»
Причина переноса: возвращена после исправлений.

[Phil Smith]

Вот набросок мой в Excel. Указанны все исходные данные (кроме зависемости плотности атмосферы от высоты и коэффициентов аэродинамического сопротивления от скорости):
http://**invalid link**/f/853/17686082.gif/
Вышезаданный вопрос остается открытым.

[GAA]

Сразу оговорю, что о новых версиях MS Excel я ничего не знаю. Если я ничего не путаю, то в старых версиях (до 2003 включительно) не было стандартной надстройки даже для решения систем дифференциальных уравнений. Знатоки, пожалуйста, поправьте меня, если я ошибаюсь.

Если действительно требуется реализовать расчеты в MS Excel, то сделать это можно запрограммировав на VBA какой-нибудь численный метод. Если в рекомендованной литературе не указаны публикации, к которым следует обратиться, то в первую очередь нужно поискать такие работы в Интернете и в реферативных журналах, по-видимому, в механике. [Со ссылками я не помогу, на соответствующие семинары не ходил.]

Если нужно просто попробовать найти «как бы решение», то я открыл бы учебник по численным методам (на форуме была тема Методы вычислений (литература)), выбрал бы метод попроще, записал бы на форуме формулы и привел бы реализацию на VBA.

[Утундрий]

Слишком не вникал, но складывается парадоксальное впечатление. Как если бы некто озаботился переложением своих результатов в понятном для "дрессированного ученного" виде, но не учёл степени деградации последнего. К примеру, мы видим здесь явно заданную и даже заботливо кусочно аппроксимированную программу изменения направления вектора тяги, которая вообще говоря получается из решения некоего диффура в отношении так называемых сопряжённых переменных, о коих тут даже не упомянуто. В общем, перед нами классическое "подставляй и считай". И стоит ли упрекать автора медотики в том, что он (ориентируясь, видимо, на какой-то свой минимум-миниморум) не посчитал нужным объяснить заодно и аксиомы арифметики?

[Phil Smith]

Вот примерно что получаеся у меня. Приведу пример на основе перового уравнения системы:

$$\int\limits_{\mu_{i-1}}^{\mu_i} dV_x = \int\limits_{\mu_{i-1}}^{\mu_i} (A_1\cos\varphi)d\mu_i} - \int\limits_{\mu_{i-1}}^{\mu_i} (A_2C_x(M)\frac {V_x} {V}) d\mu_i - \int\limits_{\mu_{i-1}}^{\mu_i} (A_3\frac {g_0R^2}{r^3}x) d\mu_i$$

где $i$ - шаг интергирования.

Если мои изложения верны, далее прелступлю к численному решению данных интегралов, что в последстивии использую в EXCEL.