Математическое моделирование: движения ракеты-носителя, движения ракеты-носителя с регулированием кажущейся скорости

метки: Математический, Ракетостроение, Модель, Моделирование, Вектор, Объект, Матрица, Орбита

Федеральное агентство по образованию

Государственное образовательное учреждение высшего профессионального образования «Московский государственный университет леса»

Кафедра САУ

Контрольная работа

по дисциплине «Математические модели ЛА»

тема: Математическое моделирование: движения ракеты-носителя, движения ракеты-носителя с регулированием кажущейся скорости

Москва 2012 год

Математическое моделирование движения ракеты-носителя

Задание:

Разработать программу и провести расчеты в среде MatLab + Simulink по исследованию математической модели движения ракеты-носителя в плотных слоях атмосферы при выведении полезного груза на орбиту. В качестве исходных данных для моделирования принять материалы из Википедии (Интернет) по ракете-носителю «Союз-2»:

Версия

Индекс

ПН на НОО, кг

ПН на ССО, кг

ПН на ГПО, кг

Масса РН, т

ДУ 1 ст

ДУ 2 ст

ДУ 3 ст

Тип РБ

СК

Примечание

1а

14А14

Плесецк 5900-6830, Байконур 5500-7020[3]

312

РД-107А

РД-108А

РД-0110

Фрегат

Плесецк 43/4, Байконур 31

1б

14А14

Плесецк 6900-7835, Байконур 6500-8250 [3]

Плесецк 4900 [4]

312

РД-107А

РД-108А

РД-0124

Фрегат

Плесецк 43/4, Байконур 31

1в

131КС

Плесецк 2800, Байконур 2850[5]

1500 кг (c БВ Волга)[6]

157-160

НК-33, РД-0110Р

РД-0124

нет

Волга

Плесецк 43/4

СТА

372РН21

4230 [7]

2850 [7]

312

РД-107А

РД-108А

РД-0110

Фрегат

ГКЦ

СТБ

372РН21

9000-9200 [8]

4900 [7]

3240 [7]

312

РД-107А

РД-108А

РД-0124

Фрегат

ГКЦ

Союз-2 (1а, 1б, СТА, СТБ)
«Союз-2» с КА «Меридиан» в МИКе
Общие сведения
	Россия
Семейство	Р-7
Индекс	14А14 (1а, 1б), 372РН21 (СТА, СТБ)
Назначение	ракета-носитель
Разработчик	ЦСКБ-Прогресс
Изготовитель	ЦСКБ-Прогресс
Основные характеристики
Количество ступеней	3
Длина	51,1 м
Диаметр	10,3 м
Стартовая масса	313 000 кг
История запусков
Состояние	действующая
Места запуска	Байконур; Плесецк; Куру(строится)
Число запусков	9 (6 — 1а, 3 — 1б)
— успешных	9
— неудачных	0
Первый запуск	1а: 8 ноября 2004 1б: 27 декабря 2006
Последний запуск	2 ноября 2010
Первая ступень — Блоки Б, В, Г, Д
Маршевые двигатели	4 x РД-107А
Тяга	85,6 тс на уровне моря 104 тс в пустоте [9]
Удельный импульс	263,3 с на уровне моря 320,2 с в пустоте [9]
Время работы	118 с
Горючее	керосин Т-1
Окислитель	жидкий кислород
Вторая ступень — Блок А
Маршевый двигатель	РД-108А
Тяга	80,8 тс на уровне моря 94 тс в пустоте [9]
Удельный импульс	257,7 с на уровне моря 320,6 с в пустоте [9]
Время работы	285-320 с
Горючее	керосин Т-1
Окислитель	жидкий кислород
Третья ступень — Блок И
Маршевый двигатель	РД-0110(1а) или РД-0124(1б)
Тяга	30,38тс(1а) или 30,00тс(1б)
Удельный импульс	326с(1а) или 359с(1б)
Время работы	300с
Горючее	керосин Т-1 (1а) или керосин РГ-1 (1б)
Окислитель	жидкий кислород

Решение задания:

В среде MatLab + Simulink разработана программа, схема которой (mdl- file) представлена на рисунке 1.

Рисунок 1 — Схема программы моделирования движения РН в плотных слоях атмосфере (mdl- file)

