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

Задание:

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

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

В среде 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 - График изменения тяги двигателей на участке выведения

Задание 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 - результатов моделирования параметров орбиты выведения