Математическое моделирование движения ЛА основывается на интегрировании системы дифференциальных уравнений, описывающих законы изменения основных параметров ЛА при движении в плотных слоях атмосферы в системах координат OXYZ, O1X1Y1Z1 и O1XvYvZv. OXYZ — начальная стартовая система координат с началом координат O > на поверхности Земли в точке (0,0), где 0 и 0 — широта и долгота точки старта, плоскость OXZ > в плоскости местного горизонта, ось OX > по направлению расчетной точки выведения (точки падения, взлета и т.п.) ЛА, ось OY > по направлению местной вертикали. O1X1Y1Z1 — связанная система координат с началом координат O1 > в центре масс ЛА, ось O1X1 > по направлению продольной оси ЛА, ось O1Y1 > в поперечной плоскости по направлению оси симметрии аппарата. O1XvYvZv — скоростная система координат с осью OvXv > по направлению воздушного потока. Угловая ориентация систем координат (OXYZ) и (O1X1Y1Z1) определяется углами тангажа — , рысканья — и крена — . Угловая ориентация систем координат (O1X1Y1Z1) и (O1XvYvZv) определяется углами атаки — , скольжения — и аэродинамического угла крена — .

Переход от начальной стартовой системы координат (OXYZ) к связанной системе координат (O1X1Y1Z1) может быть осуществлён с помощью матрицы:

следующим образом: y1= А·y, где ш, х, г — углы рысканья, тангажа и крена. А=с(г)·в(х)·а(ш), где с,в,а матрицы последовательных поворотов вектора y на углы ш, х, г. y и y1 — вектора положения в системах координат OXYZ и O1X1Y1Z1, соответственно.

|cos(ш)0-sin(ш)|

а = |010 |

|sin(ш)0 cos(ш)| ,

| cos(х)sin(х)0|

b = |- sin(х)cos(х)0|

|001| ,

| 100 |

c = | 1cos(г)sin(г) |

|0 — sin(г) cos(г) | .

Уравнения движения центра масс ЛА представлены в следующем виде:

m я = AТ(T + AvRa) + G,

где я — вектор ускорения ЛА y={x,y,z}, T — вектор тяги маршевых двигателей ЛА, Ra={Xa,Ya,Za} — вектор аэродинамических сил в поточной с.к.( Xa=Сxa·Sm·q, Ya=Сya·Sm·q, Za=Сza·Sm·q, q=сV2/2 — скоростной напор, с — плотность, Sm — площадь миделевого сечения, Сxa, Сya, Сza — аэродинамические коэффициенты, V — модуль вектора воздушной скорости), G — вектор силы тяготения, AТ — матрица перехода от связанной системы координат (с.к.) к начальной стартовой с.к., определяемая углами , и , Av — матрица перехода от поточной с.к. к связанной с.к., определяемая углами и , m = m0 — dm t (m0 — начальная масса ЛА, dm — массовый секундный расход, t — текущее время полета).

Вектор воздушной скорости ЛА, используемый при определении углов и , вычисляется следующим образом Vv =V + W. V — вектор скорости ЛА, определяемый, как и текущий вектор положения ЛА — r, в результате интегрирования дифференциальных уравнений линейного движения ЛА, W — текущий вектор скорости ветра.

Вектор угловой скорости объекта можно представить в виде суммы следующих слагаемых: щ. Проектируя это уравнение на оси связанной системы координат, получим кинематические соотношения:

Уравнения вращательного движения:

J·ю=Ma +Mu,

где Mu и Ma моменты сил управления и аэродинамических сил, определяемых в т.ч. текущими отклонениями центров масс и давления ЛА. В расчётах принята экспоненциальная модель атмосферы: = 0е-h, где — градиент изменения плотности по высоте полёта ЛА- h.

Для заданной программы движения ЛА в плотных слоях атмосферы на рисунке 2 представлен график изменения тяги двигателей на участке выведения.

Рисунок 2 — График изменения тяги двигателей на участке выведения

Математическое моделирование параметров орбиты выведения

Задание:

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

Решение задания:

Параметры орбиты выведения определяются терминальными параметрами движения ракеты-носителя, определяемыми следующими соотношениями:

(Представить основные уравнения эллиптического движения и зависимости элементов орбиты с начальными условиями — терминальными (конечными) параметрами выведения ракеты-носителя).

В среде MatLab + Simulink разработана блок программы, m- file которой представлен на рисунке 3.

function [f]=fnomin2(q);

• Rz=6371000;muZ=3.98602*10^14;

• l=0.0;fi=1.0;az=0.0;

• q0(1,1)=q(1);q0(2,1)=q(2)+Rz;q0(3,1)=q(3);

• q0v(1,1)=q(4);q0v(2,1)=q(5);q0v(3,1)=q(6);

• u(1,1)=0;u(1,2)=0;u(1,3)=-1;

• u(2,1)=0;u(2,2)=1;u(2,3)=0;

• u(3,1)=1;u(3,2)=0;u(3,3)=0;

• a(1,1)=1;a(1,2)=0;a(1,3)=0;

• a(2,1)=0;a(2,2)=sin(1);a(2,3)=cos(1);

• a(3,1)=0;a(3,2)=-cos(1);a(3,3)=sin(1);

• b(1,1)=cos(fi);b(1,2)=sin(fi);b(1,3)=0;

• b(2,1)=-sin(fi);b(2,2)=cos(fi);b(2,3)=0;

• b(3,1)=0;b(3,2)=0;b(3,3)=1;

• c(1,1)=cos(az);c(1,2)=0;c(1,3)=-sin(az);

• c(2,1)=0;c(2,2)=1;c(2,3)=0;

• c(3,1)=sin(az);c(3,2)=0;c(3,3)=cos(az);

• uabc=u*a*b*c;qg=uabc*q0;
• qgv=uabc*q0v;

• rg=sqrt(qg(1,1)^2+qg(2,1)^2+qg(3,1)^2);vg=sqrt(qgv(1,1)^2+qgv(2,1)^2+qgv(3,1)^2);vkr=sqrt(muZ/rg);

• alf=acos((qg(1,1)*qgv(1,1)+qg(2,1)*qgv(2,1)+qg(3,1)*qgv(3,1))/rg/vg);

• c1=qg(2,1)*qgv(3,1)-qg(3,1)*qgv(2,1);c2=qg(3,1)*qgv(1,1)-qg(1,1)*qgv(3,1);c3=qg(1,1)*qgv(2,1)-qg(2,1)*qgv(1,1);

• f1=-muZ*qg(1,1)/rg+c3*qgv(2,1)-c2*qgv(3,1);f2=-muZ*qg(2,1)/rg+c1*qgv(3,1)-c3*qgv(1,1);f3=-muZ*qg(3,1)/rg+c2*qgv(1,1)-c1*qgv(2,1);

• cs=sqrt(c1^2+c2^2+c3^2);fs=sqrt(f1^2+f2^2+f3^2);

• p=cs^2/muZ;e=fs/muZ;ra=p;rp=p;

• if (e<1) ra=p/(1-e);rp=p/(1+e);end;

• ap=(rp+ra)/2;
• Tz=2*pi*sqrt(ap^3/muZ);
• in=acos(c3/cs);
• om=0;if in>0 kom=-c2/(cs*sin(in));som=c1/(cs*sin(in));om=acos(kom);end;

• lr=qg(1,1)/rg;mr=qg(2,1)/rg;nr=qg(3,1)/rg;lv=qgv(1,1)/vg;mv=qgv(2,1)/vg;nv=qgv(3,1)/vg;

• f(1)=ra;
• f(2)=rp;
• f(3)=p;
• f(4)=e;
• f(5)=in;
• f(6)=om;
• f(7)=vkr;

Таблица 1 — результатов моделирования параметров орбиты выведения

Параметры движения в конце участка выведения	Орбитальные параметры для условий старта: по широте = 60, долготе = 60, азимуту плоскости выведения =0
X0=847600	ra=7042000
Y0=226900	rp=2349000
Z0=-1280000	e=0.4997
Vx0=3226	i=2.541
Vy0=-288.9	=0.843
Vzo=-4723	Tзв=3202

Математическое моделирование движения ракеты-носителя с регулированием кажущейся скорости

ракета движение орбита скорость

Задание:

Доработать, программу моделирования движения ракеты-носителя программные блоки для решения задачи регулирования кажущейся скорости (РКС).

Представить в табличном или графическом виде данные по отклонениям параметров траектории выведения от номинальной орбиты выведения для системы управления с РКС для различных значениях возмущающих факторов (разбросы стартовой массы, удельной тяги, суммарного расхода рабочего тела и др.).

Решение задания:

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

Решение задач выведения ракеты-носителя осуществляется на принципах:

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

Системы управления, построенные по методу «жестких» траекторий отличаются простыми алгоритмами и аппаратурной реализацией с использованием трёх независимых каналов: управления боковым, нормальным и продольным движением, а также канала выключения двигательных установок. В основе принципа действия подсистемы РКС лежит управление тягой маршевых двигателей в зависимости от отклонений параметров продольного движения ракеты-носителя от их номинальных значений: д = — К·?W, где К — коэффициент регулятора, ·?W — приращение кажущейся скорости.

Канал РКС предназначен для выдерживания программного значения продольной составляющей кажущейся скорости путём формирования или дросселирования основной ДУ — с целью уменьшения разброса параметров движения центра масс ракеты, что необходимо для упрощения алгоритма выключения ДУ. (РКС может отсутствовать, если можно обеспечить требования по точности без РКС).

Исполнительным органом м. б. привод винта редуктора двигателя для изменения секундного расхода топлива в камерах сгорания, что приводит к изменению давления в камерах сгорания и, следовательно, тяги ДУ. Тяга возрастает, если фактическая скорость меньше программной — замкнутая автоматическая система РКС. Для улучшения динамических свойств возможно введение внутренней обратной связи по давлению в камерах сгорания.

Рисунок 3

Рисунок 4 — mdl- file блока программы моделирования движения ракеты-носителя c решением задачи регулирования кажущейся скорости (РКС).

Математическое моделирование линейных систем

Исследования и реализация задач управления бортовыми системами управления на основе указанных методов осуществляется в большинстве случаев с использованием линеаризованных уравнений возмущённого движения РН на участке выведения.

В этих уравнениях n- вектор x (матрица-столбец nЧ1) имеет своими составляющими переменные состояния x1,x2,…,xn и называется вектором состояния объекта, r — вектор y, составляющие которого — выходные сигналы объекта, называется выходным вектором, а m- вектор u представляет собой воздействия, которые могут быть поданы на объект со стороны. Матрица объекта A, матрица управления B и матрица выходного сигнала C имеют соответственно размеры nЧ n, nЧm, rЧn.

Условия управляемости линейных стационарных систем зависят от ранга матрицы управляемости Qy = [B,AB,A2B].

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

Полная наблюдаемость системы

означает возможность определения начального состояния х0 этой системы по её выходному сигналу у(t), известному на некотором конечном интервалу t (началу интервала соответствует х0).

Если для этой системы ранг матрицы наблюдаемости

Qн = [CT|АT CT| (АT)2 CT|…(АT)n-1 CT]

равен порядку n системы, то система полностью наблюдаема.

Если ранг Qн меньше порядка объекта, то это означает, что по измеренному выходному сигналу у(t) можно оценить не все, а лишь некоторую часть переменных состояния объекта.

Задание:

1. объект:

(колебания ЛА относительно оси рысканья и уравнение рулевой машины).

2. Объект:

Решение заданий:

1. Приняв для переменных состояния объекта обозначения

можно записать систему в виде одного матричного уравнения

Приняв для элементов матриц числовые значения, удобные при расчётах

Матрица управляемости объекта

Qy = [B,AB,A2B] = имеет ранг, равный порядку объекта.

2. Определитель матрицы Qн = в2. Следовательно, ранг Qн равен порядку объекта n = 4. Т.о. измеряя y = x3 , можно при помощи наблюдающего устройства восстановить полный вектор состояния объекта [x1, x2, x3, x4], т.е. объект полностью наблюдаем.

Динамические и частотные характеристики САУ

Задана передаточная функция САУ

w=(2*p^2-3*p)/(9*p^3+3*p^2+2*p+4)

Найдем ее динамические и частотные характеристики.

1. Создадим LTI-объект с именем w, для этого выполним:

• >>
• p=tf(‘p’)

>> w=(2*p^2-3*p)/(9*p^3+3*p^2+2*p+4)

Transfer function:

2 p^2 — 3 p

9 p^3 + 3 p^2 + 2 p + 4

2. Найдем полюса и нули передаточной функции с использованием команд pole, zero.

3. Построим переходную функцию командой step(w).

Результат ее выполнения приведен на рис. 5.

Рисунок 5 — Переходная функция h(t)

4. Построим импульсную переходную функцию командой impulse(w).

Результат показан на рис. 6.

Рисунок 6 — Импульсная переходная функция

5. Диаграмму Боде получим, используя команду bode(w) — рис. 7.

Рисунок 7 — Логарифмические частотные характеристики

6. Определим частотный годограф Найквиста, выполнив команду nyquist(w) — рис. 8.

Рисунок 8 — Частотный годограф

Список использованных источников

[Электронный ресурс]//URL: https://drprom.ru/kontrolnaya/matematicheskaya-model-raketostroeniya/

1. Ю.Г. Сихарулидзе Баллистика летательных аппаратов. Москва, “Наука”, 1982.

2. Л.Н. Лысенко Наведение и навигация баллистических ракет. Издательство МГТУ им. Н.Э. Баумана, 2007 г.

3. А.П. Разыграев Основы управления полётом космических аппаратов и кораблей.